# Sea Ice Rheology Experiment (SIREx): 1. Scaling and Statistical Properties of Sea-Ice Deformation Fields

Amélie Bouchat, Nils Hutter, Jérôme Chanut, Frédéric Dupont, Dmitry Dukhovskoy, Gilles Garric, Younjoo J. Lee, Jean-François Lemieux, Camille Lique, Martin Losch, et al.

► To cite this version:

Amélie Bouchat, Nils Hutter, Jérôme Chanut, Frédéric Dupont, Dmitry Dukhovskoy, et al. Sea Ice Rheology Experiment (SIREx): 1. Scaling and Statistical Properties of Sea-Ice Deformation Fields. Journal of Geophysical Research. Oceans, 2022, 127, 10.1029/2021JC017667. insu-03683297

HAL Id: insu-03683297

https://insu.hal.science/insu-03683297

Submitted on 31 May 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

10.1029/2021JC017667

This article is a companion to Hutter et al. (2022a), https://doi.org/10.1029/2021JC017666.

**Key Points:**

* Power law scaling and multi-fractality of deformations in space and time can be achieved by both plastic and brittle sea-ice rheologies

* Scaling statistics of simulated sea-ice deformation fields depend on the model configuration and physical parameterizations

* Finite-difference plastic models need to be run at higher resolution than observations to agree with the observed deformation statistics

**Correspondence to:**

A. Bouchat,
amelie.bouchat@mail.mcgill.ca

**Citation:**

Bouchat, A., Hutter, N., Chanut, J., Dupont, F., Dukhovskoy, D., Garric, G., et al. (2022). Sea Ice Rheology Experiment (SIREX): 1. Scaling and statistical properties of sea-ice deformation fields. *Journal of Geophysical Research: Oceans*, 127, e2021JC017667. https://doi.org/10.1029/2021JC017667

Received 15 JUN 2021
Accepted 4 MAR 2022

**Author Contributions:**

**Conceptualization:** Amélie Bouchat, Nils Hutter

**Data curation:** Amélie Bouchat

**Investigation:** Amélie Bouchat

**Methodology:** Amélie Bouchat, Nils Hutter

**Project Administration:** Amélie Bouchat, Nils Hutter

**Resources:** Amélie Bouchat, Nils Hutter, Jérôme Chanut, Frédéric Dupont, Dmitry Dukhovskoy, Gilles Garric, Younjoo J. Lee, Jean-François Lemieux, Camille Lique, Martin Losch, Wieslaw Maslowski, Paul G. Myers, Einar Ólason, Pierre Rampal, Till Rasmussen, Claude Talandier, Bruno Tremblay, Qiang Wang

© 2022. The Authors.

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

Sea Ice Rheology Experiment (SIREx): 1. Scaling and Statistical Properties of Sea-Ice Deformation Fields

<figure><img src="image_2.png" /></figure>

Amélie Bouchat¹, Nils Hutter², Jérôme Chanut³, Frédéric Dupont⁴, Dmitry Dukhovskoy⁵, Gilles Garric³, Younjoo J. Lee⁶, Jean-François Lemieux⁷, Camille Lique⁸, Martin Losch², Wieslaw Maslowski⁶, Paul G. Myers⁹, Einar Ólason¹⁰, Pierre Rampal¹¹, Till Rasmussen¹², Claude Talandier⁸, Bruno Tremblay¹, and Qiang Wang²

¹Department of Atmospheric and Oceanic Sciences, McGill University, Montréal, QC, Canada, ²Alfred-Wegener-Institut, Helmholtz-Zentrum für Polar- und Meeresforschung, Bremerhaven, Germany, ³Mercator Ocean International, Ramonville-Saint-Agne, France, ⁴Service Météorologique Canadien, Environnement et Changement Climatique Canada, Dorval, QC, Canada, ⁵Center for Ocean-Atmospheric Prediction Studies, Florida State University, Tallahassee, FL, USA, ⁶Department of Oceanography, Naval Postgraduate School, Monterey, CA, USA, ⁷Recherche en Prévision Numérique Environnementale, Environnement et Changement Climatique Canada, Dorval, QC, Canada, ⁸Laboratoire d'Océanographie Physique et Spatiale (LOPS), University of Brest, CNRS, IRD, Ifremer, IUEM, Brest, France, ⁹Department of Earth and Atmospheric Sciences, University of Alberta, Edmonton, AB, Canada, ¹⁰Nansen Environmental and Remote Sensing Centre, and Bjerknes Centre for Climate Research, Bergen, Norway, ¹¹Institut de Géophysique de l'Environnement, CNRS, Grenoble, France, ¹²Danish Meteorological Institute, Copenhagen, Denmark

**Abstract** As the sea-ice modeling community is shifting to advanced numerical frameworks, developing new sea-ice rheologies, and increasing model spatial resolution, ubiquitous deformation features in the Arctic sea ice are now being resolved by sea-ice models. Initiated at the Forum for Arctic Modeling and Observational Synthesis, the Sea Ice Rheology Experiment (SIREX) aims at evaluating state-of-the-art sea-ice models using existing and new metrics to understand how the simulated deformation fields are affected by different representations of sea-ice physics (rheology) and by model configuration. Part 1 of the SIREX analysis is concerned with evaluation of the statistical distribution and scaling properties of sea-ice deformation fields from 35 different simulations against those from the RADARSAT Geophysical Processor System (RGPS). For the first time, the viscous-plastic (and the elastic-viscous-plastic variant), elastic-anisotropic-plastic, and Maxwell-elasto-brittle rheologies are compared in a single study. We find that both plastic and brittle sea-ice rheologies have the potential to reproduce the observed RGPS deformation statistics, including multi-fractality. Model configuration (e.g., numerical convergence, atmospheric representation, spatial resolution) and physical parameterizations (e.g., ice strength parameters and ice thickness distribution) both have effects as important as the choice of sea-ice rheology on the deformation statistics. It is therefore not straightforward to attribute model performance to a specific rheological framework using current deformation metrics. In light of these results, we further evaluate the statistical properties of simulated Linear Kinematic Features in a SIREX Part 2 companion paper.

**Plain Language Summary** The ice in the Arctic Ocean is not continuous: it is broken into individual pieces of ice (floes). As the winds and ocean currents continually move these ice floes, they get piled up together or pushed away from each other, forming regions of increased ice thickness (ridges) or regions of open water (leads). These leads and ridges (ice deformations) are important features of the Arctic pack ice because they control the amount of energy that can be exchanged between the atmosphere and the ocean. Current climate models cannot simulate individual ice floes and their deformations. Instead, various methods are used to represent the movement and deformation of the Arctic sea-ice cover. The goal of the Sea Ice Rheology Experiment (SIREX) is to compare these different methods and evaluate the ability of a large number of sea-ice models to reproduce observed sea-ice deformations from satellite imagery. SIREX is divided in two parts. In Part 1 (this study), we evaluate how the intensity of ice deformations varies in space and time. In Part 2 (companion paper), we track and evaluate the occurrence of specific deformation features. With this work, we show how to improve sea-ice models for realistic simulations of sea-ice deformations.
**Writing – original draft:** Amélie Bouchat

**Writing – review & editing:** Amélie Bouchat, Nils Hutter, Jérôme Chanut, Frédéric Dupont, Dmitry Dukhovskoy, Gilles Garric, Younjoo J. Lee, Jean-François Lemieux, Camille Lique, Martin Losch, Wieslaw Maslowski, Paul G. Myers, Einar Ólason, Pierre Rampal, Till Rasmussen, Claude Talandier, Bruno Tremblay, Qiang Wang

# 1. Introduction

Statistical properties of small-scale sea-ice dynamics derived from buoy records and synthetic aperture radar (SAR) imagery in the Arctic Ocean have been extensively documented in the last two decades. Observations from the RADARSAT Geophysical Processor System (RGPS) show that deformations in shear and divergence (positive and negative) define highly localized Linear Kinematic Features (LKFs—e.g., Kwok, 2001) and complex scaling laws describe their localization over a wide range of spatial and temporal scales (Marsan et al., 2004; Marsan & Weiss, 2010; Rampal et al., 2008; Stern & Lindsay, 2009). Specifically, the mean total deformation rates follow a power law with increasing spatial and temporal scales and the scaling exponent of this power law increases nonlinearly when considering higher moments of the deformation distribution, suggesting that very large deformation rates significantly affect the mean deformation statistics (Rampal et al., 2019; Weiss & Dansereau, 2017). These properties are reminiscent of fully turbulent flows, which also exhibit strong heterogeneity and intermittency and are characterized as multi-fractal processes (e.g., Benzi et al., 1984; Schmitt et al., 1994). As such, the observed sea-ice deformation characteristics might provide meaningful information about the underlying mechanisms governing the sea-ice mechanics. For example, the highly localized LKFs have been hypothesized to result from brittle compressive shear faulting (Schulson, 2004), while the sea-ice deformation multi-fractality and scaling laws are sometimes associated with the presence of a threshold/trigger, stress relaxation, and damage/healing mechanisms (Dansereau et al., 2016; Marsan & Weiss, 2010; Weiss & Dansereau, 2017).

In sea-ice dynamical models, a rheology describes the relation between the applied load and resulting deformation, effectively representing the sea-ice mechanical response to the external forcing. The viscous-plastic (VP) rheology with elliptical yield curve (Hibler, 1979) and its elastic-viscous-plastic (EVP) variant (Hunke & Dukowicz, 1997, 2002) are the most widely used in regional and Global Climate Models (see e.g., Stroeve et al., 2014). In the standard VP rheology, sea ice is assumed to deform as a plastic material when the mechanical stresses reach prescribed critical loads in compression, shear, and tension (as defined by the elliptical yield curve), and as a creeping, highly viscous fluid for smaller stresses. The EVP variant assumes the same physical concepts but uses damped artificial elastic waves that allow for an explicit numerical implementation of the dynamical equations. In this sense, the EVP approach can be considered as an alternative numerical solver for the VP rheology. Since its formulation, extensive work has been done on improving the speed and stability of the numerical schemes used for solving the (E)VP equations (e.g., Bouillon et al., 2013; Kimmritz et al., 2016; Lemieux et al., 2010, 2008), but parallel work has also pointed out inconsistencies in its basic physical assumptions (e.g., Coon et al., 2007). This has led to reconsideration of the classical (E)VP rheology by, among others, adding tensile strength (König Beatty & Holland, 2010; Zhang & Rothrock, 2005) and developing sea-ice rheologies based on different physical assumptions. Of these, the elastic-plastic-anisotropic (EAP)—Tsamados et al., 2013; Wilchinsky & Feltham, 2006) builds upon the artificial elastic closure of the EVP approach, but represents anisotropy of the ice stress by parameterizing the interactions of diamond-shaped floes. Long-range elastic interactions have also been explicitly included in the elasto-brittle (EB) and Maxwell-elasto-brittle (MEB) rheologies, in which the classical plastic response of the ice was traded in favor of a brittle parameterization accounting for fracturing and sliding of ice along fault planes (Bouillon & Rampal, 2015b; Dansereau et al., 2016; Girard et al., 2011).

Sea-ice models (and sea-ice rheologies) have traditionally been evaluated by estimating the error between the simulated and observed large-scale features such as sea-ice drift, thickness, concentration, and extent (e.g., Flato & Hibler, 1992; Ip et al., 1991; Kreyscher et al., 2000; Ungermann et al., 2017; Zhang & Rothrock, 2005). Given that these large-scale error metrics can generally be minimized by tuning the model thermodynamics, the sea-ice modeling community has recently introduced additional metrics that specifically evaluate the small-scale deformation statistics with the goal of better discriminating/calibrating the different sea-ice rheologies. Rheology and deformation metrics are of particular interest for modeling applications requiring accurate small-scale deformation statistics (e.g., short-term drift forecasting for navigation), but also potentially for climate projections since sea-ice deformations affect ice production, vertical heat and moisture fluxes, and salt rejection to the surface ocean. Using the observed sea-ice deformation statistics has now become common practice to validate or constrain sea-ice rheologies (e.g., Bouchat & Tremblay, 2017; Bouillon & Rampal, 2015b; Girard et al., 2009; Hutter et al., 2018; Sreen et al., 2017). Specifically, the observed strain rate probability density functions (PDFs) decay exponent and the spatio-temporal scaling exponents of the total deformation rates are used as metrics to assess the ability of sea-ice rheologies and models to reproduce large deformation events and their localization and multi-fractal properties.
The application of these deformation metrics resulted in a debate about the ability of the VP sea-ice rheology to reproduce the observed deformation statistics, justifying the need for the new EB/MEB rheology (Girard et al., 2011, 2009; Rampal et al., 2016). It has since been shown that the VP rheology is able to reproduce similar deformation characteristics as the EB/MEB rheology based on the same deformation metrics (Bouchat & Tremblay, 2017; Hutter & Losch, 2020; Hutter et al., 2018; Sreen et al., 2017), leaving open the question as to whether those metrics can be used to robustly discriminate between sea-ice rheologies, and if the deformation metrics accurately capture differences in the underlying deformation statistics. Additionally, deformation fields of other sea-ice rheologies (e.g., EAP) have not been thoroughly evaluated using the new set of deformation metrics as for VP and MEB rheologies. A comprehensive assessment of the ability of different sea-ice models and rheologies to reproduce the observed deformation statistics and the sensitivity of the deformation metrics to model parameterizations was therefore identified by the sea-ice modeling working group at the Forum for Arctic Modeling and Observational Synthesis Annual Meeting 2017 as a priority for the sea-ice modeling community.

To this end, the Sea Ice Rheology Experiment (SIREx) model intercomparison project was devised with the goal of evaluating if the sea-ice deformation metrics, as currently applied, are useful to discriminate between the different sea-ice models/rheologies and formulate more appropriate metrics if found necessary, and (b) determining how different physical parameterizations and model configuration can impact the simulated deformation and LKFs statistics to formulate recommendations for improving the representation of sea-ice deformations in future model developments. SIREx takes the form of a diagnostic model intercomparison project in which participating models are not constrained by the same configuration. This practical framework allows for low-level entry requirements, ensuring that a large number of modeling groups participate in the study. Most of the runs collected and analyzed here are sea-ice simulations that were in fact already performed for other studies, rather than mandated for SIREx specifically. This allows us to apply the deformation metrics to a wide range of sea-ice simulations to better determine their usefulness in a broader context, but it also limits the extent of what conclusions can and cannot be firmly assessed by comparing such a large number of unconstrained simulations with limited output information. For example, given that only two winters of simulations are analyzed in SIREx, our investigation cannot distinguish if disagreement between models and observations (or in between models) originates simply from the models internal variability or if it truly represents shortcomings of the analyzed simulations. In this sense, the scope of SIREx as a model intercomparison project is to identify and quantify commonalities in the ability of a large number of sea-ice models to reproduce the observed sea-ice deformation statistics in the Arctic, but it cannot extend to explain all specific drivers of the underlying mechanisms governing the observed and simulated sea-ice dynamics. This would require extensive model output that was not asked of participating simulations (see e.g., the output priorities that were identified by the Sea Ice Model Intercomparison Project for identification of internal variability drivers—Notz et al., 2016). The work performed in SIREx is nonetheless a necessary first step in order to further our understanding of how to evaluate and compare simulated sea-ice deformation fields, which helps formulate specific questions for future inter-model rheological comparison work with informed constraints.

The SIREx analysis is divided in two parts. First, and the subject of the present publication, the statistical distributions (PDFs) and the spatio-temporal scaling properties of the deformation fields are analyzed. Second, a feature-based comparison of the sea-ice deformation fields is performed using a recent LKF detection and tracking algorithm (Hutter et al., 2019) and is presented in a companion SIREx publication (Hutter et al., 2022a). In both parts, simulated sea-ice deformations for two winters (1997 and 2008) are compared to the full RGPS observational records, and all sea-ice rheologies commonly used in current sea-ice models (i.e., (E)VP, EAP, and MEB) are compared for the first time within the same framework.

The present article is organized as follows. The model specifications and observations used in this study are presented in Section 2. The methods used to obtain the simulated and observed deformation fields, as well as the deformation statistics and metrics used for comparison are detailed in Section 3. Results are presented in Section 4, followed by a discussion and recommendations for model development in Section 5. Finally, a summary and concluding remarks are presented in Section 6.
## 2. Models and Observations

A total of 35 simulations from 11 different models contributed to SIREx Part 1. The participating models were not constrained to use the same atmospheric/oceanic representation (e.g., reanalysis with different spatial and temporal resolution vs. interactively coupled atmospheric/oceanic models), and they also vary in their spatial and temporal resolution, grid type (e.g., Eulerian vs. Lagrangian), physical parameterizations (e.g., ice strength parameters and ice thickness distribution), numerical convergence criterion, etc. A list of all simulations and key sensitivity parameters are given in Table 1 (see also Table A1 in Appendix A for a list of model, configuration, and reanalysis acronyms). Note that the MERCATOR model participates only to SIREx Part 1, while other runs from Part 2 are not analyzed here (i.e., FESOM2 model, and the MITgcm 2 km, ITD—1997 run). For more information about the models, the reader can refer to the respective references in Table 1.

Daily sea-ice velocity, thickness, and concentration fields for January-February-March of 1997 and 2008 were requested from all participating models. Only two years (1997 and 2008) are considered in this study to allow for low-level entry participation, while at the same time sampling the available observational records over an extended period of time during which an accelerated sea-ice retreat and changing mechanical conditions could have led to different ice dynamical conditions (e.g., accelerated drift speeds and larger deformation rates—Rampal et al., 2009). Note however that the observed inter-annual variability in RGPS scaling statistics during this period is small (Bouchat & Tremblay, 2020). The sea-ice velocity fields are used to compute the deformation estimates (see below and Section 3.2), while the ice thickness and concentration fields were investigated to evaluate the possible impacts of simulating different ice conditions on the deformation statistics (not shown). No clear correspondence was found between models showing good agreement with RGPS in terms of deformation statistics and models showing smaller bias with respect to sea-ice extent observations and with respect to the Arctic sea-ice volume estimated by the Pan-Arctic Ice Ocean Modeling and Assimilation System (Zhang & Rothrock, 2003). Additionally, conclusions of the evaluation of simulated deformation statistics were noted to be independent of initialization with different ice volume conditions (Bouchat & Tremblay, 2017). Details of the simulated sea-ice thickness and concentration fields are therefore not expected to play a major role in the comparison of the simulated deformation statistics presented here. In the following, we specifically analyze the effects of sea-ice rheology jointly with spatial resolution (Section 4.1), ice strength (Section 4.2.1), ice thickness distribution parameterization (Section 4.2.2), and atmospheric influence (Section 4.2.3) on the deformation statistics. As discussed in Section 5, we do not address the possible oceanic influence on the simulated deformation fields, and this should be evaluated in a dedicated experiment.

All participating models provided daily output on an Eulerian grid, except for neXtSIM which uses a Lagrangian grid and provided output directly as Lagrangian trajectories. While spatial scaling can be studied using either Eulerian or Lagrangian deformation fields, temporal scaling requires the deformation history of tracked elements and therefore needs to be performed in a Lagrangian framework. We, therefore, construct offline Lagrangian trajectories from the daily Eulerian model output, which are then used to compute the deformation statistics (see details in Section 3.2). The choice of daily output for performing the Lagrangian trajectory integration from Eulerian model output was done for practical reasons: we did not ask new runs to be made for SIREx, but rather collected simulations that already had output at this frequency. This implies that sub-daily processes influencing the simulated sea-ice drift are not resolved in our Lagrangian trajectory integration for Eulerian models. We have however verified, using a subset of the participating models, that the Lagrangian deformation statistics presented below are robust to using higher-frequency output for integrating the Lagrangian trajectories (i.e., hourly means vs. daily means) for both low- and high-resolution model runs (not shown). Note also that most model output were provided as daily means, but some groups provided daily snapshots at an instant in time such as noon. The deformation statistics are also robust to the choice of temporal averaging of the model output (i.e., snapshots or daily means—not shown).

The simulated deformation statistics are compared with those derived from the RGPS Lagrangian motion data set. The RGPS Lagrangian motion data set is given as a list of trajectories (time and position) for a 10 × 10 km grid that is initialized at the beginning of the winter season over the central Arctic Ocean for different satellite passes (i.e., streams), tracked using sequential SAR images (Kwok, 1998). The nominal spatio-temporal resolution of the RGPS Lagrangian data set is $T^* = 3$ days and $L^* = 10$ km, however sampling of the RGPS Lagrangian data set is non-uniform given that trajectories are not always updated on the same days or at the same times or can be missing if the tracking on the SAR images was unsuccessful. For this reason, a pre-processing of
<table><caption>**Table 1**
Key Parameters of Model Runs Participating in SIReX Part 1</caption><thead><tr><th>Model/Configuration (Group)</th><th>Year</th><th>Grid spacing, Time step</th><th>Grid</th><th>Rheology (solver/nt, it)</th><th>Ice Strength Parameters</th><th>ITD #</th><th>Atm. Model or Reanalysis (Δx, Δt)</th><th>Reference</th></tr></thead><tbody><tr><td colspan="9"><strong>MITTgcm (Alfred-Wegener-Institute - AWI)</strong></td></tr><tr><td>MITTgcm (2 km, ITD)</td><td>2008</td><td>2 km, 120 s</td><td>E</td><td>VP (LSR)</td><td>P*, e = 22.64, 2.0</td><td>5</td><td>JRA55 (~60 km, 3 hr)</td><td>(Hutter & Losch, 2020)</td></tr><tr><td>MITTgcm (2 km)</td><td>1997, 2008</td><td>"</td><td>"</td><td>"</td><td>"</td><td>2</td><td>"</td><td>"</td></tr><tr><td>MITTgcm (2 km, e = 1, ↓P)</td><td>1997, 2008</td><td>"</td><td>"</td><td>"</td><td>P*, e = 9.6, 1.0</td><td>"</td><td>"</td><td>"</td></tr><tr><td>MITTgcm (2 km, e = 0.7, ↓P)</td><td>1997, 2008</td><td>"</td><td>"</td><td>"</td><td>P*, e = 9.6, 0.7</td><td>"</td><td>"</td><td>"</td></tr><tr><td>MITTgcm (4.5 km)</td><td>2008</td><td>4.5 km, 240 s</td><td>"</td><td>"</td><td>P*, e = 22.64, 2.0</td><td>"</td><td>ERA-Interim (~80 km, 6 hr)</td><td>(Mohammadi-Aragh et al., 2018)</td></tr><tr><td colspan="9"><strong>McGill-SIM (McGill University)</strong></td></tr><tr><td>McGill (e = 2)</td><td>1997</td><td>10 km, 3600 s</td><td>E</td><td>VP (JFNK)</td><td>P*, e = 27.5, 2.0</td><td>2</td><td>NCEP/NCAR (2.5°, 6 hr)</td><td>(Bouchat & Tremblay, 2017)</td></tr><tr><td>McGill (e = 1, ↓P)</td><td>1997</td><td>"</td><td>"</td><td>"</td><td>P*, e = 13.8, 1.0</td><td>"</td><td>"</td><td>"</td></tr><tr><td>McGill (e = 0.7, ↓P)</td><td>1997</td><td>"</td><td>"</td><td>"</td><td>P*, e = 9.6, 0.7</td><td>"</td><td>"</td><td>"</td></tr><tr><td>McGill (e = 1, ↑S)</td><td>1997</td><td>"</td><td>"</td><td>"</td><td>P*, e = 27.5, 1.0</td><td>"</td><td>"</td><td>"</td></tr><tr><td colspan="9"><strong>NEMO-LIM3/CREG4 (Institut Français de Recherche pour l'Exploitation de la Mer - IFREMER)</strong></td></tr><tr><td>IFREMER (e = 2)</td><td>1997, 2008</td><td>12.4 km, 720 s</td><td>E</td><td>EVP (120)</td><td>P*, e = 20.0, 2.0</td><td>5</td><td>DFS 5.2 (~0.7°, 3 hr)</td><td>(Mulwijk et al., 2019)</td></tr><tr><td>IFREMER (e = 1)</td><td>1997, 2008</td><td>"</td><td>"</td><td>"</td><td>P*, e = 13.8, 1.0</td><td>"</td><td>"</td><td>"</td></tr><tr><td colspan="9"><strong>HYCOM-CICE4 (Florida State University - FSU)</strong></td></tr><tr><td>HYCOM-CICE (FSU)</td><td>1997, 2008</td><td>3.6 km, 360 s</td><td>E</td><td>EVP (120)</td><td>C<sub>pr</sub> e = 19, 2.0</td><td>5</td><td>CFSR/CFSv2 (~38 km, 1 hr)</td><td>(Dukhovskoy et al., 2019)</td></tr><tr><td colspan="9"><strong>HYCOM-CICE4 (Danish Meteorological Institute - DMI)</strong></td></tr><tr><td>DMI</td><td>2008</td><td>9.7 km, 180 s</td><td>E</td><td>EVP (120)</td><td>P*, e = 27.5, 2.0</td><td>5</td><td>ERA-Interim (80 km, 3 hr)</td><td>(Madsen et al., 2016)</td></tr><tr><td colspan="9"><strong>MERCATOR/CREG12 (Mercator Ocean)</strong></td></tr><tr><td>MERCATOR</td><td>1997, 2008</td><td>4.1 km, 900 s</td><td>E</td><td>EVP (500)</td><td>P*, e = 27.5, 1.5</td><td>5</td><td>ERA-Interim (80 km, 3 hr)</td><td>(Dupont et al., 2015)</td></tr><tr><td colspan="9"><strong>NEMO-LIM2/ANHA12 (University of Alberta)</strong></td></tr><tr><td>ANHA (4 km)</td><td>2008</td><td>4.1 km, 180 s</td><td>E</td><td>EVP (120)</td><td>P*, e = 23.4, 2.0</td><td>2</td><td>CGRF (~35 km, 1 hr)</td><td>(Hu et al., 2018)</td></tr><tr><td colspan="9"><strong>NEMO-LIM2/ANHA4 (University of Alberta)</strong></td></tr><tr><td>ANHA (12 km)</td><td>1997, 2008</td><td>12.4 km, 1080 s</td><td>E</td><td>EVP (150)</td><td>P*, e = 23.4, 2.0</td><td>2</td><td>CORE (~200 km, 6 hr)</td><td>(Courtois et al., 2017)</td></tr><tr><td colspan="9"><strong>RIOPS/CREG12-H08 (Environment and Climate Change Canada - ECCC)</strong></td></tr><tr><td>RIOPS</td><td>2008</td><td>4.1 km, 180 s</td><td>E</td><td>EVP (900)</td><td>P*, e = 27.5, 1.5</td><td>10</td><td>CGRF (~35 km, 3 hr)</td><td>(Dupont et al., 2015)</td></tr><tr><td colspan="9"><strong>NEMO-LIM2/ANHA12 (University of Alberta)</strong></td></tr><tr><td>FESOM (Alfred-Wegener-Institute - AWI)</td><td>1997, 2008</td><td>5.1 km, 600 s</td><td>U</td><td>EVP (800)</td><td>P*, e = 27.5, 2.0</td><td>2</td><td>NCEP/NCAR (~1.9°, 24 hr)</td><td>(Wang et al., 2016)</td></tr><tr><td colspan="9"><strong>RASM - Fully Coupled (Naval Postgraduate School - NPS)</strong></td></tr><tr><td>RASM-WRF (EVP)</td><td>1997, 2008</td><td>9.1 km, 1200 s</td><td>E</td><td>EVP (600)</td><td>C<sub>pr</sub> e = 21.3, 2.0</td><td>5</td><td>WRF Model (50 km, 20 min)</td><td>"</td></tr><tr><td>RASM-WRF (EAP)</td><td>1997, 2008</td><td>9.1 km, 1200 s</td><td>E</td><td>EAP (600)</td><td>"</td><td>"</td><td>"</td><td>"</td></tr><tr><td colspan="9"><strong>RASM-CORE2 (Naval Postgraduate School - NPS)</strong></td></tr><tr><td>RASM-CORE2 (EAP)</td><td>1997, 2008</td><td>9.1 km, 1200 s</td><td>E</td><td>EAP (120)</td><td>C<sub>pr</sub> e = 21.3, 2.0</td><td>5</td><td>CORE2 (~110 km, 24 hr)</td><td>"</td></tr><tr><td colspan="9"><strong>neXISIM - VI(2018) (Nansen Environmental and Remote Sensing Center - NERSC)</strong></td></tr><tr><td>neXISIM</td><td>1997, 2008</td><td>10 km, 200 s</td><td>L</td><td>MEB</td><td>T*, P* = 21, 75</td><td>3</td><td>CFSR/CFSv2 (0.5°, 6 hr)</td><td>(Rampal et al., 2019)</td></tr></tbody></table>

Note. We use the following abbreviations (see also Table A1 for the definition of acronyms relating to models, configurations, and atmospheric reanalysis); grid types are E: Eulerian, L: Lagrangian, U: unstructured; ice strength parameters are P*: compressive strength parameter (kPa), T*: isotropic tensile strength parameter (kPa), e: ellipse aspect ratio, C*: frictional energy dissipation parameter, ITD # is the number of ice-thickness categories in the ice thickness distribution; and (solver/nt, it) is the numerical solver used to solve the VP momentum equations (LSR: Line Successive Relaxation solver; JFNK: Jacobian-Free Newton-Krylov solver) or the number of elastic iterations performed to solve the EVP/EAP equations. The grid spacing is given by the mean horizontal grid spacing within the Arctic Ocean. For unstructured grids it refers to the mean node spacing. Run labels including the symbols ↓P or ↑S indicate that the VP elliptical yield curve was modified by reducing the ice strength in isotropic compression, or by increasing the ice shear strength.
the trajectories (see Section 3.1) is necessary to eliminate temporal inconsistencies that can affect the resulting sea-ice deformation statistics (e.g., Bouchat & Tremblay, 2020).

# 3. Methods

## 3.1. Pre-Processing of RGPS Lagrangian Trajectories

To ensure temporal consistency of the RGPS deformation field, we use the Weighted-Average pre-processing method (Bouchat & Tremblay, 2020; Hutter & Losch, 2020), which consists in keeping only trajectories forming cells that have (a) simultaneous (±3 hr) start and end times for all four corners, (b) an average time resolution for all corners that corresponds to the nominal temporal resolution of $T^* = 3$ days, and (c) an area corresponding to the nominal spatial resolution of $L^* = 10$ km. We also require that all corner positions remain at least 100 km away from land for the present analysis. The remaining trajectories are then used to compute the Lagrangian strain rates (see e.g., Equations 1–5 below) and the resulting time series of strain rates for each cell are then averaged in regular 3-day intervals (using a weighted-average of contributing strain rate estimates) starting on 1 January each year. For more information on the pre-processing of the RGPS Lagrangian trajectories and the resulting observed strain rate data set, we refer the reader to Bouchat and Tremblay (2020).

## 3.2. Constructing Simulated Lagrangian Trajectories and Deformation Fields

To construct simulated Lagrangian trajectories and deformation fields from Eulerian model output, we track artificial Lagrangian quadrangle cells that are initialized with the 10 km RGPS Lagrangian positions on 1 January 1997 and 2008. Model trajectories are integrated in their respective grid projection using 1 hr time increments to prevent trajectories from crossing multiple grid cells during one integration step. At the beginning of each hour, the daily mean or snapshot sea-ice velocity field ($u, v$) is first linearly interpolated in time to the current integration time, and then spatially interpolated onto the trajectories' positions using a Great-Circle distance-weighted linear interpolation of the four nearest velocity components (e.g., Madec et al., 2019). Trajectories are integrated until 31 March, unless they drift to within 100 km of the model landmask in which case they are terminated. When the Lagrangian integration is done, the hourly model trajectories are sampled at the beginning of the same regular 3-day intervals as for the RGPS Weighted-Average data set. In the case of data gaps in the RGPS data set, we remove the corresponding simulated trajectory to minimize the effects of the different spatio-temporal coverage on the deformation statistics. Note that the neXtSIM trajectories used in this article were also initialized with the same RGPS cell positions on 1 January 1997 and 2008 and are sampled and masked to the corresponding RGPS coverage following the same constraints as described above for other model trajectories. The neXtSIM trajectories are however integrated within the model's native Lagrangian scheme. In both cases, the initialization of the model trajectories with RGPS positions and the 3-day temporal sampling ensure that the nominal spatial and temporal resolutions of the simulated Lagrangian deformation fields are the same as for the RGPS observations (i.e., $L^* = 10$ km and $T^* = 3$ days), regardless of the original resolution of the model.

The strain rates (velocity gradients) and cell area $A$ are then computed for each cell using the line integral approximations (e.g., Lindsay & Stern, 2003):

$$ \frac{\partial u}{\partial x} = \frac{1}{A} \sum_{k=1}^{4} \frac{1}{2} (u_{k+1} + u_k) (y_{k+1} - y_k), \quad (1) $$

$$ \frac{\partial u}{\partial y} = \frac{-1}{A} \sum_{k=1}^{4} \frac{1}{2} (u_{k+1} + u_k) (x_{k+1} - x_k), \quad (2) $$

$$ \frac{\partial v}{\partial x} = \frac{1}{A} \sum_{k=1}^{4} \frac{1}{2} (v_{k+1} + v_k) (y_{k+1} - y_k), \quad (3) $$

$$ \frac{\partial v}{\partial y} = \frac{-1}{A} \sum_{k=1}^{4} \frac{1}{2} (v_{k+1} + v_k) (x_{k+1} - x_k), \quad (4) $$
with,

$$A = \frac{1}{2} \sum_{k=1}^{4} (x_k y_{k+1} - x_{k+1} y_k), \qquad (5)$$

where $(x_k, y_k)$ is the position of the cell vertex $k$ at time $t$ ($k = 1, 2, 3, 4$; counterclockwise with $x_5 = x_1$ and similar cyclical identities for $y_5$, $u_5$, and $v_5$) and $(u_k, v_k) = \left( \frac{\Delta x_k}{\Delta t}, \frac{\Delta y_k}{\Delta t} \right)$, their approximate velocity during the time interval $\Delta t$. The spatial scale of the strain rate estimate is $L = \sqrt{A}$, and its temporal scale is $T = \Delta t = T^* = 3$ days. Following Bouchat and Tremblay (2020), all cells where $A \le 50$, or $A \ge 200$ km² are removed in order to only keep cells that are representative of the nominal spatial scale ($L^* = 10$ km).

The strain rate invariants (i.e., divergence $\dot{\epsilon}_I$, and shear $\dot{\epsilon}_{II}$) and total deformation rates ($\dot{\epsilon}_{tot}$) are obtained as:

$$\dot{\epsilon}_I = \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y}, \qquad (6)$$

$$\dot{\epsilon}_{II} = \left[ \left( \frac{\partial u}{\partial x} - \frac{\partial v}{\partial y} \right)^2 + \left( \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \right)^2 \right]^{1/2}, \qquad (7)$$

$$\dot{\epsilon}_{tot} = \sqrt{\dot{\epsilon}_{I}^{2} + \dot{\epsilon}_{II}^{2}}. \qquad (8)$$

### 3.3. Deformation Statistics and Associated Metrics

We detail below the deformation statistics used in this study (i.e., PDFs of shear and divergence, spatio-temporal scaling of the mean total deformation rates, and multi-fractal scaling analysis) along with their usual comparison metrics.

#### 3.3.1. Probability Density Functions (PDFs)

The observed RGPS PDFs of strain rate invariants are heavy-tailed and decay approximately linearly in log-log plots, suggesting that large deformations significantly affect the mean deformation statistics (e.g., Bouchat & Tremblay, 2017; Bouillon & Rampal, 2015b; Girard et al., 2011, 2009; Rampal et al., 2019). We investigate the ability of sea-ice models to reproduce the observed large deformation rates by comparing the RGPS and simulated PDFs of shear and divergence. We separate the PDFs of positive divergence (opening) and negative divergence (closing or convergence), as opposed to considering only the absolute divergence PDF as done in previous model evaluation studies (Girard et al., 2011; Lemieux et al., 2012; Rampal et al., 2019; Sreen et al., 2017). We show in Section 4.1.1 that compensation of errors in the positive and negative distributions can lead to a misleading agreement of the simulated absolute divergence PDFs with observations, such that separating the positive and negative divergence PDFs is a preferable way to evaluate the model's ability to reproduce large deformation events in both opening and closing as observed from RGPS.

The PDFs are obtained using logarithmic bins and the typical metric used to compare the observed and simulated PDFs is the decay exponent of the tail, obtained as the slope of a least squares linear fit in log-log space. We do not fix the fitting interval, but rather use an interval of one order of magnitude ending on the largest deformation bin available (or use the maximum available fitting interval if the PDFs spans less than one order of magnitude). We do this because models do not necessarily reproduce deformation rates as large as in the RGPS distributions. Note that, even if they decay approximately linearly in log-log space, it has recently been shown that the tail of the observed RGPS PDFs cannot be considered to follow power law distributions with sufficient evidence, based on a power law maximum likelihood estimator (MLE) fit and a goodness-of-fit test (Bouchat & Tremblay, 2020). Many heavy-tailed distributions that are not power laws can in fact show a tail that approximately follows a straight line in log-log space (e.g., log-normal, exponential, stretched exponential, etc.—Clauset et al., 2009). Trying to find the appropriate distributions representing the PDFs tails (using more rigorous or robust statistical methods such as MLE fits with goodness-of-fit tests) is out of the scope of this study. For this reason, we consider
the decay exponents found here with linear fits in log-log plot as general decay exponents that are only rough approximations of the tails decay, but that are sufficient to draw general conclusions on the usefulness of using a decay exponent as a metric to compare the observed and simulated PDFs.

### 3.3.2. Spatio-Temporal Scaling Analysis of the Mean Total Deformation Rates

We use the coarse-graining procedure with data-quality weights described in Bouchat and Tremblay (2020) to generate deformation fields at larger scales and investigate the spatial and temporal scaling of the mean total deformation rates, that is:

$$ \langle \dot{\epsilon}_{tot}(L, T) \rangle \sim L^{-\beta}, \qquad (9) $$

$$ \langle \dot{\epsilon}_{tot}(L, T) \rangle \sim T^{-\alpha} \qquad (10) $$

where $\langle \cdot \rangle$ denotes the distribution weighted average, L and T are the spatial and temporal scales of the coarse-grained deformation estimates, and $\beta$ and $\alpha$ are the spatial and temporal scaling exponents.

The values of $\beta$ and $\alpha$ characterize the degree of spatial and temporal localization of the mean total deformation rates, and are used as metrics to compare the observed and simulated spatio-temporal scaling. The spatial scaling exponent $\beta$ varies between 0 (deformation field is homogeneous in space) and 2 (deformations are highly localized in space), while the temporal scaling exponent $\alpha$ varies between 0 (deformations field is homogeneous in time) and 1 (deformations are highly localized in time). Here, $\beta$ and $\alpha$ are obtained using least squares power law fits on the average total deformation rates $\langle \dot{\epsilon}_{tot}(L, T) \rangle$ for a given range of spatial and temporal scales. We restrict the spatial scaling to scales $10 \le L \le 600$ km and the temporal scaling to scales $3 \le T \le 30$ days to minimize the effects of the reduced spatio-temporal coverage at larger scales.

Using data-quality weights to obtain the distribution average at each scale results in giving more weight to the tail of the distribution where interesting sea-ice dynamical features are represented (e.g., LKFs) and where deformation rates have smaller relative errors. As discussed in Section 4.1.2, this improves the interpretation of the scaling exponent metrics as a measure of localization of deformations when applied to simulated deformation fields. The data quality is defined by the signal-to-noise ratios which are obtained by estimating the Lagrangian trajectory errors. The details of the signal-to-noise ratio calculations for RGPS and simulated Lagrangian trajectories can be found in Appendix B.

### 3.3.3. Multi-Fractal Analysis

The spatio-temporal scaling analysis described for the mean total deformation rate above is repeated for higher moments $q$ to construct $\beta(q)$ and $\alpha(q)$, the spatial and temporal structure functions, that is:

$$ \langle \dot{\epsilon}_{tot}^q(L, T) \rangle \sim L^{-\beta(q)}, \qquad (11) $$

$$ \langle \dot{\epsilon}_{tot}^q(L, T) \rangle \sim T^{-\alpha(q)}. \qquad (12) $$

It has usually been assumed that the structure functions $\beta(q)$ and $\alpha(q)$ for sea-ice total deformation rates are quadratic, for example, $\beta(q) = aq^2 + bq$, where $a$ has been interpreted as the degree of multi-fractality of the scaling (e.g., Bouchat & Tremblay, 2020; Bouillon & Rampal, 2015b; Hutter et al., 2018; Rampal et al., 2019). However, following the universal multi-fractal formalism, the structure functions are not required to be quadratic and can have a varying degree of nonlinearity, which is then more correctly interpreted as the degree of multi-fractality (Lovejoy & Schertzer, 2007, 2013). Here, we do not assume a fixed degree of multi-fractality and instead find a general least squares fit for the structure functions of the following form (in full agreement with the universal multi-fractal formalism—e.g., Lovejoy & Schertzer, 1995, 2007; Weiss, 2008):

$$ \beta(q) = q - \zeta(q), \qquad (13) $$

with,

$$ \zeta(q) = qH - K(q), \qquad (14) $$

and
$$K(q) = \frac{C_1}{\mu - 1} (q^{\mu} - q), \qquad (15)$$

such that we can write:

$$\beta(q) = \left( \frac{C_1}{\mu - 1} \right) q^{\mu} + \left( 1 - \left( H + \frac{C_1}{\mu - 1} \right) \right) q, \qquad (16)$$

where H is a fluctuation exponent, K(q) is the universal multi-fractal moment scaling function, C₁ (0 ≤ C₁ ≤ 2) characterizes the degree of heterogeneity (or sparseness) of the field, and μ (0 ≤ μ ≤ 2) is the degree of multi-fractality (μ = 0 for a mono-fractal process, and μ = 2 for a log-normal multiplicative model with maximal degree of multi-fractality). An equivalent formulation applies for the temporal structure function α(q).

In the following, the values of the three multi-fractal parameters H, C₁, μ are used as metrics to compare the observed and simulated multi-fractal structure functions.

# 4. Results

Results for low-resolution model runs (Δx = 9–12 km) are presented separately from high-resolution model runs (Δx = 2–5 km), even if their Lagrangian deformation fields are reconstructed at the same nominal spatial scale of L* = 10 km. In fact, Eulerian models with finite-difference schemes will resolve the sea-ice dynamics with different levels of complexity as their spatial resolution changes (e.g., Sreen et al., 2017; Williams & Tremblay, 2018). It is therefore expected that higher resolution runs will resolve finer deformation features in their Lagrangian deformation fields, affecting the result of the deformation metrics. For instance, consider the observed sea-ice deformation field sampled at L* = 10 km. The deformation statistics at this scale are the result of underlying dynamics occurring at much finer scales (e.g., fractures at the sub-km scales). The observed deformation fields sampled at L* = 10 km are therefore much more rich in information than model deformation fields that are generated (rather than sampled) at the same nominal spatial scale, unless sub-grid parameterization are used and calibrated. Degrading the observed deformation fields to larger spatial scales could help minimizing this discrepancy when comparing the observed and simulated deformation statistics, but only if the degraded spatial scales are much larger than the nominal spatial scales at which models are run (e.g., observations at L ~50–100 km vs. models at L ~10 km), in which case the range of scales available for determining the observed statistical characteristics (e.g., spatio-temporal localization) becomes too small. Note that we also consider atmosphere-ice-ocean coupled model simulations with an interactive atmospheric model (i.e., RASM-WRF) separately from coupled ice-ocean models (or stand-alone ice models) forced with atmospheric reanalyses (see Section 4.2.3).

In the following sections, the agreement between models and observations is interpreted in terms of the RGPS interannual variability. That is, metrics are first obtained for all years in the RGPS record and, unless stated otherwise, the full RGPS distribution is used as a range defining a good agreement between models and observations.

## 4.1. Effects of Sea-Ice Rheology

### 4.1.1. Probability Density Functions

Most of the simulated PDFs of shear and divergence decay approximately linearly in log-log plot, with a wide range of simulated decay exponents (Figures 1a–1c, 2a–2c, 3a–3c and 4a–4c). We note that very different distributions can lead to very similar decay exponents, suggesting that this metric does not adequately capture differences in the deformation fields (e.g., compare RGPS with HYCOM-CICE (FSU) in Figure 3a, or with IFREMER (e = 1) in Figure 2a). We, therefore, define a new metric as the sum of the absolute difference between the simulated and observed PDFs in logarithmic scale, divided by the number of bins spanned by the simulated PDF. Dividing by the number of bins ensures that the metric penalizes models that do not simulate sufficiently large deformation rates and have a smaller number of bins. An advantage of this metric is that differences in the tail of the PDFs (i.e., where probabilities are small, but represent larger deformation rates that are likely to affect climate interactions or operational applications) are given more weight by using a logarithmic scale. To interpret the value of the metric, we compute its interannual variability for all available RGPS observations, using either
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Probability density functions (PDFs) of (a): maximum shear rate, (b): negative divergence, and (c) positive divergence, for low-resolution runs ($\Delta x \simeq 10$ km) and RADARSAT Geophysical Processor System (RGPS) observations (black) at $L = 10$ km and $T = 3$ days in January-February-March 1997. Insets give the decay exponents ($k_{\epsilon_{11}}$, $k_{\epsilon_{10}}$, $k_{\epsilon_{10}}$), and each run is identified in panel (d) by a different color corresponding to its sea-ice rheology (Blue/Purple: VP, Pink/Red: EVP, Yellow: EAP, Green: MEB); (d-f): difference (per bin) between the logarithm of models and RGPS PDFs. The insets give the average absolute difference per bin, where bold font marks values that are equal to or less than the RGPS interannual average obtained using all other RGPS years in comparison to 1997.</figcaption></figure>

the RGPS PDFs of 1997 or 2008 as the reference and computing the difference metric with all other years in the RGPS data set. We then use the mean value of the RGPS PDFs difference metric for each comparison year (one value for 1997 and another for 2008) as an upper threshold defining a reasonable agreement between models and observations. These references mean values, as well as the mean difference per bin (in logarithmic scale) for the RGPS data set are shown in Figures 1d-1f, 2d-2f, 3d-3f, and 4d-4f for comparison.

Considering this new metric, we find that no low-resolution run is able to reproduce a reasonable agreement with the RGPS negative divergence PDFs, whereas only neXtSIM and McGill ($e = 1, \uparrow S$) agree with either the shear or positive divergence PDFs (Figures 1d-1f and 2d-2f). This reflects a clear underestimation of the range over which the PDFs extend (i.e., smaller number of bins), along with a drop in probabilities in the respective last bins. Modifying the (E)VP plastic elliptical yield curve parameters at low resolution helps increasing the range over which the PDFs extend and also reduces the drop in the tail (Figure 1—see also Section 4.2.1).

Increasing the spatial resolution of the models generally improves the agreement of simulated PDFs with observations, particularly in negative divergence (Figures 3 and 4). This can be attributed in part to a refinement of LKFs in which deformation rates increase at higher resolution (e.g., Spreen et al., 2017; Williams & Tremblay, 2018), and in part to an increased spatio-temporal LKF density (Hutter et al., 2022a). We find that high-resolution runs with the (E)VP rheology (i.e., the only rheology represented by the high-resolution runs) can reproduce PDFs that simultaneously agree with the RGPS shear and positive and/or negative divergence. In fact, the MITgcm (2 km) run agrees with all three RGPS PDFs (i.e., in shear, positive divergence, and negative divergence) for at least one year (2008—Figure 4). However, some EVP runs at high-resolution still poorly agree with the RGPS PDFs, even if the range of the simulated PDFs is improved compared to low-resolution EVP runs (see e.g., HYCOM-CICE (FSU) and ANHA 4 km). We hypothesize that this reflects a numerical artifact originating from insufficient subcycling with the EVP approach. In the EVP equations, an artificial elastic strain is added to the VP rheology to allow explicit solving of the momentum equations. Within each advective time step, small iteration steps (subcycling) are used to explicitly advance the solution, while damping the artificial elastic waves in order to recover a solution that approximates a VP solution. When using too few subcycles with the EVP solver, the solution is noisy with large residual errors, and the probability of simulating large deformation rates is significantly reduced (Kimmritz et al., 2015; Lemieux et al., 2012). While we cannot explicitly check their numerical convergence, we
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Same as Figure 1 for low-resolution runs ($\Delta x \approx 10$ km) in January-February-March 2008.</figcaption></figure>

note that noise is present in the EVP deformation fields analyzed here (results not shown), and that the EVP runs poorly agreeing with the observed RGPS PDF consistently use a small number of subcycles (and vice-versa: EVP runs showing a good PDF agreement also use a large number of subcycles—see e.g., ANHA 4 km: 120 subcycles, vs. RIOPS: 900 subcycles in Figure 4). We, therefore, hypothesize that the high-resolution EVP runs showing a poor PDF performance here are also affected by large residual errors originating from undamped elastic waves and too few subcycles. This could explain the lower performance of EVP compared to VP for low-resolution runs as well, but it remains to be validated with further experiments.

We lastly note that the conclusions of the PDF analysis could be significantly different if the PDFs of absolute divergence were compared, as opposed to separating the negative and positive divergence PDFs as done here. For instance, simulations that closely agree with the RGPS PDF of absolute divergence may hide errors in their positive and negative divergence PDFs that compensate each other, suggesting that the conditions or processes giving rise to large opening vs. closing events are not well represented in the model and tend to overestimate one in favor of the other (see e.g., neXtSIM and MITgcm (2 km, $e = 0.7$, ↓P) in Figure 5). Similarly, simulations that

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Same as Figure 1 for high-resolution runs ($\Delta x \approx 2–5$ km) in January-February-March 1997.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Same as Figure 1 for high-resolution runs ($\Delta x \approx 2–5$ km) in January-February-March 2008.</figcaption></figure>

closely match the observed RGPS PDF in either positive or negative divergence but not in the other (see e.g., McGill, *e* = 1, ↑S in Figure 5) will be strongly penalized when comparing the absolute divergence PDF and good model performance may be missed. We argue that model evaluation or comparison should therefore be done by investigating both the negative and positive divergence PDFs, in order to better evaluate these subtle differences that are likely affected by the representation of sea-ice rheology (i.e., ice strength, yield curve, and flow rule).

### 4.1.2. Spatio-Temporal Scaling and Coupling

The spatio-temporal scaling analysis of simulated deformation rates has typically been investigated without using data-quality weights (Bouchat & Tremblay, 2017; Bouillon & Rampal, 2015a; Girard et al., 2009; Hutter et al., 2018; Sreen et al., 2017; Rampal et al., 2019). In this case, the spatio-temporal scaling exponents depend on the LKF density, such that runs with fewer LKFs have lower scaling exponents (Hutter & Losch, 2020). However, lower scaling exponents are also expected for diffuse deformation fields (see Section 3.3.2), complicating the interpretation of low scaling exponents within an inter-model comparison analysis. For example, comparing runs with few LKFs but highly localized deformations (e.g., FESOM—Figure 6a) with runs showing obviously more diffuse deformations (e.g., HYCOM-CICE [FSU]—Figure 6b), we find here that the slope of the spatial scaling (i.e., the spatial scaling exponents) are comparable when the data quality (signal-to-noise ratio) is not used to weight the deformation estimates in the scaling analysis (Figure 6c). In contrast, using the signal-to-noise ratios to weight the simulated deformation distribution helps to separate both cases, as the slope increases for simulations with highly localized deformation features, while it remains low for more diffuse deformation fields (Figure 6d). Implementing the scaling analysis with signal-to-noise ratio weights to compare observations and models therefore improves the interpretation of the scaling exponent metric as a measure of the localization of the deformation fields. It also allows us to investigate the presence of a spatio-temporal coupling of the spatial and temporal scaling exponents (i.e., a logarithmic decay of $\beta$ and $\alpha$ when increasing *T* and *L*, respectively—Marsan & Weiss, 2010), which is otherwise absent for the observed RGPS mean total deformation rates when using weights equal to one (Bouchat & Tremblay, 2020). In the following, the scaling analysis is therefore performed with the signal-to-noise ratio weighting method. Note that, while the coupling and scaling exponents are affected, we have verified that finding a power law scaling in space or time does not depend on the weights used to average the deformation distribution (i.e., signal-to-noise ratio weights vs. weights equal to one as in previous studies).

We find that all sea-ice rheologies produce a power law spatial scaling of the total deformation rates holding over ~1.5 orders of magnitude (i.e., $10 \le L \le 600$ km—Figures 7a, 7b and 8a, 8b). However, the simulated spatial scaling exponent $\beta$ (i.e., the slope of the power law decay in log-log space) varies largely from run to run (Figures 7c, 7d and 8c, 8d). We note that the only runs showing a spatial scaling exponent large enough to be within the
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Probability density functions (PDFs) of (a): negative divergence, (b): positive divergence, and (c) absolute divergence, for selected runs and RGPS observations (black) at $L = 10$ km and $T = 3$ days in January-February-March 1997; (d-f): difference (per bin) between the logarithm of models and RGPS PDFs. The insets give the average absolute difference per bin, where bold font marks values that are equal to or less than the RGPS interannual average obtained using all other RGPS years in comparison to 1997.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Top: Total deformation field snapshot for 28–30 January 1997 (in day⁻¹) for (a): FESOM, and (b): HYCOM-CICE (FSU); bottom: spatial scaling of total deformation rates estimated at $T = 3$ days in January-February-March 1997, using either (c): weights equal to one, or (d): weights equal to the signal-to-noise ratios (data quality) of the deformation estimates.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 7.** (a, b): Spatial scaling for total deformation rates estimated at *T* = 3 days in January-February-March 1997. (c, d): Spatial scaling exponent *β* as a function of the temporal scale *T* at which the mean total deformation rates are estimated. (e, f): Coupling coefficient *c*<sub>β</sub> obtained from least squares logarithmic fits *β* ~ -*c*<sub>β</sub> ln(*T*) for 3 ≤ *T* ≤ 30 days. Dashed lines are the least squares power law fits used to obtain *β*. The solid black lines, dark gray, and light gray shaded areas are the mean, standard deviation, and min/max for the entire RGPS data set. Model results are separated with low-resolution runs in top panels, and high-resolution runs in bottom panels.</figcaption></figure>

observed RGPS interannual variability or larger (i.e., neXtSIM, RIOPS, FESOM, and MITgcm—2 km runs) also showed a reasonable agreement in their PDFs of deformations. The presence of large deformation rates therefore appears as a necessary condition for also having a large degree of spatial localization. It is not sufficient however, since it is the spatial organization of these large deformation rates along well-defined features (i.e., LKFs) that is responsible for the spatial scaling (e.g., Marsan et al., 2004; Stern & Lindsay, 2009).

For low-resolution runs, the largest spatial scaling exponents are obtained with the MEB rheology (neXtSIM). While the neXtSIM deformation fields do show highly localized LKFs (Figure 9), this model uses an adaptive Lagrangian mesh as opposed to a static Eulerian grid as in all other runs. It is therefore not straightforward to attribute this stronger spatial localization of deformation to the rheology alone since moving meshes are known to be very efficient at capturing and preserving singularities or discontinuities in the solution (e.g., Ceniceros & Hou, 2001).

The lowest spatial scaling exponents are obtained with the EVP rheology, in both low- and high-resolution runs (Figures 7c, 7d and 8c, 8d). The deformation fields for these runs (i.e., DMI, IFREMER, HYCOM-CICE (FSU), ANHA 4 and 12 km) clearly underestimate the presence of well-defined deformation features (Figures 9 and 10 and Hutter et al., 2022a). We again hypothesize that this could be due to insufficient damping of the artificial elastic wave and small numbers of subcycling steps, but the effects of the numerical convergence on the scaling statistics need to be further evaluated. We note, however, that more iterations to obtain more accurate VP and EVP solutions leads to additional lines of deformations in the solution (Bouchat & Tremblay, 2014; Koldunov et al., 2019; Lemieux et al., 2012; Wang et al., 2016), which should increase the spatial scaling exponent (or spatial localization).

We note that the RASM-WRF (EAP) runs simulate more localized deformations than its RASM-WRF (EVP) counterparts and all other low-resolution (E)VP runs, as indicated by the larger spatial scaling exponents and the presence of well-defined LKFs (Figures 7a, 7c, 8a, 8c and 10). While the larger number of elastic subcycles in the RASM-WRF runs may play a role in the improved scaling statistics as mentioned above, this may also hint
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Same as Figure 7 for total deformation rates in January-February-March 2008.</figcaption></figure>

at an advantage of the anisotropic sub-grid parameterization for a better localization of deformations. However, grid aliasing (i.e., orthogonal alignment of LKF along the RASM atmospheric or ocean grid) is also present in the RASM deformation fields (Figure 10) and the impacts of such a preferred alignment of LKF on the scaling statistics have not yet been investigated. The improved spatial localization for RASM-WRF (EAP) nonetheless suggests that a more detailed analysis of the potential advantages of using the EAP rheology compared to the classic (E)VP rheology would be a welcome contribution in the future.

Ignoring the runs with a smaller number of subcycles mentioned above, the spatial scaling exponent for (E)VP runs generally increases as the grid is refined. This is consistent with a refinement of the spatial localization of deformation lines with increasing spatial resolution in Eulerian plastic sea-ice models (Williams & Tremblay, 2018). In contrast, the spatial scaling exponent was shown to be approximately resolution-independent for neXtSIM (MEB) when tested on a range of spatial resolutions from 30 to 7.5 km (Rampal et al., 2019). It is still unclear whether this is a consequence of using a Lagrangian mesh that better adapts to discontinuities in the solution (regardless of the resolution), or of using a brittle rheology. We can however conclude that a large spatial localization of deformation is possible for both visco-plastic (E)VP and brittle visco-elastic (MEB) rheologies, as long as Eulerian sea-ice models are run at high spatial resolution. Modifying the ice strength parameters and coupling with an atmospheric model can also have a large effect on increasing the scaling exponents as discussed later in Sections 4.2.1 and 4.2.3.

Interestingly, both low- and high-resolution runs span a similar range of temporal scaling exponents that overlaps with the RGPS interannual variability, showing that a strong degree of temporal localization of deformations is reproduced by all models, at least for the range of temporal scales considered in this study (i.e., [3–30] days— Figures 11 and 12). Here, it is important not to confuse strong temporal localization with strong intermittency. A field can be highly localized in time, but it is the change of localization within the data set (or with changing deformation magnitude) that reflects the intermittency (or heterogeneity). The intermittency of the deformation field is indicated by its (nonlinear) moment scaling function (or *structure function*) which is investigated in Section 4.1.3. Temporal scaling (or localization in time) of the deformation rates is assumed to originate from the presence of long-ranged temporal correlations in the time series of deformations. We have verified that when randomly re-ordering the times series of deformation, the power law temporal scaling is lost for both
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Total deformation rate snapshots (in day⁻¹) for selected runs for the period of 21-22-23 February 2008.</figcaption></figure>

RGPS observations and simulated deformation fields (results not shown). This is analogous to the presence of long-ranged spatial correlation (for instance, LKFs) giving rise to the spatial scaling. The origin of these temporal correlations in models and observations remains to be identified. We note however that a larger simulated temporal scaling exponent does not necessarily correlate with the use of a smaller advective time-step, nor with higher spatio-temporal resolution of the atmospheric reanalysis. Preliminary analysis with the MEB rheology (not shown) also shows that the choice of damage propagation scheme can significantly affect the spatio-temporal scaling and could be used to tune this rheology against observations.

Finally, a logarithmic reduction in the spatial and temporal scaling exponents when increasing the temporal and spatial scales of the deformation estimates (i.e., $\beta \sim -c_{\beta} \ln(T)$ and $\alpha \sim -c_{\alpha} \ln(L)$, the so-called *space-time coupling*) is achieved by all sea-ice rheologies, regardless of the original spatio-temporal resolution of the model runs (Figures 7, 8, 11, 12c and 12d). This indicates that the simulated deformation fields appear less and less localized as the spatial and temporal scales are increased, consistent with the smoothing of deformation features when averaged at larger and larger scales. The strength of the observed RGPS coupling, evaluated by the coupling constants $c_{\beta}$ and $c_{\alpha}$ (i.e., absolute value of the slope in semi-log plot), is also well reproduced by all rheologies at both low- and high-resolution Figures 7, 8, 11, 12e and 12f). Runs for which the space-time coupling is systematically absent or very weak (e.g., IFREMER $e = 2$ and ANHA 12 km) are low-resolution EVP runs and already
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Total deformation rate snapshots (in day⁻¹) for selected runs for the period of 10-11-12 January 1997.</figcaption></figure>

have smoother deformation fields to start with. Marsan and Weiss (2010) suggested that a space-time coupling of sea-ice deformation scaling can emerge from brittle dynamics and a possible chain-triggering deformation mechanism similar to that observed for earthquakes. We show here that sea-ice rheologies that do not assume brittle parameterizations also reproduce such a coupling.

4.1.3. Multi-Fractal Analysis

As the moment *q* of the total deformation distribution increases, the scaling exponents *β*(*q*) and *α*(*q*) also increase, given that the scaling still holds. For mono-fractal systems, the increase in localization is linear with increasing moment, while for multi-fractal systems, the increase in localization with increasing moment deviates from linearity. Multi-fractality then reflects a large variability of the scaling exponent within the field. For sea-ice deformation fields, multi-fractality can be interpreted as larger deformation rates being more localized (in space and time) than smaller deformation rates (Rampal et al., 2019; Weiss & Dansereau, 2017).

Using the universal multi-fractal formalism, the nonlinear multi-fractal structure functions are described by three variables: the degree of multi-fractality *μ*, the degree of heterogeneity $C_1$, and the fluctuation exponent *H* (see Equation 16). The spatial scaling exponent of the mean total deformation rates evaluated in the previous section
<figure><img src="image_14.png" /></figure>

**Figure 11.** (a, b): Temporal scaling for total deformation rates estimated at *L* = 10 km in January-February-March 1997. (c, d): Temporal scaling exponent *α* as a function of the spatial scale *L* at which the mean total deformation rates are estimated. (e, f): Coupling coefficient *c*<sub>α</sub> obtained from least squares logarithmic fits *α* ~ *c*<sub>α</sub> ln(*L*) for 10 ≤ *L* ≤ 300 km. Dashed lines are the least squares power law fits used to obtain *α*. The solid black lines, dark gray, and light gray shaded areas are the mean, standard deviation, and min/max for the entire RGPS data set. Model results are separated with low-resolution runs in top panels, and high-resolution runs in bottom panels.

<figure><img src="image_15.png" /></figure>

**Figure 12.** Same as Figure 11 for total deformation rates in January-February-March 2008.
is equal to $\beta(1) = 1 - H$ and therefore, the larger the $H$, the smoother (or less localized) the field appears. Interpretation of the effects of $\mu$ and $C_1$ on the observable fields are less intuitive. Generally, a larger value of $\mu$ characterizes a field dominated by singularities of larger values, and a larger $C_1$ indicates that these singularities are more sparsely grouped (Lovejoy & Schertzer, 2007, 2013). However, for the same values of $\mu$, $C_1$, and $H$ the field-to-field variability can be large (Lovejoy & Schertzer, 2013) and it is not straightforward to visually distinguish the effects of the different parameters. We can nonetheless identify a few general points below regarding the use of the structure functions and the multi-fractal parameters as deformation metrics for evaluating sea-ice models. Note that the universal multi-fractal formalism requires $0 \le C_1 \le 2$ and $0 \le \mu \le 2$, but there is no theoretical constraint on $H$. However, in practice, the way that the fluctuations of the analyzed field are obtained at different scales (e.g., using first-order differences, or wavelets analysis, etc.) will limit the range of $H$ over which multi-fractality can be detected (see e.g., Lovejoy & Schertzer, 2013). For velocity derivatives obtained with first-order differences such as in our case, this range is $0 < H < 1$. Given that all scaling exponents found for models and observations correspond to an equivalent $H$ within this range, we can therefore expect to detect multi-fractality, if it applies. As such, we do not constrain the least squares fits used to obtain the multi-fractal parameters to return values in those theoretical intervals, allowing us to evaluate the validity of the multi-fractal hypothesis for the observed and simulated deformation fields.

All sea-ice rheologies reproduce nonlinear structure functions in space and time, suggesting that multi-fractality (i.e., $\mu \neq 0$) and heterogeneity ($C_1 \neq 0$) are not exclusive to a specific rheology assumption (Figures 13 and 14 and Figures 15 and 16). In general, the conclusions of the previous section based on the scaling of the mean ($q=1$) total deformation rates also apply to $q > 1$. These conclusions include: higher scaling exponent for MEB and high-resolution models, lower scaling exponents for EVP runs with fewer subcycles, larger variability of spatial scaling exponents compared to temporal scaling exponents. We note however that agreement with the RGPS interannual variability for $q=1$ does not necessarily carry over to higher moments (see e.g., $\alpha(q)$ for RIOPS in Figure 14). In fact, models agreeing with the RGPS distribution for the fluctuation exponent $H$ (i.e., for the scaling of the mean, or $q=1$) do not necessarily agree in the other multi-fractal parameters describing the structure functions, and vice-versa (Figures 15 and 16). Models agreeing with the RGPS distribution of $\beta(q)$ or $\alpha(q)$ for all $q$'s and for at least one year (i.e., McGill ($e=0.7, \downarrow P$), IFREMER ($e=1$), neXtSIM, MITgcm (2 km), MITgcm (2 km, $e=1, \downarrow P$), and RIOPS), also generally agree in their multi-fractal parameters (i.e., $\mu, C_1$, and $H$). Other models with all three multi-fractal parameters within the RGPS inter-annual variability for at least one year include McGill ($e=1, \uparrow S$), RASM-CORE2 (EAP), and MERCATOR. However, we note that the spatial and temporal multi-fractality hypothesis for RGPS observations is not robust since the distribution of the fitted degree of multi-fractality ($\mu$) reaches values outside the theoretical range, which complicates the comparison and interpretation of the observed and simulated multi-fractal parameters (e.g., in 2008—Figure 16). In this case, the usefulness of the multi-fractal structure functions to evaluate sea-ice deformation fields is not clear and more work is required to better understand why the multi-fractal hypothesis is not valid for certain years. Nevertheless, we note that the degree of multi-fractality ($\mu$) for other years of the RGPS records is generally not quadratic (i.e., $\mu \neq 2$). This confirms that all three multi-fractal parameters should be used as metrics for the structure functions, as opposed to considering a fixed (e.g., quadratic) degree of multi-fractality and using only the degree of heterogeneity as a metric.

## 4.2. Effects of Model Configuration and Other Parameterizations

Results from the previous section show that deformation statistics have a run-to-run variability that can be as large or larger than the effects of the choice of a given sea-ice rheology. In the present section, we explore model parameterizations that could explain part of this variability.

### 4.2.1. Ice Strength Parameters

In the classical two ice-categories (E)VP rheology, the ice strength is parameterized using an elliptical yield curve and a compressive ice strength parameter $P^*$, which defines the maximum isotropic compressive stress that can be supported by ice for a given thickness and concentration (Hibler, 1979). The elliptical yield curve then implicitly defines the shear strength parameter $S^*$ of the ice through the ratio of the major to minor axes, that is, the ellipse ratio $e$ (Bouchat & Tremblay, 2017). Calibration of the ellipse ratio and compressive ice strength parameter have usually been performed by minimizing the drift and/or thickness errors (e.g., Hibler & Walsh, 1982; Miller et al., 2006; Ungermann et al., 2017). However, the PDFs of sea-ice deformation rates are sensitive to
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Left: RGPS (black) and simulated (colors) spatial structure functions $\beta(q)$ for total deformation rates estimated at $T = 3$ days in January-February-March 1997. Right: temporal structure functions $\alpha(q)$ for total deformation rates estimated at $L = 10$ km in January-February-March 1997. Dashed lines are the least squares fit to Equation 16 used to derive the degree of multi-fractality $\mu$, the degree of heterogeneity $C_1$, and the fluctuation exponent $H$. The solid black lines, dark gray, and light gray shaded areas are the mean, standard deviation, and min/max for the entire RGPS data set. Model results are separated with low-resolution runs in top panels, and high-resolution runs in bottom panels.</figcaption></figure>

independent changes of $P^*$ or $S^*$, and therefore it has been suggested that observed RGPS PDFs of deformation rates can be used to calibrate the ice strength parameters in sea-ice models (Bouchat & Tremblay, 2017). Specifically, increasing the ratio of shear-to-compressive ice strength parameters (i.e., reducing the ellipse ratio from $e = 2$, $P^* = 27.5 \text{ kNm}^{-2}$ to $e = 1$, $P^* = 13.8 \text{ kNm}^{-1}$) significantly improved the agreement between observed and simulated PDFs of deformation rates for VP gridded deformation fields at low (10 km) resolution. Other studies using the (E)VP rheology with a reduced ellipse ratio (i.e., $0.7 \le e \le 1.8$) at low resolution also showed improved landfast ice and ice bridges simulation, as well as reduced ice thickness bias (Dumont et al., 2009; Lemieux et al., 2016; Miller et al., 2005). Whether these conclusions are configuration-dependent (e.g., resolution, forcing, ridging scheme, etc.) has however not been tested.

We revisit the McGill runs (with the same parameters as in Bouchat & Tremblay, 2017) in order to investigate the sensitivity of the deformation statistics to the ice strength parameters with our updated deformation metrics, which now include temporal scaling, multi-fractal structure functions, and the new PDF-difference metric. We also extend this analysis to the IFREMER runs (low-resolution) and MITgcm 2 km runs (high-resolution), where only the compressive ice strength parameter $P^*$ and the ellipse aspect ratio $e$ were modified. At low resolution (McGill and IFREMER runs), the results confirm that increasing the ratio of shear-to-compressive strength parameter can improve the agreement of all simulated deformation statistics with RGPS observations, independently of the model configuration. The PDF-difference metric reveals that reducing the ice strength in compression even lower than suggested in Bouchat and Tremblay (2017) provides a better agreement with the RGPS distributions, at least in shear and positive divergence (see e.g., McGill $e = 0.7, \downarrow P$ in Figure 1). We also note that the spatio-temporal scaling analysis of Lagrangian trajectories with signal-to-noise ratios as weights is more conclusive than the gridded scaling analysis in Bouchat and Tremblay (2017). The results show that an
<figure><img src="image_17.png" /><figcaption>**Figure 14.** Same as Figure 13 for total deformation rates in January-February-March 2008.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 15.** Multi-fractal parameters ($\mu$, $C_1$, $H$—see Equation 16) for the spatial structure function $\beta(q)$, and for the temporal structure function $\alpha(q)$, for runs in 1997 and RGPS inter-annual variability (boxplots). Dashed areas represent parameters outside the valid range predicted by the multi-fractal formalism. Model results are separated with low-resolution runs in top panels, and high-resolution runs in bottom panels.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Figure 16.** Same as Figure 15, for runs 2008.</figcaption></figure>

increase in the shear-to-compressive strength ratio (either by reducing *P*\* or increasing *S*\*) systematically leads to spatial and temporal scaling exponents closer to those for the RGPS observations for both the McGill and IFREMER low-resolution runs (Figures 7, 8, 11, and 12). The analysis of the structure functions also reveals that the degree of heterogeneity and intermittency (*C*₁) in space and time is sensitive to changes in the shear-to-compressive strength ratio with the (E)VP rheology (Figures 15 and 16).

At high-resolution (MITgcm 2 km runs), increasing the shear-to-compressive strength ratio can also improve the sea-ice deformation statistics (Figures 3, 4, 7, 8, 11, 12, 13 and 14). However, the effects on the degree of heterogeneity and multi-fractality are less clear than at low resolution (Figures 15 and 16). We also note that the same changes in *e* and *P*\* or *S*\* that improve the positive divergence PDFs in low-resolution runs instead result in a clear overestimation of positive divergence in high-resolution runs (Figures 3 and 4). In this sense, it is not clear that the combination of ice strength parameters that provides the best model-observation agreement for low-resolution runs are also indicated for high-resolution runs. These results point at the need to better understand how sea-ice parameterizations should (or not) change with changing model resolution given the intricate links that exist between the (E)VP yield curve (and flow rule), the deformation fields, the energy dissipation, and the model resolution.

Finally, we note that the yielding shear, compressive, and tensile strength are much larger for the Mohr-Coulomb yield curve in the neXtSIM model than for typical plastic elliptical yield curves (see Table 1). In VP models at low-resolution, a higher shear strength allows the stress level to increase within the ice and to be relieved along well-defined and less frequent (more intermittent) deformation features, which helps improving the simulated deformation statistics (Bouchat & Tremblay, 2017). Whether this is also the case in the MEB rheology and could also partly explain the better deformation statistics of the neXtSIM model at low resolution compared to other (E) VP models remains to be verified.

### 4.2.2. Ice Thickness Distribution

The simplest way to represent the presence of ice in a continuum sea-ice model is to use two categories of ice thickness: thick ice, and thin or no ice. The ice is then characterized by its mean thickness (*h*) and concentration (A) per grid cell, and the ice strength is typically assumed to depend linearly on *h* (Hibler, 1979). However, as multiple sub-grid scale processes in the Arctic climate system are affected by the local presence of thick vs. thin ice (e.g., albedo, conductive heat fluxes, etc.), it is now common practice to use an ice-thickness distribution
(ITD) with more than two thickness categories (Thorndike et al., 1975). In this case, the ice strength can instead be parameterized as a function of the change in potential energy during the ridging process (Rothrock, 1975), which explicitly depends on the thinnest ice category and on the local distribution of ice in the different thickness categories. This change in the ice strength formulation was shown to increase the spatial heterogeneity of the simulated ice strength and to significantly increase the deformation rates in convergence for thick multiyear ice in a very-low resolution (Δx ~ 36 km) coupled ice-ocean model (Ungermann et al., 2017). On the other hand, ITD runs with the simpler ice strength parameterization of Hibler (1979) perform better than runs with the strength parameterization of Rothrock (1975) at reproducing large-scale observations of sea-ice concentration, drift, and thickness (Ungermann et al., 2017). Hutter and Losch (2020) recently showed that using the ice strength parameterization of Rothrock (1975) with a multiple-category ITD results in a larger number of LKFs in high-resolution runs compared to using a two-category model, however, there was no comparison to an ITD run with the strength parameterization of Hibler (1979). A thorough discussion of the effects of using an ITD and of the different ice strength parameterizations on the sea-ice deformation statistics is therefore still needed.

Here, the 2008 MITgcm 2 km runs (one with two thickness categories and the other with an ITD) allow us to investigate the effects of introducing multiple sub-grid ice categories on the deformation statistics within the same model, in light of the new PDF-difference metric introduced in Section 4.1.1 and the updated scaling analysis with signal-to-noise ratio weights. However, we still cannot firmly comment on the effects of the different ice strength parameterizations because the MITgcm (2 km, ITD) run uses the Rothrock (1975) strength formulation and the two-category MITgcm (2 km) run uses the Hibler (1979) formulation and no other model provided ITD runs with both strength formulations. On the one hand, we find that there is no clear improvement in the agreement of the simulated PDFs of shear rates and divergence with observations when introducing an ITD at high resolution (Figure 4). On the other hand, using a multiple-category ITD significantly increases the spatial scaling exponent for the mean total deformation rate (Figure 8), apparently because there are more LKFs that can be initiated in the thicker pack ice (Hutter et al., 2022a). We note however that the temporal scaling exponent remains unchanged by the introduction of multiple categories in the ITD (Figure 12), suggesting that the local sub-grid redistribution of ice in the ITD that can initiate the formation of new LKFs does not affect the long-range temporal correlations giving rise to the temporal scaling, or at least that the temporal effects of this process are not resolved at the 3-day scale. We also note that the spatial scaling exponents for both runs are more similar for larger moments q (Figure 14). This indicates that the multiple-category ITD mostly increases the spatial localization of smaller deformation rates. We, therefore, hypothesize that the effects of the ITD on the deformation statistics might be more important at lower resolution since strain rates are smaller to start with, but this remains to be verified.

We finally note that the use of an ITD in itself does not guarantee a better spatio-temporal localization of deformations. For instance, the HYCOM-CICE (FSU) runs have a five-category ITD, but the localization of the simulated deformation fields remains low compared to other high-resolution runs. In this case, too few EVP subcycles and large residual errors on the solution may again partly explain the poor localization of deformations.

### 4.2.3. Atmospheric Influence

Kwok (2001) showed that LKF patterns in the observed RGPS deformation fields can remain very similar for long periods of time (~months) suggesting that pack ice deformations occur independently of variability in the wind patterns. However, the majority of LKFs are active on much shorter time scales and LKF lifetimes show an exponential tail (Hutter et al., 2019). Thus, one can question the importance of the atmosphere in setting the observed and simulated small-scale deformation statistics. On the one hand, the majority of the energy input that sets the ice cover in motion originates from the atmospheric wind stress (e.g., Bouchat & Tremblay, 2014; Steele et al., 1997). It could therefore be expected that the simulated deformation scaling statistics are inherited from the turbulent/multi-fractal scaling properties of the atmosphere (e.g., Schmitt et al., 1994). For example, Hutter (2015) showed that the spatial scaling exponent in idealized numerical experiments depends on the spatial resolution of the reanalysis wind forcing, suggesting that the simulated small-scale deformation statistics are, in part, limited by the representation of the atmosphere-ice interactions. On the other hand, the observed scaling properties of sea-ice deformations were shown to hold down to temporal scales much smaller than the atmospheric mesoscale or synoptic temporal scale using ship-based radar observations (Oikkonen et al., 2017). Weiss (2017) suggests this to be a confirmation that the mechanical response of the ice cover is not controlled by the atmosphere, at least not at the mesoscale or synoptic temporal scale.
Here, we note that the degree of temporal multi-fractality and heterogeneity for turbulent wind (i.e., $\mu = 1.45 \pm 0.1$, $C_1 = 0.25 \pm 0.1$; Schmitt et al., 1994) is close to that for RGPS deformation rates (see e.g., Figures 15a and 15b). While this does not confirm that the observed multi-fractality of RGPS deformation rates originates from that of the wind stress, it nonetheless shows that we cannot assume a specific lowest cutoff scale for the atmospheric influence, such that observed sea-ice deformation scaling statistics could well be influenced by the atmosphere below the mesoscale and synoptic scale. We further note that the deformation statistics in the fully coupled atmosphere-ice-ocean RASM-WRF (EAP) runs with higher spatial and temporal resolution of the atmospheric conditions are closer to RGPS observations compared to runs with the same model but forced with an atmospheric reanalysis (i.e., RASM-CORE2 (EAP)—see Figures 1, 2, 7, 8, 11 and 12). Larger deformation rates appear in the PDFs (especially in shear, where the PDF difference metric reduces by ~50%), and the spatio-temporal scaling exponents for the mean total deformation rate also increase. However, we cannot firmly attribute these improvements to the increased complexity of the atmosphere-ice-ocean interactions only, since the fully coupled runs also have an increased number of elastic subcycles (i.e., smaller subcycling time step for the same advective time step) which suggest a better numerical convergence of their solution, although this is not directly quantifiable with the numerical implementation of the EAP rheology.

## 5. Discussion

In the previous sections, both plastic and brittle sea-ice rheologies have shown the potential for reproducing the observed RGPS deformation scaling statistics, even if plastic rheologies do not use specific assumptions that were hypothesized to give rise to the observed scaling of sea-ice deformations (e.g., long-range elastic interactions, damage and healing mechanism, etc.—Weiss & Dansereau, 2017). In particular, a non-zero temporal scaling, intermittency, and temporal multi-fractality are observed for practically all sea-ice models, independently of their spatial scaling. It has previously been assumed that the temporal correlations (or a certain form of memory resulting in time clustering of deformations) giving rise to the temporal scaling and intermittency of deformations should be inherent to the imposed sea-ice mechanical behavior (e.g., Hutter et al., 2018; Weiss & Dansereau, 2017). For instance, Weiss and Dansereau (2017) suggested that plastic sea-ice rheologies cannot reproduce temporal scaling because they do not include stress relaxation, such that temporal correlations cannot develop in their deformation fields. Well-defined LKFs in high-resolution models could also provide such a “memory” via local weakening and divergence of the ice along LKFs (Hutter et al., 2018). Here, we show that plastic sea-ice rheologies, even those without well-defined LKFs, do reproduce a strong temporal localization of deformations and a degree of temporal multi-fractality and intermittency similar to that of the observed RGPS deformation fields. The origin of the multi-fractal temporal scaling in both observed and simulated deformation fields remains to be identified. We note however that we find no significant correlation between simulated temporal scaling exponents and LKFs growth rates or lifetimes (not shown). We hypothesize that temporal correlations in the simulated deformation field could emerge from persistent synoptic atmospheric wind forcing at the basin scale, loading the ice and re-opening recently frozen leads, keeping the ice pack active for several days at a time followed by periods of rest. This is in agreement with RGPS observations showing that deformation of the multiyear ice pack is accommodated by long-lasting LKFs (e.g., Coon et al., 2007). This hypothesis remains to be tested in future work.

We also note that the role of the ocean in setting the observed and simulated deformation statistics has not yet been fully evaluated, even though eddies (e.g., Cassianides et al., 2021), tides (e.g., Heil et al., 2008), and ocean circulation patterns (e.g., Wang et al., 2021; Willmes & Heinemann, 2016) are known to impact sea-ice dynamics. A major difficulty in assessing the role of varying ocean model complexity in coupled ice-ocean simulations resides in the fact that the spatial (and temporal) resolution of the sea-ice model is closely tied, if not the same, to that of the ocean model. In this case, an increased ocean complexity (e.g., eddies) is only resolved together with an increased complexity in sea-ice dynamics (e.g., heterogeneity of the ice strength), such that it becomes difficult to unambiguously separate the effects of refining the sea-ice model resolution from the effects of using a higher-resolution oceanic model without studying the ocean fields in more detail. This is outside the scope of the present study and we, therefore, refrain from providing a quantitative analysis related to the oceanic influence. We note however that the dynamical impacts of eddies on the sea ice is generally limited to the marginal ice zone, where low ice concentration and a highly fractured ice pack allows for a reduced importance of the rheology in the sea-ice dynamics (e.g., Manucharyan & Thompson, 2017). In this sense, we hypothesize that
resolving ocean eddies (or not) would not lead to major differences in the simulated sea-ice deformation statistics evaluated in the Central Arctic Ocean consolidated winter pack ice, such as presented in SIREx. On the other hand, the ocean configurations in current sea-ice models vary widely in their choice of other numerical and physical parameterizations (e.g., choice of ocean coordinate system, mixing schemes, inclusion of tides, etc.) and the effects of these parameterizations on the sea-ice deformation fields also remain to be evaluated. From our experience, however, we argue that the effects on the deformation fields should be small as most of the ocean parameterizations affect thermodynamics more than dynamics. An exception may be the ice-ocean drag param- eterization, because ice-ocean drag removes a large portion of the energy input by the atmosphere (Bouchat & Tremblay, 2014), but this happens on different (larger) spatial scales than the deformation. The stress divergence (rheology) term may still be more important on short spatial and temporal scales. A detailed evaluation of the role of the ocean on the simulated sea-ice deformation statistics using dedicated numerical experiments is undoubt- edly needed in future studies.

In light of the results presented in this first part of the SIREx analysis, a few recommendations for model devel-
opment and implementation emerge for improving the representation of sea-ice deformation statistics by sea-ice
models. First, a spatial resolution of Eulerian models higher than that of the observations is required in order
to better localize the deformations and capture their heterogeneity at the observation scale. In Eulerian models,
several grid cells are always required to represent a velocity discontinuity (e.g., a lead opening or a shear fracture
line). Specifically, in VP finite-difference models, the number of grid points required to resolve a discontinuity
forming under the same forcing conditions remains approximately the same with increasing model spatial reso-
lution (5-7 grid points; Williams & Tremblay, 2018), leading to a spatial refinement of LKFs and an increased
spatial localization of deformations with increasing resolution. The spatial resolution of Eulerian models should
therefore be at least 5-7 times that of the observations for a fair comparison of their deformation field. Note that
as the spatial resolution increases (Δx ≤ 100 km), the continuum assumption (requiring the presence of a large
number of ice floes within one grid cell) is technically no longer valid. However, current sea-ice models remain
able to capture the observed deformation statistics because the simulated deformations are shown here to be
scale-independent.

We further note that it is not expected that models (Eulerian or Lagrangian) reproduce the observed deformation statistics when run at the same nominal scale as the RGPS observation scale. The observed Lagrangian deformation fields are obtained from the motion of tracers at a 10 km spatial scale, but displacement at this scale is closely tied to processes acting on much finer scales that can act as initiation for larger-scale deformations (e.g., micro-fractures, thermal cracking and bending, etc.). These fine-scale processes are sub-grid-scale processes and are usually not resolved or parameterized by sea-ice models, with the exception of neXtSIM which uses a damage parameterization that can represent sub-grid brittle fracturing to some extent. Developing and including well-tuned parameterizations of sub-grid-scale mechanical processes could help with the representation of larger-scale sea-ice deformations in future model developments. For example, the inclusion of a multi-category ice thickness distribution improves the simulated deformation scaling statistics and can also partly improve the LKFs statistics (see also Hutter et al., 2022a). A sub-grid fracturing parameterization (e.g., using a damage formulation) can also be implemented in plastic rheologies, which will help better understand its role on the simulated deformation statistics. New developments in discrete-element models (or hybrid discrete-continuum) and applications of Machine Learning to sea-ice observations may also provide an interesting avenue for the definition and calibration of new sub-grid parameterizations for continuum sea-ice models.

Second, calibrating the yield curve parameters proves to be an efficient solution to improve the deformation statistics, even if sea-ice models are not run at very-high resolution or do not include sub-grid scale mechanical parameterizations. Specifically, we find that increasing the ratio of shear-to-compressive strength provides a better agreement with observed RGPS deformation statistics for both the VP and EVP rheologies. We provided here a new quantitative metric, the sum of the absolute difference of PDFs in logarithmic scale, that is useful for such a calibration of the yield curve parameters. Development and implementation of new yield curves or flow rules (e.g., non-normal flow rule, or adding a time-dependent granular dilatancy—Tremblay & Mysak, 1997; Ringeisen et al., 2019, 2021) may also allow for a better representation of the observed distribution of deformation rates. The spatio-temporal scaling exponents of the mean total deformation rates could also be used for further calibrating new or existing yield curves, however, the usefulness of the scaling of higher moments of the
deformation distribution (i.e., the structure functions) is not clear since the multi-fractality assumption is not
robust for all years in the RGPS records.

Third and finally, ensuring a numerically converged solution without remaining noise appears to be critical for
the small-scale deformation statistics when using an explicit numerical solver such as originally designed in the
EVP and EAP rheologies, although this could not be directly assessed with the available runs. Results nonethe-
less suggest that using an increased number of iterations in the numerical solver along with a small dynamical
time step (i.e., reducing the subcycling time step) improves the EVP deformation scaling statistics. The impact of
the numerical convergence and noise with the EVP solver (but also in VP, EAP, and MEB rheologies), and the
impact of using the modified or adaptative EVP solvers (i.e., mEVP or aEVP—Bouillon et al., 2013; Kimmritz
et al., 2016; Lemieux et al., 2012) remains to be further evaluated. We also note that recent findings suggest that
details in the numerical implementation of sea-ice models (e.g., grid discretization) can have significant impacts
on LKFs patterns in idealized simulations (Mehlmann et al., 2021). The impacts of using different numerical
schemes and numerical implementation of sea-ice models on the deformation statistics in the context of pan-Arc-
tic simulations also needs to be further evaluated and considered for future model developments.

**6. Concluding Remarks**

The first part of the SIREx, with a total of 11 different models, 35 simulations, three different sea-ice rheologies ((E)VP, EAP, and MEB) and a wide range of other model parameterizations, allowed us to investigate how different sea-ice representations affect the deformation statistics using existing and new deformation metrics, namely, the sum of the absolute difference of observed and simulated PDFs of deformation rates, the spatio-temporal scaling exponents, and the multi-fractal parameters describing the structure functions. It is found that the sea-ice rheology, as well as the model configuration (e.g., resolution, atmospheric coupling, numerical convergence, etc.) and physical parameterizations (e.g., ITD and ice strength parameters) can affect the deformation statistics to a similar extent. For this reason, we argue that the aforementioned deformation metrics do not only evaluate the effect of the sea-ice rheology, and that it is important to analyze both the effects of the model configuration or parameterizations along with the effects of the rheological parameters in order to discuss the appropriateness of a given sea-ice rheology in terms of deformation statistics.

We find that a power law scaling and multi-fractality of deformations in both space and time can be achieved
by all sea-ice rheologies evaluated in this study, showing that these metrics are not sufficient to favor the use of
a given rheology, and closing the debate on whether plastic rheologies can reproduce the observed deformation
properties. However, the VP/EVP rheologies implemented in an Eulerian framework need to be run at higher
resolution than that of the observations to yield spatial scaling exponents as high as those observed, because 5-7
grid cells are necessary to spatially resolve discontinuities with such a numerical scheme. It is also expected that
spatial scaling exponents in agreement with the RGPS distribution could be obtained with the EAP rheology at
very-high spatial resolution, given that its spatial scaling exponents are on the same order as for VP/EVP simu-
lations at high-resolution. On the other hand, the spatial localization of MEB (brittle) simulations is larger than
for the plastic rheologies when run at the same resolution as observations. Since these simulations (neXtSIM) are
performed on a Lagrangian mesh that can better localize and follow discontinuities, it is not clear if the higher
spatial scaling exponents are attributable only to the difference in sea-ice rheology.

Interestingly, a strong temporal scaling is better resolved by all rheologies compared to the spatial scaling, inde-
pendently of the models' temporal resolution. While the origin of the observed and simulated temporal scal-
ing remains to be identified, this confirms that there is not only one set of specific rheological assumptions
that can give rise to strong temporal correlations in the deformation fields. We further note that increasing the
shear-to-compressive strength ratio of the ice in elliptical plastic rheologies significantly increases the scaling
exponents, while the addition of multiple ice categories in the ITD does not have a large influence on the temporal
scaling. Coupling the ice model with an atmospheric model instead of forcing with a reanalysis also appears to
significantly affect the temporal (and spatial) multi-fractal parameters and scaling. However, due to a different
number of elastic subcycles in the runs with these variations (likely leading to a difference in numerical conver-
gence of the solution), we cannot firmly attribute this only to a change in the atmospheric forcing/coupling
resolution.
The present study also allowed us to evaluate the usefulness of the scaling metrics to discriminate between differ-
ent sea-ice models, as per SIREx's goal. First, we showed that the decay exponent of the tail of the deformation
PDFs does not efficiently characterize departure from reference PDFs and cannot be used to extract information
on the agreement of the simulated PDFs with observations. We, therefore, introduced a new quantitative metric
that evaluates the sum of the bin-wise absolute differences between the observed and simulated PDFs in loga-
rithmic scale. This metric better characterizes the ability of models to reproduce deformations as large as in
RGPS observations since the logarithmic scale puts more weight on differences in the tail of the PDFs. We also
showed that PDFs of both positive and negative divergence should be preferred to PDFs of absolute divergence
when evaluating model deformation fields, since errors in the negative and positive divergence distributions
may compensate each other to show a (misleading) good agreement with observations in the absolute diver-
gence PDFs. Second, we showed that the spatio-temporal scaling of the mean total deformation rates as usually
implemented does not capture differences in localization of deformations when the density of LKFs also changes
between different simulations. For example, simulated deformation fields with few, but highly localized LKFs
return similarly low scaling exponents as more diffuse deformation fields. We showed that using the signal-to-
noise ratios as weights in the scaling analysis (as introduced by Bouchat & Tremblay, 2020) helps to distinguish
both cases and improves the interpretation of the scaling exponents as a measure of localization of deformations.
This also allows the space-time coupling of the scaling exponents for the mean (q = 1) total deformation rates
to emerge in RGPS observations (Bouchat & Tremblay, 2020) and to be used as an additional metric to evaluate
the simulated deformation fields. Third, we found that the degree of multi-fractality for observed and simulated
deformation fields is generally not quadratic as previously assumed, and that the multi-fractality hypothesis is not
robust for all years of the RGPS records. Our results also show that multi-fractality in both space and time can be
achieved without assuming specific “cascade-like” models for the deformation of the sea-ice cover, which leaves
open the question of what physical/mechanical parameterizations common to all the tested sea-ice models are
critical in producing the multi-fractality. In this sense, it is unclear whether the multi-fractal analysis is appropri-
ate to calibrate or evaluate sea-ice rheologies, since the observed deformation multi-fractality could emerge from
parameterizations other than the rheology (e.g., atmospheric turbulent momentum transfer).

Keeping in mind that the MEB and EAP rheologies are under-represented in the participating sea-ice models, the conclusions presented here should be tested using a larger number of experiments including more MEB and EAP runs, or ideally, by running a unique model configuration with different sea-ice rheologies. Specifically, to eliminate the potential differences associated with using a Lagrangian mesh, the deformation statistics of MEB runs implemented on a Eulerian grid (as recently done by Plante et al., 2020) should be evaluated. Nevertheless, this study shows that the (E)VP rheology—used in a majority of climate models—does generate large deformation rates that are highly localized in space and time, albeit by using a higher spatial resolution than currently used in GCMs and CMIP-type climate models. Generating large, localized deformation rates is a necessary condition for sea-ice models to achieve before their effect on the Arctic climate system can be assessed. While a thorough study of the impacts of sea-ice deformations and rheology in Global Climate Model runs remains to be performed, the analysis of LKFs statistics (and their link to ice thickness and concentration anomalies) presented in the second part of the SIREx analysis offers a complementary step to the present analysis toward improving the representation of sea ice in climate projections.

**Appendix A: List of Model, Configuration, and Reanalysis Acronyms**

Table A1 details the meaning of most acronyms relating to models, configurations, and reanalyses used in the
runs participating in this study.
<caption>Table A1</caption>

Definition of Acronyms Relating to Models, Configurations, and Reanalyses Found in Table 1

<table><thead><tr><th>Acronym</th><th>Meaning</th></tr></thead><tbody><tr><td colspan="2"><strong>Models</strong></td></tr><tr><td>MITgcm</td><td>Massachusetts Institute of Technology General Circulation Model</td></tr><tr><td>McGill-SIM</td><td>McGill University Sea Ice Model</td></tr><tr><td>NEMO</td><td>Nucleus for European Modeling of the Ocean Model</td></tr><tr><td>LIM3</td><td>Louvain-la-Neuve Sea Ice Model v3</td></tr><tr><td>LIM2</td><td>Louvain-la-Neuve Sea Ice Model v2</td></tr><tr><td>HYCOM</td><td>Hybrid Coordinate Ocean Model</td></tr><tr><td>CICE4</td><td>Los Alamos Sea Ice Model v4</td></tr><tr><td>RIOPS</td><td>Regional Ice Ocean Prediction System Model</td></tr><tr><td>FESOM</td><td>Finite-Element/volume Sea ice-Ocean Model</td></tr><tr><td>RASM</td><td>Regional Arctic System Model</td></tr><tr><td>WRF</td><td>Weather Research and Forecasting Model</td></tr><tr><td colspan="2"><strong>Configurations</strong></td></tr><tr><td>CREG4</td><td>CONCEPTS Regional configuration—1/4°</td></tr><tr><td>CREG12</td><td>CONCEPTS Regional configuration—1/12°</td></tr><tr><td>ANHA4</td><td>Arctic and Northern Hemisphere Atlantic configuration—1/4°</td></tr><tr><td>ANHA12</td><td>Arctic and Northern Hemisphere Atlantic configuration—1/12°</td></tr><tr><td colspan="2"><strong>Reanalyses</strong></td></tr><tr><td>JRA55</td><td>Japanese 55 yr Reanalysis</td></tr><tr><td>ERA-Interim</td><td>ECMWF Re-Analysis, third-generation</td></tr><tr><td>NCEP/NCAR</td><td>National Centers for Environmental Prediction/National Center for Atmospheric Research</td></tr><tr><td>DFS</td><td>Drakkar Forcing Sets</td></tr><tr><td>CFSR/CFSv2</td><td>Climate Forecast System Reanalysis v2</td></tr><tr><td>CGRF</td><td>Canadian Meteorological Centre's Global Deterministic Prediction System Reforecasts</td></tr><tr><td>CORE</td><td>Common Ocean–Ice Reference Experiment</td></tr><tr><td>CORE2</td><td>Common Ocean–Ice Reference Experiment v2</td></tr></tbody></table>

## Appendix B: Strain Rate Error Estimation

Trajectory errors and boundary-definition errors affect both the observed and simulated Lagrangian deformation estimates. Following Bouchat and Tremblay (2020), we consider only the trajectory errors to compute the signal-to-noise ratio of the deformation estimates and use this ratio as weight when averaging the deformation distribution for the scaling analysis. Trajectory errors result from uncertainty on the Lagrangian position of the cell vertices used to compute the strain rates ($u_x, u_y, v_x, v_y$). When using the line integral approximations of Equations 1–4 to evaluate the strain rates between time $t$ and $t + \Delta t$, the trajectory error on the strain rates can be approximated using the propagation of uncertainty as in Bouchat and Tremblay (2020):

$$ \sigma_{u_x}^2 = u_x^2 \left( \frac{\sigma_A^2}{A^2} \right) + \sum_{k=1}^{4} \left( \frac{(y_{k+1} - y_{k-1})^2}{4A^2 T^2} \right) (\sigma_x^2 + \sigma_{x'}^2) + \sum_{k=1}^{4} \left( \frac{(u_{k-1} - u_{k+1})^2}{4A^2} \right) \sigma_x^2, \quad (B1) $$

$$ \sigma_{u_y}^2 = u_y^2 \left( \frac{\sigma_A^2}{A^2} \right) + \sum_{k=1}^{4} \left( \frac{(x_{k+1} - x_{k-1})^2}{4A^2 T^2} \right) (\sigma_x^2 + \sigma_{x'}^2) + \sum_{k=1}^{4} \left( \frac{(u_{k-1} - u_{k+1})^2}{4A^2} \right) \sigma_x^2, \quad (B2) $$

$$ \sigma_{v_x}^2 = v_x^2 \left( \frac{\sigma_A^2}{A^2} \right) + \sum_{k=1}^{4} \left( \frac{(y_{k+1} - y_{k-1})^2}{4A^2 T^2} \right) (\sigma_x^2 + \sigma_{x'}^2) + \sum_{k=1}^{4} \left( \frac{(v_{k-1} - v_{k+1})^2}{4A^2} \right) \sigma_x^2, \quad (B3) $$
$$ \sigma_{v_y}^2 = v_y^2 \left( \frac{\sigma_A^2}{A^2} \right) + \sum_{k=1}^{4} \left( \frac{(x_{k+1} - x_{k-1})^2}{4A^2T^2} \right) (\sigma_x^2 + \sigma_{x'}^2) + \sum_{k=1}^{4} \left( \frac{(v_{k-1} - v_{k+1})^2}{4A^2} \right) \sigma_x^2 \quad (B4) $$

and $\sigma_A$, the error on the cell area A at time t, is also derived from the propagation of uncertainty as (e.g., Lindsay & Stern, 2003):

$$ \sigma_A^2 = \frac{1}{4} \sum_{k=1}^{4} \left[ (y_{k+1} - y_{k-1})^2 + (x_{k-1} - x_{k+1})^2 \right] \sigma_x^2 \quad (B5) $$

where $\sigma_x$ and $\sigma_{x'}$ are the position errors at time t and $t + \Delta t$ respectively, and $(x_k, y_k)$, $(u_k, v_k)$ are the position and velocity of the cell vertex k at time t. The error on the total deformation rate ($\dot{e}_{tot}$) and strain rate invariants ($\dot{e}_I$, $\dot{e}_{II}$) is also derived from the propagation of uncertainty and is proportional to Equations B1–B4 (see e.g., Bouchat & Tremblay, 2020). Note that we have ignored timing uncertainties (i.e., $\sigma_t$) in Equations B1–B5 above.

The position errors $\sigma_x$ and $\sigma_{x'}$ contributing to Equations B1–B5 can originate from (a) geolocation errors ($\sigma_{geo}$) that are due to uncertainty of the recording instrument or acquisition method, and/or (b) tracking errors ($\sigma_{track}$) that occur when the position of tracked features on images are misidentified at the pixel level on satellite images. For RGPS strain rates derived from the tracking of ice features in consecutive SAR images, we can assume that the geolocation error is zero and that the position of a tracked feature on the first SAR image at time t is always known exactly, that is, $\sigma_x^{\text{RGPS}} = \sigma_{\text{geo}}^{\text{RGPS}} = 0$ (see e.g., Bouchat & Tremblay, 2020; Dierking et al., 2020). The position of that feature on the second image at time $t + \Delta t$ is however affected by a tracking error of one pixel in the SAR images, that is, $\sigma_{x'}^{\text{RGPS}} = \sigma_{\text{track}}^{\text{RGPS}} = 100 \text{ m}$ (Lindsay & Stern, 2003). In this case, $\sigma_A^{\text{RGPS}} = 0$ and the uncertainty on the RGPS strain rates is reduced to a single summation, for example:

$$ \sigma_{u_x}^2 = \sum_{k=1}^{4} \left( \frac{(y_{k+1} - y_{k-1})^2}{4A^2T^2} \right) \sigma_{\text{track}}^2 \quad (B6) $$

For the reconstructed model Lagrangian trajectories, no tracking is done, but tracers are instead advected using the model velocity fields. Model tracking errors are therefore zero (i.e., $\sigma_{\text{track}}^{\text{model}} = 0$), but geolocation errors accumulate in time with every Lagrangian integration step n due to uncertainty on the model velocity fields, such that $\sigma_{x_n}^{\text{model}} = \sigma_{\text{geo},n}^{\text{model}} \neq 0$. To see this, consider the case where the initial position of a tracer at time $t_0 = 0$ is known perfectly (i.e., $\sigma_{x_0}^{\text{model}} = \sigma_{\text{geo},0}^{\text{model}} = 0$). At time $t_1 = \Delta t$, the position of the advected tracer is $x_1 = x_0 + U_0 \Delta t$, where $U_0$ is the model velocity in the x-direction at $t_0$. At time $t_2 = 2\Delta t$, the position is $x_2 = x_1 + U_1 \Delta t = x_0 + (U_0 + U_1) \Delta t$ and similarly, at any number n of subsequent integration steps $\Delta t$, we have:

$$ x_n = x_0 + (U_0 + U_1 + \dots + U_{n-1}) \Delta t. \quad (B7) $$

Using the propagation of uncertainty and again neglecting timing uncertainties, the uncertainty $\sigma_{x_n}^{\text{model}}$ on the model Lagrangian trajectory position at time $t_n$, is therefore given by:

$$ \sigma_{x_n}^{\text{model}} = \sigma_{\text{geo},n}^{\text{model}} = \sqrt{n} \sigma_U \Delta t, \quad (B8) $$

where we assume that the uncertainty on the model velocity remains the same in time (i.e., $\sigma_{U_n} = \sigma_U$ for all n). The error on the model Lagrangian trajectory positions therefore grows with the square-root of the number of integration steps. Assuming that the error on the model velocity in the y-direction is the same as in the x-direction, it is also straightforward to show that $\sigma_{y_n}^{\text{model}} = \sigma_{x_n}^{\text{model}}$.

Here, we conservatively assume that all the points on the model trajectories have the largest error possible, that is, the error of the last point after the full integration is done. Other than for simplicity in our calculations, this also ensures that earlier deformation estimates (e.g., in January) will not be weighted more in the model scaling analysis than later ones (e.g., in March). This allows for a better comparison with the RGPS signal-to-noise weighted scaling analysis, given that the uncertainty on the RGPS deformation estimates does not explicitly depend on time (see e.g., Equation B6). For $\Delta t = 1 \text{ hr time steps}$ and a 90-day integration period (i.e., $n = 2,160$ steps—1 January to 30 March), we therefore fix $\sqrt{n}\Delta t = 1.7 \times 10^5$ seconds, and for any point along the trajectory we have:

$$ \sigma_x^{\text{model}} = \sigma_{x'}^{\text{model}} = \sigma_{x_{2160}}^{\text{model}} = (1.7 \times 10^5) \sigma_U. \quad (B9) $$
Acknowledgments

The authors thank the anonymous
reviewers for their dedicated review
and constructive comments that helped
improve both SIREx manuscripts. D.
Dukhovskoy was funded by the DOE
(award DE-SC0014378) and HYCOM
NOPP (award N00014-19-1-2674).
The HYCOM-CICE simulations were
supported by a grant of computer
time from the DoD High-Performance
Computing Modernization Program at
NRL SSC. The daily fields from the 0.08
HYCOM-CICE experiment are available
at the HYCOM data server ftp://ftp.
hycom.org/datasets/ARCc0.08/expt_11.0/
data/. T. Rasmussen was funded by the
Danish State through the National centre
for Climate Research and by the SALI-
ENSEAS project part of the ERA4CS
program, which is co-funded by the Inno-
vation Fund Denmark and the Horizon
2020 Framework Programme of the Euro-
pean Union (Grant 690462). P. Myers
was funded by the Natural Sciences and
Engineering Research Council (NSERC)
of Canada (RGPIN 04357 and RGPCC
433898). Experiments were run and are
archived using facilities provided by
Compute Canada (www.computecanada.
ca). For more details on the ANHA
configuration, visit http://knossos.eas.
ualberta.ca/anha/anhatable.php. Q.
Wang was supported by the German
Helmholtz Climate Initiative REKLIM
(Regional Climate Change). B. Tremblay
was funded by the Natural Science and
Engineering and Research Council
(NSERC) Discovery Program and by
the Environment and Climate Change
Canada Grants & Contributions program.
This work is also a contribution to the
research program of Québec-Océan. Y.
Lee and W. Maslowski were supported
by the grants from the U.S. Department
of Energy Regional and Global Model
Analysis program, the Office of Naval
Research Arctic and Global Prediction
program, and the National Science Foun-
dation Arctic System Science program.
The computer resources for the RASM
simulations were provided by the Frontier
project from the Department of Defense
High-Performance Computing Modern-
ization Program. Open access funding
enabled and organized by Projekt DEAL.

The error on the ice velocity σ<sub>U</sub> is due to an interpolation error of gridded model velocity fields to the trajectory positions, as well as to the numerical error on the dynamical solution resolved by the sea-ice models. The latter source of error depends on the model time step and spatial resolution, the choice of numerical solver and number of iterations performed to solve the nonlinear dynamical equations (i.e., convergence of the solution), on the numerical regularization methods and parameterization schemes used, etc (e.g., Bouillon et al., 2013; Kimmritz et al., 2015, 2017; Lemieux et al., 2012, 2010, 2008; Plante et al., 2020). The values of σ<sub>U</sub> are therefore expected to vary within the participating simulations, however, those values are unknown and a complete convergence/ error analysis is outside the scope of the present study. We, therefore, assume an upper bound of σ<sub>U</sub> = 0.006 m/s for all simulations regardless of their specific configurations and parameterizations, which corresponds to a typi- cal velocity error for high-resolution EVP simulations with a default number (120) of elastic subcycles and a time step of 20 min (Lemieux et al., 2012), and should also largely encompass the interpolation error. This corresponds to a position uncertainty of σ<sub>x</sub><sup>model</sup> = σ<sub>x'</sub><sup>model</sup> ≈ 1000 m for simulated Lagrangian trajectories. Note that this error is especially overestimated for very-high resolution models which generally have a much smaller time step and a larger number of elastic subcycles.

Data Availability Statement

All published code and data products originating from the SIREx is organized in the SIREx data collection: https://zenodo.org/communities/sirex/: this includes the model output in netCDF format (Bouchat et al., 2022a), the model Lagrangian trajectories and deformation data (Bouchat et al., 2022b), and the RGPS Lagrangian composite deformation data (Bouchat and Hutter, 2022), as well as the LKF data (Hutter et al., 2022b).

References

Benzi, R., Paladin, G., Parisi, G., & Vulpiani, A. (1984). On the multifractal nature of fully developed turbulence and chaotic systems. *Journal of Physics A: Mathematical and General*, 17(18), 3521-3531. https://doi.org/10.1088/0305-4470/17/18/021

Bouchat, A., Hutter, N., Chanut, J., Dupont, F., Dukhovskoy, D., Garric, G., et al. (2022a). Sea Ice Rheology Experiment (SIREx)-Model output data (Version 1.0) [Data set]. Zenodo. https://doi.org/10.5281/ZENODO.5555329

Bouchat, A., Hutter, N., Chanut, J., Dupont, F., Dukhovskoy, D., Garric, G., et al. (2022b). Model Lagrangian trajectories and deformation data analyzed in the Sea Ice Rheology Experiment—Part I (Version 1.0) [Data set]. Zenodo. https://doi.org/10.5281/zenodo.6321323

Bouchat, A., & Hutter, N. (2022). Deformation composite of the RADARSAT Geophysical Processor System (RGPS) Lagrangian motion data (1.0) [Data set]. Zenodo. https://doi.org/10.5281/zenodo.6321327

Bouchat, A., & Tremblay, B. (2014). Energy dissipation in viscous-plastic sea-ice models. *Journal of Geophysical Research: Oceans*, 119(2), 976-994. https://doi.org/10.1002/2013JC009436

Bouchat, A., & Tremblay, B. (2017). Using sea-ice deformation fields to constrain the mechanical strength parameters of geophysical sea ice. *Journal of Geophysical Research: Oceans*, 122, 5802–5825. https://doi.org/10.1002/2017JC013020

Bouchat, A., & Tremblay, B. (2020). Reassessing the quality of sea-ice deformation estimates derived from the RADARSAT geophysical processor system and its impact on the spatio-temporal scaling statistics. *Journal of Geophysical Research: Oceans*, 125(8), e2019JC015944. https://doi.org/10.1029/2019JC015944

Bouillon, S., Fichefet, T., Legat, V., & Madec, G. (2013). The elastic-viscous-plastic method revisited. *Ocean Modeling*, 71, 2–12. https://doi.org/10.1016/j.ocemod.2013.05.013

Bouillon, S., & Rampal, P. (2015a). On producing sea ice deformation data sets from SAR-derived sea ice motion. *The Cryosphere*, 9(2), 663-673. https://doi.org/10.5194/tc-9-663-2015

Bouillon, S., & Rampal, P. (2015b). Presentation of the dynamical core of neXtSIM, a new sea ice model. *Ocean Modeling*, 91, 23–37. https://doi.org/10.1016/j.ocemod.2015.04.005

Cassianides, A., Lique, C., & Korosov, A. (2021). Ocean eddy signature on sar-derived sea ice drift and vorticity. *Geophysical Research Letters*, 48(6). e2020GL092066. https://doi.org/10.1029/2020GL092066

Ceniceros, H. D., & Hou, T. Y. (2001). An efficient dynamically adaptive mesh for potentially singular solutions. *Journal of Computational Physics*, 172(2), 609–639. https://doi.org/10.1006/jcph.2001.6844

Clauset, A., Shalizi, C. R., & Newman, M. E. J. (2009). Power-law distributions in empirical data. *SIAM Review*, 51(4), 661–703. https://doi.org/10.1137/070710111

Coon, M., Kwok, R., Levy, G., Pruis, M., Schreyer, H., & Sulsky, D. (2007). Arctic Ice Dynamics Joint Experiment (AIDJEX) assumptions revisited and found inadequate. *Journal of Geophysical Research: Oceans*, 112(C11). https://doi.org/10.1029/2005JC003393

Courtois, P., Hu, X., Pennelly, C., Spence, P., & Myers, P. G. (2017). Mixed layer depth calculation in deep convection regions in ocean numerical models. *Ocean Modeling*, 120, 60–78. https://doi.org/10.1016/j.ocemod.2017.10.007

Dansereau, V., Weiss, J., Saramito, P., & Lattes, P. (2016). A Maxwell elasto-brittle rheology for sea ice modeling. *The Cryosphere*, 10(3), 1339–1359. https://doi.org/10.5194/tc-10-1339-2016

Dierking, W., Stern, H. L., & Hutchings, J. K. (2020). Estimating statistical errors in retrievals of ice velocity and deformation parameters from satellite images and buoy arrays. *The Cryosphere*, 14(9), 2999–3016. https://doi.org/10.5194/tc-14-2999-2020

Dukhovskoy, D. S., Yashayaev, I., Proshutinsky, A., Bamber, J. L., Bashmachnikov, I. L., Chassignet, E. P., et al. (2019). Role of Greenland freshwater anomaly in the recent freshening of the subpolar North Atlantic. *Journal of Geophysical Research: Oceans*, 124(5), 3333–3360. https://doi.org/10.1029/2018JC014686
Dumont, D., Gratton, Y., & Arbetter, T. E. (2009). Modeling the dynamics of the North water polynya ice bridge. *Journal of Physical Oceanography*, 39(6), 1448–1461. https://doi.org/10.1175/2008JPO3965.1

Dupont, F., Higginson, S., Bourdallé-Badie, R., Lu, Y., Roy, F., Smith, G. C., et al. (2015). A high-resolution ocean and sea-ice modeling system for the Arctic and the North Atlantic oceans. *Geoscientific Model Development*, 8, 1577–1594. https://doi.org/10.5194/gmd-8-1577-2015

Flato, G. M., & Hibler, W. D. (1992). Modeling pack ice as a cavitating fluid. *Journal of Physical Oceanography*, 22(6), 626–651. https://doi.org/10.1175/1520-0485(1992)022<0626:mpiaac>2.0.co;2

Girard, L., Bouillon, S., Weiss, J., Amitrano, D., Fichefet, T., & Legat, V. (2011). A new modeling framework for sea-ice mechanics based on elasto-brittle rheology. *Anals of Glaciology*, 52(57), 123–132. https://doi.org/10.3189/172756411795931499

Girard, L., Weiss, J., Molines, J. M., Barnier, B., & Bouillon, S. (2009). Evaluation of high-resolution sea ice models on the basis of statistical and scaling properties of Arctic sea ice drift and deformation. *Journal of Geophysical Research: Oceans*, 114(C8). https://doi.org/10.1029/2008JC005182

Heil, P., Hutchings, J. K., Worby, A. P., Johansson, M., Launiainen, J., Haas, C., & Hibler, W. D. (2008). Tidal forcing on sea-ice drift and deformation in the western Weddell sea in early austral summer, 2004. *Deep Sea Research Part II: Topical Studies in Oceanography*, 55(8), 943–962. https://doi.org/10.1016/j.dsr2.2007.12.026

Hibler, W. D. (1979). A dynamic thermodynamic sea ice model. *Journal of Physical Oceanography*, 9, 815–946. https://doi.org/10.1175/1520-0485(1979)009<0815:adtsim>2.0.co;2

Hibler, W. D., & Walsh, J. E. (1982). On modeling seasonal and interannual fluctuations of Arctic sea ice. *Journal of Physical Oceanography*, 12(12), 1514–1523. https://doi.org/10.1175/1520-0485(1982)012<1514:omsaif>2.0.co;2

Hu, X., Sun, J., Chan, T., & Myers, P. G. (2018). Thermodynamic and dynamic ice thickness contributions in the Canadian Arctic Archipelago in NEMO-LIM2 numerical simulations. *The Cryosphere*, 12, 1233–1247. https://doi.org/10.5194/tc-12-1233-2018

Hunke, E. C., & Dukowicz, J. K. (1997). An elastic-viscous-plastic model for sea ice dynamics. *Journal of Physical Oceanography*, 27(9), 1849–1867. https://doi.org/10.1175/1520-0485(1997)027<1849:aevpmf>2.0.co;2

Hunke, E. C., & Dukowicz, J. K. (2002). The elastic-viscous-plastic sea ice dynamics model in general orthogonal curvilinear coordinates on a sphere-incorporation of metric terms. *Monthly Weather Review*, 130, 1848–1865. https://doi.org/10.1175/1520-0493(2002)130<1848:tevpsi>2.0.co;2

Hutter, N. (2015). *Viscous plastic sea ice models at very high-resolution, (M.Sc. thesis)*. University of Bremen, Alfred Wegener Institute. Retrieved from http://epic.awi.de/38870/

Hutter, N., Bouchat, A., Dupont, F., Dukhovskoy, D., Koldunov, N., Lee, Y. J., et al. (2022a). Sea Ice Rheology Experiment (SIREx): 2. Evaluating linear kinematic features in high-resolution sea ice simulations. *Journal of Geophysical Research: Oceans*, 127, e2021JC017666. https://doi.org/10.1029/2021JC017666

Hutter, N., Bouchat, A., Dupont, F., Dukhovskoy, D., Koldunov, N., Lee, Y., et al. (2022b). Linear kinematic feature detected and tracked in sea-ice deformation simulation by all models participating in the Sea Ice Rheology Experiment and from RGPS. *Zenodo*. https://doi.org/10.5281/zenodo.6315226

Hutter, N., & Losch, M. (2020). Feature-based comparison of sea ice deformation in lead-permitting sea ice simulations. *The Cryosphere*, 14(1), 93–113. https://doi.org/10.5194/tc-14-93-2020

Hutter, N., Losch, M., & Menemenlis, D. (2018). Scaling properties of Arctic sea ice deformation in a high-resolution viscous-plastic sea ice model and in satellite observations. *Journal of Geophysical Research: Oceans*, 123(1), 672–687. https://doi.org/10.1002/2017JC013119

Hutter, N., Zampieri, L., & Losch, M. (2019). Leads and ridges in Arctic sea ice from RGPS data and a new tracking algorithm. *The Cryosphere*, 13(2), 627–645. https://doi.org/10.5194/tc-13-627-2019

Ip, C. F., Hibler, W. D. III, & Flato, G. M. (1991). On the effect of rheology on seasonal sea-ice simulations. *Annals of Glaciology*, 15, 17–25. https://doi.org/10.1017/s0260305500009496

Kimmritz, M., Danilov, S., & Losch, M. (2015). On the convergence of the modified elastic-viscous-plastic method for solving the sea ice momentum equation. *Journal of Computational Physics*, 296, 90–100. https://doi.org/10.1016/j.jcp.2015.04.051

Kimmritz, M., Danilov, S., & Losch, M. (2016). The adaptive EVP method for solving the sea ice momentum equation. *Ocean Modeling*, 101, 59–67. https://doi.org/10.1016/j.ocemod.2016.03.004

Kimmritz, M., Losch, M., & Danilov, S. (2017). A comparison of viscous-plastic sea ice solvers with and without replacement pressure. *Ocean Modeling*, 115, 59–69. https://doi.org/10.1016/j.ocemod.2017.05.006

Koldunov, N. V., Danilov, S., Sidorenko, D., Hutter, N., Losch, M., Goessling, H., et al. (2019). Fast EVP solutions in a high-resolution sea ice model. *Journal of Advances in Modeling Earth Systems*, 11(5), 1269–1284. https://doi.org/10.1029/2018MS001485

König Beatty, C., & Holland, D. M. (2010). Modeling landfast sea ice by adding tensile strength. *Journal of Physical Oceanography*, 40(1), 185–198. https://doi.org/10.1175/2009JPO4105.1

Kreyscher, M., Harder, M., Lemke, P., & Flato, G. M. (2000). Results of the Sea Ice Model Intercomparison Project: Evaluation of sea ice rheology schemes for use in climate simulations. *Journal of Geophysical Research: Oceans*, 105(C5). https://doi.org/10.1029/1999JC000016

Kwok, R. (1998). The RADARSAT geophysical processor system. In C. Tsatsoulis, & R. Kwok (Eds.), *Analysis of SAR data of the polar oceans* (pp. 235–257). Springer-Verlag. https://doi.org/10.1007/978-3-642-60282-5_11

Kwok, R. (2001). Deformation of the Arctic Ocean sea ice cover between November 1996 and April 1997: A qualitative survey. In *IUTAM Symposium on scaling laws in ice mechanics and ice dynamics* (pp. 315–322). Springer Netherlands. https://doi.org/10.1007/978-94-015-9735-7_26

Lemieux, J.-F., Dupont, F., Blain, P., Roy, F., Smith, G. C., & Flato, G. M. (2016). Improving the simulation of landfast ice by combining tensile strength and a parameterization for grounded ridges. *Journal of Geophysical Research: Oceans*, 121(10), 7354–7368. https://doi.org/10.1002/2016JC012006

Lemieux, J.-F., Knoll, D. A., Tremblay, B., Holland, D. M., & Losch, M. (2012). A comparison of the Jacobian-free Newton-Krylov method and the EVP model for solving the sea ice momentum equation with a viscous-plastic formulation: A serial algorithm study. *Journal of Computational Physics*, 231(17), 5926–5944. https://doi.org/10.1016/j.jcp.2012.05.024

Lemieux, J.-F., Tremblay, B., Sedláček, J., Tupper, P., Thomas, S., Huard, D., & Auclair, J.-P. (2010). Improving the numerical convergence of viscous-plastic sea ice models with the Jacobian-free Newton-Krylov method. *Journal of Computational Physics*, 229, 2840–2852. https://doi.org/10.1016/j.jcp.2009.12.011

Lemieux, J.-F., Tremblay, B., Thomas, S., Sedláček, J., & Mysak, L. A. (2008). Using the preconditioned Generalized Minimum RESidual (GMRES) method to solve the sea-ice momentum equation. *Journal of Geophysical Research*, 113(C10004). https://doi.org/10.1029/2007JC004680

Lindsay, R. W., & Stern, H. L. (2003). The RADARSAT geophysical processor system: Quality of sea ice trajectory and deformation estimates. *Journal of Atmospheric and Oceanic Technology*, 20(9), 1333–1347. https://doi.org/10.1175/1520-0426(2003)020<1333:trgpsq>2.0.co;2

Lovejoy, S., & Schertzer, D. (1995). How bright is the coast of Brittany? *Fractals in geoscience and remote sensing* (pp. 102–151).
Lovejoy, S., & Schertzer, D. (2007). Scale, scaling, and multifractals in geophysics: Twenty years on. In *Nonlinear dynamics in geosciences* (pp. 311–337). Springer.

Lovejoy, S., & Schertzer, D. (2013). *The weather and climate: Emergent laws and multifractal cascades*. Cambridge University Press.

Madec, G., Bourdallé-Badie, R., Bouttier, P.-A., Bricaud, C., Bruciaferri, D., Calvert, D., et al. (2019). NEMO ocean engine (computer software manual No. 27). Zenodo. https://doi.org/10.5281/zenodo.1464816

Madsen, K. S., Rasmussen, T. A. S., Ribergaard, M. H., & Ringgaard, I. M. (2016). High-resolution sea ice modeling and validation of the Arctic with focus on south Greenlandic waters, 2004–2013. *Polarforschung*, 85. https://doi.org/10.2312/POLFOR.2016.006

Manucharyan, G. E., & Thompson, A. F. (2017). Submesoscale sea ice-ocean interactions in marginal ice zones. *Journal of Geophysical Research: Oceans*, 122(12), 9455–9475. https://doi.org/10.1002/2017JC012895

Marsan, D., Stern, H., Lindsay, R., & Weiss, J. (2004). Scale dependence and localization of the deformation of Arctic sea ice. *Physical Review Letters*, 93, 178501. https://doi.org/10.1103/PhysRevLett.93.178501

Marsan, D., & Weiss, J. (2010). Space/time coupling in brittle deformation at geophysical scales. *Earth and Planetary Science Letters*, 296(3–4), 353–359. https://doi.org/10.1016/j.epsl.2010.05.019

Mehlmann, C., Danilov, S., Losch, M., Lemieux, J. F., Hutter, N., Richter, T., et al. (2021). Simulating linear kinematic features in viscous-plastic sea ice models on quadrilateral and triangular grids with different variable staggering. *Journal of Advances in Modeling Earth Systems*, 13(11). https://doi.org/10.1029/2021ms002523

Miller, P. A., Laxon, S. W., & Feltham, D. L. (2005). Improving the spatial distribution of modeled Arctic sea ice thickness. *Geophysical Research Letters*, 32(18), L18503. https://doi.org/10.1029/2005GL023622

Miller, P. A., Laxon, S. W., Feltham, D. L., & Cresswell, D. J. (2006). Optimization of a sea ice model using basinwide observations of Arctic sea ice thickness, extent, and velocity. *Journal of Climate*, 19(7), 1089–1108. https://doi.org/10.1175/JCLI3648.1

Mohammadi-Aragh, M., Goessling, H. F., Losch, M., Hutter, N., & Jung, T. (2018). Predictability of Arctic sea ice on weather time scales. *Scientific Reports*, 8, 6514. https://doi.org/10.1038/s41598-018-24660-0

Muliwijk, M., Ilicak, M., Cornish, S. B., Danilov, S., Gelderloos, R., Gerdes, R., et al. (2019). Arctic Ocean response to Greenland sea wind anomalies in a suite of model simulations. *Journal of Geophysical Research: Oceans*, 124(8), 6286–6322. https://doi.org/10.1029/2019JC015101

Notz, D., Jahn, A., Holland, M., Hunke, E., Massonnet, F., Stroeve, J., et al. (2016). The CMIP6 Sea-Ice Model Intercomparison Project (SIMIP): Understanding sea ice through climate-model simulations. *Geoscientific Model Development*, 9(9), 3427–3446. https://doi.org/10.5194/gmd-9-3427-2016

Oikkonen, A., Haapala, J., Lensu, M., Karvonen, J., & Itkin, P. (2017). Small-scale sea ice deformation during N-ICE2015: From compact pack ice to marginal ice zone. *Journal of Geophysical Research: Oceans*, 122(6), 5105–5120. https://doi.org/10.1002/2016JC012387

Plante, M., Tremblay, B., Losch, M., & Lemieux, J.-F. (2020). Landfast sea ice material properties derived from ice bridge simulations using the Maxwell elasto-brittle rheology. *The Cryosphere*, 14(6), 2137–2157. https://doi.org/10.5194/tc-14-2137-2020

Rampal, P., Bouillon, S., Olason, E., & Morlighem, M. (2016). neXtSIM: A new Lagrangian sea ice model. *The Cryosphere*, 10(3), 1055–1073. https://doi.org/10.5194/tc-10-1055-2016

Rampal, P., Dansereau, V., Olason, E., Bouillon, S., Williams, T., Korosov, A., & Samaké, A. (2019). On the multi-fractal scaling properties of sea ice deformation. *The Cryosphere*, 13(9), 2457–2474. https://doi.org/10.5194/tc-13-2457-2019

Rampal, P., Weiss, J., & Marsan, D. (2009). Positive trend in the mean speed and deformation rate of Arctic sea ice, 1979–2007. *Journal of Geophysical Research: Oceans*, 114(C5). https://doi.org/10.1029/2008JC005066

Rampal, P., Weiss, J., Marsan, D., Lindsay, R., & Stern, H. (2008). Scaling properties of sea ice deformation from buoy dispersion analysis. *Journal of Geophysical Research: Oceans*, 113(C3). https://doi.org/10.1029/2007JC004143

Ringeisen, D., Losch, M., Tremblay, L. B., & Hutter, N. (2019). Simulating intersection angles between conjugate faults in sea ice with different viscous-plastic rheologies. *The Cryosphere*, 13(4), 1167–1186. https://doi.org/10.5194/tc-13-1167-2019

Ringeisen, D., Tremblay, L. B., & Losch, M. (2021). Non-normal flow rules affect fracture angles in sea ice viscous-plastic rheologies. *The Cryosphere*, 15(6), 2873–2888. https://doi.org/10.5194/tc-15-2873-2021

Rothrock, D. A. (1975). The energetics of the plastic deformation of pack ice by ridging. *Journal of Geophysical Research*, 80(33), 1896–1977. https://doi.org/10.1029/JC080i033p04514

Schmitt, F., Schertzer, D., Lovejoy, S., & Brunet, Y. (1994). Estimation of universal multifractal indices for atmospheric turbulent velocity fields. *In Fractals in natural sciences* (pp. 274–281). https://doi.org/10.1142/9789814503792_0025

Schulson, E. M. (2004). Compressive shear faults within arctic sea ice: Fracture on scales large and small. *Journal of Geophysical Research: Oceans*, 109(C7), C07016. https://doi.org/10.1029/2003JC002108

Sreen, G., Kwok, R., Menemenlis, D., & Nguyen, A. T. (2017). Sea-ice deformation in a coupled ocean-sea-ice model and in satellite remote sensing data. *The Cryosphere*, 11(4), 1553–1573. https://doi.org/10.5194/tc-11-1553-2017

Steele, M., Zhang, J., Rothrock, D., & Stern, H. (1997). The force balance of sea ice in a numerical model of the Arctic Ocean. *Journal of Geophysical Research: Oceans*, 102(C9), 21061–21079. https://doi.org/10.1029/97JC01454

Stern, H. L., & Lindsay, R. W. (2009). Spatial scaling of Arctic sea ice deformation. *Journal of Geophysical Research: Oceans*, 114(C10), C10017. https://doi.org/10.1029/2009JC005380

Stroeve, J., Barrett, A., Serreze, M., & Schweiger, A. (2014). Using records from submarine, aircraft and satellites to evaluate climate model simulations of Arctic sea ice thickness. *The Cryosphere*, 8(5), 1839–1854. https://doi.org/10.5194/tc-8-1839-2014

Thorndike, A. S., Rothrock, D. A., Maykut, G. A., & Colony, R. (1975). The thickness distribution of sea ice. *Journal of Geophysical Research*, 80(33), 4501–4513. https://doi.org/10.1029/JC080i033p04501

Tremblay, B., & Mysak, L. A. (1997). Modeling sea ice as a granular material, including the dilatancy effect. *Journal of Physical Oceanography*, 27, 2342–2360. https://doi.org/10.1175/1520-0485(1997)027<2342:msiaag>2.0.co;2

Tsamados, M., Feltham, D. L., & Wilchinsky, A. V. (2013). Impact of a new anisotropic rheology on simulations of Arctic sea ice. *Journal of Geophysical Research: Oceans*, 118(1), 91–107. https://doi.org/10.1029/2012JC007990

Ungermann, M., Tremblay, L. B., Martin, T., & Losch, M. (2017). Impact of the ice strength formulation on the performance of a sea ice thickness distribution model in the Arctic. *Journal of Geophysical Research: Oceans*, 122(3), 2090–2107. https://doi.org/10.1002/2016JC012128

Wang, Q., Danilov, S., Jung, T., Kaleschke, L., & Wernecke, A. (2016). Sea ice leads in the Arctic Ocean: Model assessment, interannual variability and trends. *Geophysical Research Letters*, 43(13), 7019–7027. https://doi.org/10.1002/2016GL068696

Wang, Q., Danilov, S., Mu, L., Sidorenko, D., & Wekerle, C. (2021). Lasting impact of winds on Arctic sea ice through the ocean's memory. *The Cryosphere*, 15(10), 4703–4725. https://doi.org/10.5194/tc-15-4703-2021

Weiss, J. (2008). Intermittency of principal stress directions within Arctic sea ice. *Physical Review E*, 77, 056106. https://doi.org/10.1103/PhysRevE.77.056106
Weiss, J. (2017). Exploring the “solid turbulence” of sea ice dynamics down to unprecedented small scales. *Journal of Geophysical Research: Oceans*, 122(8), 6071–6075. https://doi.org/10.1002/2017JC013236

Weiss, J., & Dansereau, V. (2017). Linking scales in sea ice mechanics. *Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 375(2086). https://doi.org/10.1098/rsta.2015.0352

Wilchinsky, A. V., & Feltham, D. L. (2006). Modeling the rheology of sea ice as a collection of diamond-shaped floes. *Journal of Non-Newtonian Fluid Mechanics*, 138(1), 22–32. https://doi.org/10.1016/j.jnnfm.2006.05.001

Williams, J., & Tremblay, L. B. (2018). The dependence of energy dissipation on spatial resolution in a viscous-plastic sea-ice model. *Ocean Modeling*, 130, 40–47. https://doi.org/10.1016/j.ocemod.2018.08.001

Willmes, S., & Heinemann, G. (2016). Sea-ice wintertime lead frequencies and regional characteristics in the Arctic, 2003–2015. *Remote Sensing, 8*(1), 4. https://doi.org/10.3390/rs8010004

Zhang, J., & Rothrock, D. A. (2003). Modeling global sea ice with a thickness and enthalpy distribution model in generalized curvilinear coordinates. *Monthly Weather Review*, 131(5), 845–861. https://doi.org/10.1175/1520-0493(2003)131<0845:mgsiwa>2.0.co;2

Zhang, J., & Rothrock, D. A. (2005). Effect of sea ice rheology in numerical investigations of climate. *Journal of Geophysical Research: Oceans*, 110(C8). https://doi.org/10.1029/2004JC002599