# Sea Ice Rheology Experiment (SIREx): 2. Evaluating Linear Kinematic Features in High-Resolution Sea Ice Simulations

Nils Hutter, Amélie Bouchat, Frédéric Dupont, Dmitry Dukhovskoy, Nikolay Koldunov, Younjoo J. Lee, Jean-François Lemieux, Camille Lique, Martin Losch, Wieslaw Maslowski, et al.

► To cite this version:

Nils Hutter, Amélie Bouchat, Frédéric Dupont, Dmitry Dukhovskoy, Nikolay Koldunov, et al. Sea Ice Rheology Experiment (SIREx): 2. Evaluating Linear Kinematic Features in High-Resolution Sea Ice Simulations. Journal of Geophysical Research. Oceans, 2022, 127, 10.1029/2021JC017666. insu-03683298

HAL Id: insu-03683298

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

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
<figure><img src="image_2.png" /></figure>

JGR Oceans

RESEARCH ARTICLE

10.1029/2021JC017666

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

**Key Points:**

* All models simulate linear kinematic features (LKF), but none accurately reproduces all LKF statistics

* Resolved LKF are affected strongest by spatial and temporal resolution of model grid and atmospheric forcing and rheology

* Accurate scaling of deformation rates is a proxy only for realistic LKF numbers but not for any other LKF static

**Correspondence to:**

N. Hutter
nils.hutter@awi.de

**Citation:**

Hutter, N., Bouchat, A., Dupont, F., Dukhovskoy, D., Koldunov, N., Lee, Y. J., et al. (2022). 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

Received 15 JUN 2021

Accepted 4 MAR 2022

**Author Contributions:**

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

**Data curation:** Amélie Bouchat

**Investigation:** Nils Hutter

**Methodology:** Nils Hutter

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

**Resources:** Nils Hutter, Amélie Bouchat, Frédéric Dupont, Dmitry Dukhovskoy, Nikolay Koldunov, 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

**Writing – original draft:** Nils Hutter

© 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.

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

Sea Ice Rheology Experiment (SIREx): 2. Evaluating Linear Kinematic Features in High-Resolution Sea Ice Simulations

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

¹Alfred-Wegener-Institute, Helmholtz Zentrum für Polar- und Meeresforschung, Bremerhaven, Germany, ²Department of Atmospheric and Oceanic Sciences, McGill University, Montréal, Montréal, QC, Canada, ³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, ⁷University of Brest, CNRS, IRD, Ifremer, Laboratoire d'Océanographie Physique et Spatiale (LOPS), IUEM, Brest, France, ⁸Department of Earth and Atmospheric Sciences, University of Alberta, Edmonton, AB, Canada, ⁹Nansen Environmental and Remote Sensing Centre, Bjerknes Centre for Climate Research, Bergen, Norway, ¹⁰Institut de Géophysique de l'Environnement, CNRS, Grenoble, France, ¹¹Danish Meteorological Institute, Copenhagen, Denmark

**Abstract** Simulating sea ice drift and deformation in the Arctic Ocean is still a challenge because of the multiscale interaction of sea ice floes that compose the Arctic Sea ice cover. The Sea Ice Rheology Experiment (SIREx) is a model intercomparison project of the Forum of Arctic Modeling and Observational Synthesis (FAMOS). In SIREx, skill metrics are designed to evaluate different recently suggested approaches for modeling linear kinematic features (LKF) to provide guidance for modeling small-scale deformation. These LKF are narrow bands of localized deformation that can be observed in satellite images and also form in high resolution sea ice simulations. In this contribution, spatial and temporal properties of LKF are assessed in 36 simulations of state-of-the-art sea ice models and compared to deformation features derived from the RADARSAT Geophysical Processor System. All simulations produce LKF, but only very few models realistically simulate at least some statistics of LKF properties such as densities, lengths, or growth rates. All SIREx models overestimate the angle of fracture between conjugate pairs of LKF and LKF lifetimes pointing to inaccurate model physics. The temporal and spatial resolution of a simulation and the spatial resolution of atmospheric boundary condition affect simulated LKF as much as the model's sea ice rheology and numerics. Only in very high resolution simulations (≤2 km) the concentration and thickness anomalies along LKF are large enough to affect air-ice-ocean interaction processes.

**Plain Language Summary** Winds and ocean currents continuously move and deform the sea ice cover of the Arctic Ocean. The deformation eventually breaks an initially closed ice cover into many individual floes, piles up floes, and creates open water. The distribution of ice floes and open water between them is important for climate research, because ice reflects more light and energy back to the atmosphere than open water, so that less ice and more open water leads to warmer oceans. Current climate models cannot simulate sea ice as individual floes. Instead, a variety of methods is used to represent the movement and deformation of the sea ice cover. The Sea Ice Rheology Experiment (SIREx) compares these different methods and assesses the deformation of sea ice in 36 numerical simulations. We identify and track deformation features in the ice cover, which are distinct narrow areas where the ice is breaking or piling up. Comparing specific spatial and temporal properties of these features, for example, the different amounts of fractured ice in specific regions, or the duration of individual deformation events, to satellite observations provides information about the realism of the simulations. From this comparison, we can learn how to improve sea ice models for more realistic simulations of sea ice deformation.

1. Introduction

Continuous and omnipresent deformation turns the Arctic Sea ice cover into a mosaic of ice floes. The deforma-
tions concentrate in narrow bands along floe boundaries, where the ice breaks and ridges in divergent, convergent,
**Writing – review & editing:** Nils Hutter, Amélie Bouchat, Frédéric Dupont, Dmitry Dukhovskoy, Nikolay Koldunov, 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

and especially shear motions. Recently, the focus of the sea ice modeling community on thermodynamics and large-scale circulation of sea ice has extended to resolving these small-scale deformation processes in sea ice dynamics. High resolution applications as well as changes to the model physics describing the material properties of ice allow to explicitly resolve deformation that is localized in narrow lines consisting of segments of leads and pressure ridges. These elongated deformation bands are referred to as Linear Kinematic Features (LKFs).

Leads and pressure ridges represent only a small fraction of the large scale ice cover, but their presence changes the interaction of sea ice with the ocean and atmosphere in the Arctic climate system substantially. The opening of the ice cover in a lead results in intensified heat and humidity exchange between the ocean and atmosphere, resulting in ice growth, brine rejection in the ocean, and convective processes in both the ocean and atmosphere (e.g., Lüpkes et al., 2008; McPhee et al., 2005). The sea ice piled along pressure ridges determine the regional surface roughness which in turn affects the atmospheric and oceanic boundary layer circulation, snow distribution, and drag forces acting on the ice (e.g., Martin et al., 2016). Currently, coarse resolution Global Climate Models resolve very poorly resolve discontinuities in the pack ice (or LKFs). Instead, the effects of leads on heat and freshwater fluxes and ultimately on the Arctic climate are modeled by sub-grid scale parameterizations such as fractional ice cover variables. To directly simulate these processes and to provide a more detailed picture of the complex Arctic climate system, we need sea ice models that explicitly resolve LKFs. A dynamical framework with strongly localized sea ice deformation is the first step toward a realistic representation of leads and pressure ridges in continuum sea ice models that will be used in climate simulations for the foreseeable future (Blockley et al., 2020; Hunke et al., 2020).

Various adjustments have been suggested to improve the representation of LKFs and to resolve leads in contin-
uum sea ice models: (a) increasing the model resolution to a horizontal grid spacing smaller than 5 km (e.g.,
Hutter et al., 2018; Wang et al., 2016), (b) modifying the yield curve (e.g., Bouchat & Tremblay, 2017), and
(c) introducing new rheological frameworks (e.g., Elasto Brittle (EB), Maxwell Elasto Brittle (MEB), Bouillon
& Rampal, 2015; Dansereau et al., 2016; Girard et al., 2011). All three approaches require the convergence of
the dynamics solver to ensure accurate solutions (Lemieux & Tremblay, 2009; Koldunov et al., 2019; Losch
et al., 2014; Wang et al., 2016). In most cases, the observed localization of deformation rates have been assessed
with multifractal scaling analyses to describe the LKF representation in space and time (Bouchat & Trem-
blay, 2017; Hutter et al., 2018; Marsan et al., 2004; Rampal et al., 2016, 2019). The scaling analysis does not
allow unambiguous discrimination between different rheologies in comparison to satellite observations (Bouchat
& Tremblay, 2017; Hutter et al., 2018; Rampal et al., 2016). This raises two questions: First, if all rheologies
perform similarly well, are the scaling analyses a sufficient tool to investigate differences between the rheolo-
gies? Second, scaling analyses give insights into the underlying material properties and deformation physics,
but is this relevant for LKF properties on climate scales (e.g., LKF size, opening times, etc.)? It is plausible that
explicit simulations of the interaction of atmosphere, ice, and ocean associated with sea ice leads may require a
realistic spatial and temporal distribution of LKFs (e.g., Ólason et al., 2021). Hutter and Losch (2020) showed
that scaling analyses alone cannot evaluate this aspect in their simulations with the viscous-plastic (VP) rheology
(Hibler, 1979). However, using a combination of scaling analysis and statistics from automated LKF detection
algorithms (Hutter et al., 2019a; Linow & Dierking, 2017) allows for a comprehensive evaluation of LKFs.

In 2017, the sea ice modeling working group of the Forum of Arctic Modeling and Observational Synthesis (FAMOS, Proshutinsky et al., 2020) launched the Sea Ice Rheology Experiment (SIREx) model intercompari-
son project with two aims: (a) To extend the current research on simulating small-scale sea ice deformation to
additional modeling frameworks, namely all rheologies used in the sea ice modeling on climate scales, and (b)
to develop, compare, and combine new and existing evaluation metrics to gauge the realism of the simulated
features. In total, 10 international groups participated with 36 simulations from 11 different models. The contrib-
uted simulations cover all rheologies commonly used in Pan-Arctic Sea ice simulations with continuum models
(viscous plastic or its elastic-viscous-plastic variation—(E)VP, elastic anisotropic plastic—EAP, and, Maxwell
elasto brittle—MEB), as well as a large range of model resolutions (1–15 km), different atmospheric boundary
conditions (reanalysis with different spatial and temporal resolution as well as interactively coupled atmospheric
models), coupling to different complex ocean fields, and different parameterizations with different effects on the
ice strength (different number of ice thickness categories or ice thickness distribution (ITD) classes and modi-
fied yield curve parameters). The analysis of this suite of simulations is structured in two parts: Part I (Bouchat
et al., 2022b) studies the heterogeneity and intermittency in the simulated deformation fields with established
scaling metrics. Combined with new uncertainty estimates (Bouchat & Tremblay, 2020), this analysis offers
insights into the physical properties of the simulated ice deformation that forms LKFs. They find that power-law
scaling and multifractality of deformations in space and time can be simulated with both plastic and brittle rheol-
ogies. In SIReX Part II—The subject of this paper—We make use of automated detection and tracking algorithms
(Hutter et al., 2019a) to study the spatial and temporal distribution and characteristics of LKFs, for example,
densities, lengths, and lifetimes. This analysis provides a comprehensive description of simulated deformation
features and allows for the evaluation of LKF properties that are highly relevant for interaction processes at the
air-ice-ocean interface.

The objective of this paper is to evaluate the spatial and temporal properties of simulated deformation features in all SIReX models with satellite observations. We use detection and tracking algorithms to extract LKFs in simulated deformation fields and compare them to the RGPS LKF data set (Hutter et al., 2019b). The comparison is made for the two winters (January, February, and March; JFM) of 1997 and 2008. The results of this comparison are interpreted in light of the different model parameters and parameterizations (e.g., rheology or spatial resolution) to assess the impact and importance of individual parameters on the quality of resolved LKFs. We link our feature-based evaluation to the scaling analysis of SIReX Part I (Bouchat et al., 2022b) and to an analysis of sea ice thickness and concentration anomalies along resolved LKFs. This forms an in-depth comparison of different dynamical modeling frameworks for sea ice and their capabilities for simulating localized deformation along floe boundaries. A special focus of this intercomparison project, besides the insights for sea ice rheology and model development, is to provide guidance for users of sea ice models in the context of coupled climate simulations in the Arctic.

## 2. Data

In this section, we introduce the different model simulations and the observational data set that we use for the
comparison.

### 2.1. RGPS LKF Data

The RGPS data set is a high spatial resolution data set of sea ice deformation that is often used for model evalua-
tion (e.g., Bouchat & Tremblay, 2017; Hutter & Losch, 2020; Rampal et al., 2019; Sreen et al., 2017). From the
RGPS deformation data set, an LKF data set was generated (Section 3.1, Hutter et al., 2019a, 2019b (data set))
that contains 165 000 detected and 36 000 tracked LKFs in the winters 1996/1997 to 2007/2008. For this study,
we subsample the RGPS LKF data set to the two SIReX winters (January, February, and March in 1997–2008)
for the model evaluation. The remaining winters (JFM 1998–2007) are used to estimate the interannual variability
of LKF statistics in the RGPS data set.

### 2.2. Model Data

11 models contributed to the second part of SIReX with a total of 22 simulations and 36 simulated winters. The participating models cover a broad variety in terms of rheology, spatial grid resolution, temporal and spatial resolution of atmospheric boundary condition, ITD classes, and grid type (structured vs. unstructured, Lagrangian vs. Eulerian). Some important details of the model simulations are summarized in Table 1; for further information we refer to the specific references provided in Table 1. Note that the MERCATOR model participating in SIReX Part I does not participate in our analysis. In addition to the simulations presented in SIReX Part I, we analyze a MITgcm (2 km, ITD) simulation for the year 1997 and FESOM2 simulations (1997 and 2008). For all simulations sea ice drift, concentration, and thickness were provided for at least one SIReX winter (14 of the simulations cover both years and 8 simulations cover either 1997 or 2008). The SIReX winters have been chosen based on the availability of existing model output and to be representative of a large period of sea ice retreat. Most of the model output is provided as daily mean fields, and only some groups provided daily snapshots. The Probability Distribution Functions (PDFs) of deformation rates are robust to the choice of output diagnostic (snapshots vs. daily mean, results not shown). For this reason, we process both types of output in the same way as in the processing steps outlined in Section 3.2.
<caption>**Table 1**</caption>

Key Parameters of the Model Runs Participating in Sea Ice Rheology Experiment Part II

<table><thead><tr><th>Model/Configuration (Group)<br>Label</th><th>Year</th><th>Grid spacing,<br>Time step</th><th>Grid</th><th>Rheology<br>(solver/nb, it)</th><th>Ice Strength<br>Parameters</th><th>ITD #</th><th>Atm. Reanalysis (A<sub>t</sub>, Δt) or<br>Atm. Model</th><th>References</th></tr></thead><tbody><tr><td colspan="9"><strong>MITgcm (Alfred-Wegener-Institute - AWI)</strong></td></tr><tr><td>MITgcm (2 km, ITD)</td><td>1997, 2008</td><td>2 km, 120s</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 and Losch, (2020)</td></tr><tr><td>MITgcm (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>2</td><td>"</td><td>"</td></tr><tr><td>MITgcm (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>2</td><td>"</td><td>"</td></tr><tr><td>MITgcm (2 km)</td><td>1997, 2008</td><td>"</td><td>"</td><td>"</td><td>P*, e = 22.64, 2.0</td><td>"</td><td>"</td><td>Hutter and Losch, (2020)</td></tr><tr><td>MITgcm (4.5 km)</td><td>2008</td><td>4.5 km, 240s</td><td>"</td><td>"</td><td>"</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, 3600s</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 and 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, 720s</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>Muilwijk 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, 360s</td><td>E</td><td>EVP (120)</td><td>C<sub>j</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, 180s</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>NEMO-LIM2/ANHA12 (University of Alberta)</strong></td></tr><tr><td>ANHA (4 km)</td><td>2008</td><td>4.1 km, 180s</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, 1080s</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, 180s</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>FESOM (Alfred-Wegener-Institute - AWI)</strong></td></tr><tr><td>FESOM</td><td>1997, 2008</td><td>5.1 km, 600s</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>FESOM M2 (Alfred-Wegener-Institute - AWI)</strong></td></tr><tr><td>FESOM M2</td><td>1997, 2008</td><td>1 km, 180s</td><td>U</td><td>mEVP (400)</td><td>P*, e = 27.5, 2.0</td><td>2</td><td>JRA55 (~60 km, 3 hr)</td><td>Wang et al. (2020)</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, 1200s</td><td>E</td><td>EVP (600)</td><td>C<sub>j</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, 1200s</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, 1200s</td><td>E</td><td>EAP (120)</td><td>C<sub>j</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>neXSIM - V1(2018) (Nansen Environmental and Remote Sensing Center - NERSC)</strong></td></tr><tr><td>neXSIM</td><td>1997, 2008</td><td>10 km, 200s</td><td>L</td><td>MEB</td><td>T*, P*, = 21, 75</td><td>3</td><td>CFSR/CFSv2 (0.5° version, 6 hr)</td><td>Rampal et al. (2019)</td></tr></tbody></table>

Note. We use the following abbreviations (see also Table C1 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<sub>j</sub>: Frictional energy dissipation parameter; ITD # is the number of ice-thickness categories in the ice thickness distribution; And (solver/nb, 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. The colors of the run labels refer to the marker and line colors used in all figures to differentiate the participating runs. 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.
## 3. Methods

In this section, we describe how the LKF detection and tracking algorithms are applied and what sampling we use for the model data.

### 3.1. LKF Detection and Tracking Algorithms

At the core of our feature-based evaluation are data sets of LKFs. These are derived from sea ice deformation fields from satellite observations and simulation data by automatic feature detection and tracking. The automatic feature detection and tracking algorithms are described in Hutter et al. (2019a). A brief summary follows.

From the original map of deformation, a binary map of LKF pixels is created that have significantly higher deformation rates than their immediate neighborhood. Then, this binary map is divided into short segments of neighboring LKF pixels. Finally, segments are reconnected based on a probability function that describes their distance, the orientation relative to each other, and the difference in the deformation rate.

The tracking algorithm uses drift information between pairs of subsequent LKF fields to advect the LKFs of the first field. The advected features are then compared to the LKFs in the second field. Tracked LKFs are identified based on the degree of overlap between advected LKFs of the first field and detected LKFs of the second field.

### 3.2. Detection and Tracking of LKFs in the Model Data

The LKF detection and tracking algorithms require the data to be on a regular grid (Hutter et al., 2019a). The output fields from the unstructured grid models were interpolated onto regular grids with spatial resolution similar to their original resolution (FESOM ice velocities to the grid of MITgcm 4.5 km, FESOM2 ice velocities and deformation rates to the grid of MITgcm 2 km, and neXtSIM ice velocities to a grid with 10 km grid spacing in Polar Stereographic projection). All other model output was processed on the native model grids.

We adapt the filtering technique of the LKF detection algorithm to account for the fact that the model and RGPS data do not have the same sampling frequencies (3 days for the RGPS and 1 day for the models). We mimic the 3-day sampling of RGPS by combining three daily deformation fields in the following way: (a) Total deformation rates are computed by finite differences from the daily drift output. (b) Pixels are flagged as LKF pixels in each daily total deformation field with a difference of the Gaussian (DoG) filter in the following way: We scale the original kernel sizes of the DoG filter (radii $r_1 = 1$ pixel and $r_2 = 5$ pixel) by a factor $f = \frac{1}{2} + \frac{\Delta x_{rgps}}{2 \Delta x_{model}}$ to take into account the fact that very high resolution simulations can have finer scale features than RGPS data. Pixels are then flagged as LKF when their total deformation exceeds the average deformation rate in the immediate neighborhood (averaged over a radius of 62.5 km) by $d_{LKF} = 0.01 \text{ days}^{-1}$. This threshold is determined to be fine enough to filter all LKFs, but still high enough to prevent spurious detection of noise in the deformation fields that is caused by a lower accuracy of the solution of the momentum equation. Only the neXtSIM simulation contains LKFs with deformation rates lower than this threshold such that we use a threshold of $d_{LKF} = 0.002 \text{ days}^{-1}$ for this simulation. (c) The three daily deformation fields are combined into one binary map where a specific pixel is flagged as LKF if any of the three daily fields are flagged as LKF at this pixel position. (d) A morphological thinning algorithm is applied to the combined binary map to reduce all LKFs to a width of one pixel. By applying the morphological thinning algorithm explicitly after combining the daily LKF maps, we ensure that LKFs that move within 3 days are not detected more than once. A detailed discussion of the comparability of the RGPS LKF data set and the derived model LKF data sets with respect to spatial and temporal resolution is included in Appendix A.

The detection routines (segment detection and reconnection of Hutter et al., 2019a) are applied to the combined and thinned binary maps. Both algorithms utilize optimized parameters for the RGPS data set (Table 1 in Hutter et al., 2019a), with the minimal length of LKFs scaled by the corresponding model resolution. The simulated ice drift fields are also used, besides for deriving deformation, to advect the LKFs over 3 day intervals (between the
3 daily records) for the tracking algorithm. In the following, we use the optimized tracking algorithm parameters for the RGPS data set (Table 2 in Hutter et al., 2019a).

### 3.3. Principles of the Model-Observation Comparison

The RGPS data set covers most of the Arctic Ocean, but the coverage varies in time depending on the available synthetic aperture radar (SAR) imagery. Gaps in the deformation data can split long LKFs into multiple smaller LKFs or inhibit tracking of LKFs. This affects derived LKF statistics leading to, for instance, fewer LKFs with long lifetimes (Hutter & Losch, 2020). In contrast, the model output fields provide deformation data for the whole study region without any gaps in time. Thus, we mask the LKFs detected in each simulation with the RGPS coverage in order to exclude any effects of varying RGPS coverage on the comparison between the model and observations. The tracking of LKFs is repeated only for the masked LKFs. We use the masked LKFs for all statistics presented in this paper because these results are directly or indirectly affected by the number, length, or lifetime of the LKFs. The original unmasked LKFs are only used in the concentration and thickness anomaly analysis, as we do not compare to RGPS data in this analysis, and in the intersection angle analysis, because in the coarse resolution simulations there are otherwise too few LKFs available for this analysis. We tested with all high-resolution simulations that masking the LKFs with the RGPS coverage does not affect the distribution of simulated LKF intersection angles (not shown).

A direct comparison of LKFs detected in RGPS and model output is not possible due to the chaotic nature of ice fractures (Coon et al., 2007; Kwok et al., 2008). Instead, statistics of the spatial and temporal properties of the simulated features can be compared to the statistics derived from observations (Hutter & Losch, 2020). This paper focuses on a comparison between the probability distribution functions (PDFs) of LKF properties related to space (length, intersection angle) and time (lifetime, growth rates). In addition we study the density of LKFs described by the regional distribution of relative LKF frequency. To properly evaluate the simulations with RGPS, we seek a quantitative way to compare the PDF of a given LKF property obtained from the model and RGPS LKF data. Most of these distributions have a heavy tail so that standard metrics, such as the Kolmogorov-Smirnov statistic or the Wasserstein distance, cannot be applied. Instead we define the skill metric as the integral of the difference between two (more specifically, between RGPS and simulation) PDFs of a given property on a logarithmic scale. We choose this metric as it is closest to the visual comparison of the PDFs and it emphasizes the tail of the distributions. For clarity and consistency, we also use the same statistic to define the misfit in intersection angles, although the data is not heavy tailed. By using the log-scale in the skill metric, all parts of the intersection angle PDFs are weighted without overemphasizing the peak of the distribution. As not all simulations contributing to SIREx are run over the same time period, we use the RGPS years consistent with each simulation in the computation of the metric functional. In this manner, we take interannual variability of the LKF statistics between both SIREx winters into account. By definition, the skill metric increases with larger misfit between models and observations. Thus, a lower skill metric value indicates better model performance.

For reference of the magnitude of the skill metric, we use the interannual variability of different LKF properties within the RGPS data set. To this end, we compute the PDFs of each property for all RGPS winters (JFM 1998–2007) individually. Comparing these PDFs with both SIREx winters, reference values are computed using the same “skill” metric as for comparison between the simulations and RGPS data. For each property, we average the computed “skill” values of all years to quantify the interannual variability that can be compared to the computed models’ skill. We note here that for quantities with little year-to-year variability, for example, the intersection angle, the interannual variability is a reasonable benchmark to assess model performance. For quantities that are thought to be affected more strongly by wind patterns and ice conditions and thus vary stronger from year to year (e.g., LKF density), the interannual variability provides a maximum of the skill metric that should not be exceeded by a model simulation to be called useful.

### 4. Simulated Deformation Features

In this section, we present the characteristics of the simulated LKFs and compare them to satellite observations. First, we describe the overall number of LKFs in each simulation and their regional distribution, followed by an analysis of the spatial properties, length and intersection angle, and the characteristics
related to time, lifetime and growth rate. We include a short metric-specific discussion of the results in this section. The relation between specific model parameters on the skill values of simulated LKFs is discussed specifically in Section 5.

## 4.1. Number of Simulated LKFs

Recent sea ice modeling studies described in a qualitative way how the number of resolved deformation features varies with model parameters (e.g., Bouillon & Rampal, 2015; Sreen et al., 2017; Wang et al., 2016), but only very few quantified these variations (Hutter & Losch, 2020; Koldunov et al., 2019). Visual analysis of deformation fields does not provide enough information to distinguish the tendencies in the number of LKFs and the total length of all LKFs, as both are proportional to the total number of grid cells associated with LKFs. Our object-based approach allows assessing both properties in a quantitative way. In Figure 1a, we present how the numbers of detected LKFs change with model resolution for the SIREx simulations. Figure 1b shows the total length of all LKFs as function of the model resolution. For the discussion, the total length of all LKFs is more interesting, because it is directly related to the total area of potential air-ice-ocean interactions, whereas the number of LKFs does not include the length, and hence, area information.

For our set of simulations we find that both the number of LKFs (Figure 1a) and the total length of all LKFs (Figure 1b) increase as the grid resolution becomes finer, even though for both metrics there is considerable variation also between simulations with the same grid resolution. Note that we cannot comment on the effect of the resolution for EAP and MEB, as there are only simulations using these rheologies in our comparison at one specific resolution. Most models underestimate both the number of LKFs and their total length. MITgcm (2 km, ITD) agrees with RGPS data in the total length of all LKFs but has too many LKFs indicating that too short LKFs are simulated (see the PDFs of LKF length in Figure 4).

Besides the effect of grid resolution, some other model parameters seem to affect the number of LKFs as well: (a) More ITD classes lead to more resolved LKFs in the simulations, and to an increase of the total length of all LKFs as seen by comparing MITgcm (2 km) and FESOM2 to MITgcm (2 km, ITD), or by comparing ANHA (4 km), FESOM, and MITgcm (4.5 km) to HYCOM-CICE (FSU) and RIOPS. All simulations listed here use an (E)VP rheology. (b) The comparison of RASM-WRF (both EAP and EVP simulations) with RASM-CORE2 (EAP) suggests that coupling the sea ice model to an interactive atmosphere rather than forcing it with atmospheric reanalysis increases the number of resolved LKFs. While this is plausible, because prescribed forcing generally suppresses internal variability, the available simulations alone do not allow to conclude if the increase in LKFs is caused by the coupling itself or the higher spatial and temporal resolution of the wind fields driving the ice. We note here that the coupled RASM simulation also use an increased number of EVP subcycles, which also likely improves the LKF representation. A coupling time step of 20 min as in RASM-WRF allows to resolve inertial oscillations leading to higher variable wind forcing that potentially initiates the formation of additional LKFs. Simulations using very high spatial and temporal resolution atmospheric forcing also tend to show better agreement with observations than simulations forced by medium resolution winds (ANHA (4 km), FESOM, and MITgcm (4.5 km) versus HYCOM-CICE (FSU) and RIOPS; Danish Meterological Institute (DMI) versus McGill). (c) The neXtSIM simulation shows results much closer to observations than other simulations with similar resolution. Its grid resolution of 10 km is in the range of all models, still this simulation stands out for the following reasons: It is the only simulation to use a brittle rheology (MEB) and it uses a Lagrangian modeling approach (i.e., an unstructured moving grid) and high resolution atmospheric forcing. This makes it difficult to clearly separate the effects of the individual parameters. The similarity of the two FESOM simulations, which also use unstructured, but fixed grids, with the MITgcm simulations of similar resolution suggests that the type of grid itself has little effect. The relatively high resolution forcing (40 km or less) will add more variability to the system, but does not guarantee good agreement with RGPS observations (see e.g., RIOPS HYCOM-CICE (FSU), DMI), so that the MEB rheology—Potentially when associated with a Lagrangian grid—Appears to be responsible for the better agreement with observations. Clearly, a direct comparison of MEB and VP rheologies in the same model framework (grid, forcing, etc.) is required to illustrate the differences between these different rheology approaches.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** (a) The average number of detected linear kinematic features (LKF's) in 3-day intervals and (b) the total length of all LKF's together in 3-day intervals in all simulations as a function of spatial grid resolution. The interannual variability of the LKF number and length for each model that contributed simulations for 2 years is shown as a thin vertical line around the multiyear mean. The thin gray lines represent the RADARSAT Geophysical Processor System value across all spatial resolution for reference. The dashed gray lines show a linear trend fitted through all model points with statistically significant nonzero slope for both metrics.</figcaption></figure>

## 4.2. Regional Distribution

There are LKF's in the entire Arctic Ocean, but particularly high densities are found along coastlines and bathymetric features (shoals and continental slopes), with slightly pronounced densities in the Beaufort Sea, and low densities in fast ice regions (Hutter et al., 2019a; Mahoney et al., 2012; Wernecke & Kaleschke, 2015; Willmes & Heinemann, 2016). The distribution of LKF's in the SIREX simulations have stronger regional variations than RGPS (Figure 2). In general, there is a large spread in modeled LKF densities and their difference to the observed LKF densities (Figure 3). Indeed, nearly all models strongly underestimate LKF densities by far more than the range of interannual variability of RGPS. As expected from the analysis of the number of LKF's, low resolution models using (E)VP or EAP rheologies tend to underestimate the LKF density more than high resolution models. The largest differences are mainly in the pack ice area (here defined as sea ice regions that are more than 150 km away from the coastline) with generally too few LKF's compared to observations. In coastal regions the distribution of LKF's is better reproduced with half of the models showing differences within the range of the RGPS interannual variability. Along the closed boundary of the coastline stress concentrates and initiates ice fracture. Our results indicate that this process forms LKF's independent of the model's grid resolution.

In the pack ice area, only the MITgcm (2 km, ITD) and the neXtSIM simulation produce overall LKF densities within the interannual variability of RGPS (both in pack ice and the entire Arctic). Both models use mechanisms to locally reduce the ice strength based on the deformation history, which initiates fractures in pack ice
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Linear kinematic feature (LKF) densities for RADARSAT Geophysical Processor System (RGPS) and all models defined as the relative frequency of LKFs within 150 × 150 km boxes in a polar stereographic projection normalized by the number of deformation observations available for the box. Cells without any LKF are indicated in faded colors. Note that the RGPS coverage for both Sea Ice Rheology Experiment years is different leading (see first row) to a slightly different masks for the simulations.</figcaption></figure>
<figure><img src="image_6.png" /></figure>

**Figure 3.** Cumulative box-wise differences of linear kinematic feature (LKF) density for all simulations compared to RADARSAT Geophysical Processor System (RGPS) separated into three regions: Entire Arctic (left), within 150 km of the coastline (middle), and pack ice area 150 km away from the coastline (right). Regions where models underestimate LKF densities are accumulated in the blue bar, while the overestimation of LKF densities is given by the red bar. The two color maps indicate the size of the over- or underestimation in each grid cell. In the bottom row, we illustrate how the RGPS LKF densities of the other RGPS years differ from the two Sea Ice Rheology Experiment winters. Two boxplots show this interannual variability of both under- and overestimation. In the boxplot, the orange center lines denote the mean, the box extends from the lower to upper quartile, and the whiskers show the 1.5 interquartile range. The dotted gray lines show the first quartile of the underestimation and the third quartile of the overestimation of LKF densities in all RGPS years as a reference range for LKF densities.

(the direct deformation feedback by damage parameterization in neXtSIM and the indirect deformation feedback on the ice strength defined in Rothrock (1975) by preferred opening in lower ITD classes). We analyzed only two model winters, but the interannual variability of LKF distribution in the RGPS data set is small, so that the underestimation of sea ice deformation in the pack ice area by nearly all models is very likely a general issue in every year.

### 4.3. LKF Length

Deformation features crisscross the sea ice cover at ranges from meters to multiple kilometers (Kwok, 2001). The distribution of LKF lengths in the RGPS data set is heavy-tailed and can be described by a stretched exponential distribution (Hutter et al., 2019a), which is consistent with the notion of scale-invariant sea ice deformation. The LKF length distributions of all SIREx models also have fat tails (Figure 4), even though most models overestimate the number of very long LKFs. This phenomenon is most pronounced for models with fewer LKFs (see Section 4.1). Some models (FESOM2, HYCOM-CICE (FSU), MITgcm (2 km, ITD), MITgcm (2 km, e = 1, ↓P), and MITgcm (2 km, e = 0.7, ↓P)) overestimate the amount of small LKFs thereby overly steepening the tail of the distribution. An obvious reason for this is that if more LKFs cover the same area, it is more likely that the intersection with other LKFs makes them shorter (Hutter & Losch, 2020). Only HYCOM-CICE (FSU) and both MITgcm (2 km, ↓e, ↓P) simulations with reduced ellipse ratios reproduce the LKF length distribution within the range of the reference, which represents the interannual variability of the RGPS data set.
<figure><img src="image_7.png" /></figure>

**Figure 4.** (a) Distribution (probability distribution functions (PDFs)) of linear kinematic feature (LKF) lengths for all models and RADARSAT Geophysical Processor System (RGPS). The color references for the models are provided in (b). The difference between the PDFs of simulated LKF lengths and RGPS LKF lengths is shown in (b) for each model. A scalar metric for LKF length is the integrated area between the LKF length PDF of the model and RGPS, given on the right hand side. As reference, we compute the differences between the LKF length PDFs of all RGPS years and the two Sea Ice Rheology Experiment years to quantify the interannual variability. The differences in all bins and all RGPS years are summarized in the box plot given in the color bar. The reference value of this metric for the interannual variability is computed and provided at the right hand side of the color bar.

4.4. Intersection Angle

The material properties of sea ice affect the intersection angle of its deformation features (Erlingsson, 1988; Wang, 2007), and intersection angles have been used to evaluate the dynamics of high resolution sea ice simulations (Heorton et al., 2018; Hutter & Losch, 2020; Ringeisen et al., 2019). Acute intersection angles ranging around 30°–50° have been reported from both satellite imagery (Cunningham et al., 1994; Schulson & Hibler, 2004; Walter & Overland, 1993; Wang, 2007) and laboratory measurements (Schulson et al., 2006). These studies focused on intersecting LKFs, so-called conjugate faults, which form simultaneously under the same compressive forcing. The distribution of intersection angles in the RGPS LKF data set also peaks in this range (see Figure 5 and Hutter et al., 2019a). For an intersecting pair of LKFs, two intersection angles can be computed, both of which add up to 180°. For conjugate faults the fracture angle of each LKF $\delta$ is measured relative to the direction of compressive stress with the intersection angle summing up to $2\delta$. Except for laboratory experiments, the direction of stress is unknown, also in the SIREX simulations and RGPS. The stress states causing the deformation, however, can be also deduced from the resulting sea ice drift (see Appendix B). In this study, we introduce a new approach using the vorticity of two intersecting LKFs to interpret the deformation behavior at the intersection.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** (a) Distributions (probability distribution functions (PDFs)) of intersection angles between pairs of linear kinematic features for RADARSAT Geophysical Processor System (RGPS) and all models. The color references for the models are provided in (b). The difference between the PDFs of simulated intersection angles and intersection angles observed with RGPS is shown in (b) for each model. A scalar metric for intersection angles is defined as the integrated area between the intersection angle PDF of the model and RGPS. For each model, the value of this metric is given at the right hand side of (b). The interannual variability in the RGPS data set is presented in the same way as in Figure 4. Note, that the intersection angles are not periodic data as we define them relative to the stress direction.</figcaption></figure>

The vorticities of two LKF s of a conjugate fault are of different signs, such that we can determine a main direction of compressive stress and choose an intersection angle (0°–180°) relative to this direction. Using this method, we determine that roughly one third of intersecting pairs of LKF s in RGPS are conjugate faults. If the vorticities of both LKF s have the same sign, there is no clear stress direction. In this case, both potential intersection angles are taken into account when computing the distribution of intersection angles but weighted down by a factor of two. The distribution of intersection angles is computed for all SIREx models (Figure 5). We do not mask the LKF s with the RGPS coverage, because some models with very few LKF s do not show intersecting LKF s within the region covered by RGPS. For all other models with intersecting LKF s within the RGPS coverage, we have tested that the region of the analysis does not affect the distribution of intersection angles significantly (not shown).

We find that the distributions of intersection angles in all models peak around 90°, in strong contrast with the peak at 45° for RGPS observations. Only for the McGill simulations with nonstandard yield curves (e = 1, ↓S; e = 0.7, ↓P; e = 1, ↓P; these values of e are outside the commonly used range, but they still follow the principle that ice is stronger in compression than in shear), there is a peak around 55°. However, these simulations also produce a second peak of similar magnitude at 180°–55° = 125°. This peak can mean that (a) the corresponding intersecting pairs are not formed under uniaxial compressive forcing, or (b) the yield curve settings allow a
wider range of intersection angles. The normal flow rule used in these models determines the fracture angle to be normal to yield a curve at the stress state of the fracture (Ringeisen et al., 2019). By decreasing the ellipse aspect ratio *e*, the fracture angle varies more strongly with the stress state assuming that shear deformation is dominant. This could cause the larger variety of intersection angles in these simulations. This interpretation is supported by the fact that the MITgcm and IFREMER simulations with reduced ellipse ratio (*e* = 0.7 and *e* = 1) simulation leads to a broader distribution of intersection angles compared to the simulations with the standard ellipse ratio (*e* = 2).

The three simulations with rheologies for which the previous reasoning does not apply (neXtSIM, RASM-CORE2 (EAP), and RASM-WRF (EAP)) also produce a peak at 90°, while compared to (E)VP simulations the peak is broader showing similar probabilities for intersection angles between 70°–110°. The strong underestimation of small intersection angles compared to RGPS, especially in the neXtSIM simulations, results in a larger overall misfit between simulated and observed PDF distributions than for (E)VP simulations.

Experiments with simple geometry and uniform forcing suggest a clear connection between rheology and intersection angles (e.g., Ringeisen et al., 2019, 2020). The SIREx simulations with realistic forcing and realistic geometry produce a large variety of stress states that eventually lead to deformation. This is manifested in the broad distributions of intersection angles (Figure 5). Some models that use atmospheric boundary conditions with low resolution have a better representation of the intersection angles compared to observations. We speculate that the stronger gradients in a high resolution atmospheric boundary condition partly imprint on the simulated deformation field, which reduces the impact of the rheology on simulated intersection angles and explains why all rheologies lead to similar angle distributions. However, observations point toward a distinct peak of intersection angles that depends less on the atmospheric conditions. Thus, the apparent connection between wind forcing and simulated intersection angles in the simulations in this study does not appear realistic and points to incorrect model physics (Ringeisen et al., 2019, 2020). We can rule out the grid geometry as the primary source of the peak at 90° because it also appears in the simulations on unstructured grids (FESOM, FESOM2, and neXtSIM). In the MITgcm (2 km) simulations, a thorough analysis did not reveal a preferred orientation of LKFs along the grid lines, either (Hutter & Losch, 2020).

## 4.5. LKF Lifetimes

The lifetime of an LKF describes the period when the LKF actively deforms. The lifetime of LKFs together with the area covered by LKFs are key factors for air-ice-ocean interaction. For instance, the overall heat loss through a lead strongly depends on its opening time. The majority of LKFs are active for less than 3 days, but the distribution of LKF lifetimes shows an exponential tail such that LKFs up to a month old can be observed as well (Hutter et al., 2019a; Kwok, 2001). For both SIREx years, RGPS recorded on average shorter LKF lifetimes compared to the other RGPS years. This is likely to be caused by exceptional wind forcing leading to enhanced LKF formation and rapidly varying satellite coverage that makes tracking of LKFs more difficult. Here, both effects reduce the LKF lifetimes. As we use atmospheric reanalysis to drive the simulations and mask the simulated LKFs with the satellite coverage for this analysis, all models are expected to reproduce these shorter lifetimes.

All SIREx models agree with RGPS in the sense that LKFs younger than 3 days are most abundant. In addition, the lifetime distributions of simulated LKFs follow an exponential tail for all SIREx models (Figure 6). The only exception is ANHA (4 km) with no LKFs older than 15 days. In both IFREMER simulations there are no LKFs lifetimes longer than 21 days.

Most models overestimate the relative frequency of long-lived LKFs and produce a too slow decay rate of the exponential tail. This is particularly the case for simulations with low absolute LKF numbers, for which LKFs concentrate in coastal regions. Coastal LKFs, like flaw leads, are a frequently occurring and stable phenomena, as they are mainly determined by coastal geometry. Thus, LKF lifetimes are likely to be overestimated in simulations that show too many coastal LKFs compared to LKFs in the pack ice area. Simulations with a better representation of LKFs in the pack ice area also reduce the overestimation of LKF lifetimes (e.g., FESOM2, RIOPS, HYCOM-CICE (FSU), MITgcm (2 km), and MITgcm (2 km, *e* = 1, ↓P)). We speculate that this remaining overestimation may have its root in the deformation ice-strength feedback initiating LKF growth in nonbrittle models: the strong deformation localized in LKFs leads to a reduction in concentration and thickness, which reduces the ice strength.
<figure><img src="image_9.png" /><figcaption>**Figure 6.** (a) Distribution (probability distribution functions (PDFs)) of Linear kinematic feature lifetimes for RADARSAT Geophysical Processor System (RGPS) and all models. The color references for the models are provided in (b). The differences between the PDFs of simulated lifetimes and lifetimes observed with RGPS are shown in (b) for each models. The scalar metric for lifetime is defined as the integrated area between the lifetime PDF of each model to RGPS normalized by the interannual variability in the RGPS data set. For each model, this value is given at the right hand side of (b). The interannual variability in the RGPS data set is presented in the same way as in Figure 4.</figcaption></figure>

In the next time step, the reduced ice strength makes deformation in the same grid cell more likely. The locally incurred thickness and concentration anomalies can only be reduced by converging ice motion or thermodynamic ice growth. The overestimation of LKF lifetimes by VP, EVP, and EAP models indicates that the memory of past deformation associated with this feedback may be too strong. Instead of this deformation ice-strength feedback, the MEB rheology of the neXtSIM model uses a subgrid-scale damage parameterization that scales the ice strength based on how strongly the stress exceeds the Mohr-Coulomb yield criterion. The memory of the damage parameter is limited by a healing time parameter that is independent of the thermodynamic parameters. Although neXtSIM overestimates the LKF lifetimes as other models in our comparison, neXtSIM simulations with a modified damage criterion (not shown) have shown better agreement with RGPS LKF lifetimes. Thus, tuning the healing parameter could probably improve the overestimation of LKF lifetimes in neXtSIM.

## 4.6. LKF Growth Rates

In this section, we study how quickly LKFs grow and shrink. We define the LKF growth rate as the change in length of an LKF between two records (for definition see Hutter & Losch, 2020). In our analysis, we combine three different growth rates: the initial growth of the LKF at formation and the growing and shrinking of persistent LKFs. We find that the distribution of growth rates shows a heavy tail for RGPS data and all models (Figure 7).
<figure><img src="image_10.png" /><figcaption>**Figure 7.** (a) Distribution (probability distribution functions (PDFs)) of linear kinematic feature growth rates for RADARSAT Geophysical Processor System (RGPS) and all models. The color references for the models are provided in (b). The differences between the PDFs of simulated growth rates and growth rates observed with RGPS are shown in (b) for each models. The scalar metric for growth rates is defined as the integrated area between the lifetime PDF of each model to RGPS. For each model this value is given at the right hand side of (b). The interannual variability in the RGPS data set is presented in the same way as in Figure 4.</figcaption></figure>

The majority of the models overestimate large growth events. We note that most of these models also produce too few LKFs (Section 4.1). Models that overestimate the number of LKFs (neXtSIM and MITgcm (2 km, ITD)) show fewer large growth rate events. We find a similar dependence on the number of LKFs for the distribution of LKF lengths and note here that the models' performance in terms of LKF length and growth rates are correlated. This correlation emerges because the temporal resolution of 3 days is not high enough to accurately observe the speed with which a fracture in the ice can propagate through the entire Arctic Ocean. Propagation speeds in the sea ice cover are linked to elastic wave speeds and range from 10 to 1000 m/s (Dempsey et al., 2012; Marsan et al., 2011; Stamoulis & Dyer, 2000). Thus, the growth rates computed for our 3-daily data are limited by the length of the LKFs. Still, all reported growth rates are within a physically reasonable range. In particular, the fact that we find very high growth rates in all simulations shows that all rheological frameworks in our comparison simulate fast fracture propagation. A temporal sampling rate of seconds to minutes would be required to avoid the limiting effect of LKF length and determine the actual fracture speed in the models.
## 4.7. Summary of LKF Properties

Values of all LKF statistics and all simulations are summarized in Figure 8, where all skill metrics are normalized with the RGPS interannual variability. A cumulative skill metric is computed from the weighted average of all model metrics. In this weighted average, both density metrics (distribution and total length) are considered and the lifetime metric is weighted by a factor two in order to take into account that these skill metrics are most relevant to interaction processes along LKFs in coupled climate simulations. All other skill metrics are weighted with a factor of one in the computation of the cumulative skill metric. We find that no model in the comparison is able to reproduce the spatial and temporal statistics of LKFs within the range of the interannual variability of the RGPS data set, as all simulations have cumulative skill metrics larger than 1. The best simulations in the comparison (MITgcm (2 km, *e* = 1, ↓P), MITgcm (2 km, *e* = 0.7, ↓P), MITgcm (2 km), HYCOM CICE (FSU), and RIOPS) show cumulative skill metric values between 1 and 2. The cumulative skill metric values of MITgcm (2 km, ITD), DMI, FESOM2, and neXtSIM range between 2 and 3 show an overall good performance but too high values for specific LKF statistics. The cumulative skill metric values larger than 3 of MITgcm (4.5 km), McGill, IFREMER, ANHA, FESOM, and RASM are caused by large skill metric values for multiple LKF statistics. A common feature within these simulations is that they strongly underestimate the number of LKFs in the pack ice area, which negatively affects other statistics.

The good performance of some simulations is related to a combination of several aspects: High resolution atmospheric boundary conditions, a high number of ITD classes, high spatial resolution, or a brittle rheology. In contrast, simulations with a large time step have higher cumulative skill metric values. In Section 5 we will discuss these dependencies in more detail.

In the summary table of all model skill metrics (Figure 8) links between different LKF statistics become apparent: (a) An overestimation of LKF density coincides with too short LKFs (MITgcm (2 km, ITD)) and vice versa (McGill, FESOM, and RASM). (b) Simulations that overestimate LKF lengths likely overestimate LKF lifetime as well. (c) The skill metrics for LKF length and LKF growth rates are correlated, because the low temporal resolution of the data does not allow to study the actual propagation speed of the ice fracture.

## 5. Relationship Between Model Configuration Parameters and Deformation Features

In this section, we study how different parameters of the model configuration affect the LKF statistics to provide advice for the configuration of sea ice models. The following configuration parameters are considered: grid spacing, time step, number of ITD classes, spatial and temporal resolution of the atmospheric boundary conditions, and rheology. The dependence of the model performance on the choice of the parameters is investigated by employing a linear regression (Figure 9). The wide spread of skill metric values for similar configuration parameters illustrates the high sensitivity of sea ice models to parameter combinations and different code implementation of model physics. Despite the large spread, some trends emerge and provide guidance on how to appropriately choose model parameters. We need to be cautious with making inferences from the intersection angle metric, as all simulations produce too large intersection angles. Some models with a slightly better skill in intersection angles use somewhat extreme parameters (e.g., 3600 s time step and very coarse resolution atmospheric boundary conditions). The implications may be misleading, also because these models usually have poor skill scores in other LKF statistics such as density and lifetime. We note that, rather than tuning configuration parameters studied in this paper, model physics need to be adapted to improve the intersection angles skill.

From all configuration parameters, the time step of the simulation and the spatial resolution of the atmospheric boundary condition have the strongest effect on the skill metrics, mainly their density, length, and growth rates (*r*² > 0.25). Increasing the spatial resolution of the model only improves the density of deformation features. Thus, the effect of a small time step on LKF length and growth rate is independent of the coupling between high resolution and short time step. Linear kinematic feature lifetimes are not significantly affected by the choice of any of considered configuration parameters. A high temporal resolution of the atmospheric boundary condition improves the simulated LKF length distribution and LKF growth rates (*r*² ≈ 0.2). Horizontal grid spacing (Hutter et al., 2018; Sreen et al., 2017), rheology and ice strength parameters (Bouchat & Tremblay, 2017; Girard et al., 2011; Hutchings et al., 2005; Rampal et al., 2016), and solver accuracy (Koldunov et al., 2019; Wang et al., 2016) have been reported to have an effect on the amount of LKFs in Pan-Arctic simulations, while
Skill metrics

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Summary table of the individual model skill metrics. Each entry provides the absolute value of the metric as a number along with a gray color code. All model values are normalized with the interannual variability for the given property in the RADARSAT Geophysical Processor System (RGPS) data set, that is, model values within the RGPS interannual variability are smaller than one. The magnitude of the underestimation and overestimation contributing to the skill metric value are shown in blue and red color codes. For density, two skill metric values are computed: The sum of over- and underestimation as a measure for the entire distribution of Linear kinematic feature (LKFs) and the difference between both as a measure for the overall area covered by LKFs. The last column provides the weighted sum of all skill metrics for each model. To put more emphasis on properties that are mostly relevant to air-ice-ocean interaction processes along LKFs, we take both density metrics (distribution and total length) into account and scale the lifetime metric by a factor of 2. All other skill metrics are weighted with a factor of 1. Details of the model configuration that help to explain the overall performance of individual models are provided by footnotes to the cumulative skill metric.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Box plots showing the dependence of the individual skill metrics on the following configuration parameters used in the simulations: (a) Spatial resolution of the model, (b) temporal resolution, (c) number of ITD classes, (d) spatial resolution of the atmospheric forcing, (e) temporal resolution of the atmospheric forcing and (f) rheology. The skill metric values on the y-axis are normalized. For each combination of skill metric and model parameter, a linear regression analysis is performed and the r²-value of all regressions is summarized at the top of each subfigure. The rheologies in (f) are ordered to maximize the correlation to the skill metrics. The box plots are faded with decreasing r²-value, where saturated colors represents skill model-parameter combinations with r² > 0.25.</figcaption></figure>

so far time stepping and atmospheric boundary conditions have not been studied in realistic model configurations. In idealized sea ice model configurations, improvements in the scaling characteristics of sea ice deformation were found for experiments driven by high-resolution atmospheric forcing fields (Hutter, 2015) in agreement with the findings presented here.
More simulations using an MEB or EAP type of rheology are needed to robustly disentangle the effect of the rheology on LKF simulations from the effect of other parameters. Both rheologies are only used by one model each. All RASM configurations show very high cumulative skill metric values, also for the EAP simulations. Comparison of the RASM simulations using the EAP rheology versus RASM-WRF (EVP) shows that LKF intersection angles are better reproduced; however, the skill of the simulation with the EAP rheology in terms of the other LKF properties is lower than that in the RASM-WRF. The neXtSIM simulation with the MEB rheology generally reproduces LKFs statistics in closer agreement with RGPS than the (E)VP simulations at the same resolution. The neXtSIM model is also the only model running on a Lagrangian grid, which therefore does not allow us to undoubtedly conclude from this study on the respective importance of the rheology versus the advection scheme. In addition, the neXtSIM simulation are driven with relatively high resolution atmospheric forcing (~40 km), which may have a positive effect on the LKF skill metrics. The combination of its brittle rheology and the use of a Lagrangian moving mesh allows neXtSIM to achieve a good representation of LKFs despite the coarser grid spacing—As it was designed to do.

Three (E)VP models (MITgcm, McGill, and IFREMER) in the comparison run sets of simulations differ only in yield curve parameters, namely the compressive ice strength parameter *P*\* and ellipse ratio *e*. In these simulations, the yield curve parameters have an impact on the LKF quality. For all models we find that simulations with a reduced ellipse ratio show improved LKF densities compared to the simulations with the standard ellipse ratio *e* = 2, as one could expect from the more localized deformation in these simulations (Bouchat & Tremblay, 2017; Bouchat et al., 2022b). For other LKF statistics the effects vary from model to model and might therefore depend on other configuration parameters that are different between the models, for instance resolution. We note, however, that three of the five models with the lowest cumulative skill metric value (<2) in our comparison use a reduced ellipse ratio (MITgcm (2 km, *e* = 1, ↓P), MITgcm (2 km, *e* = 0.7, ↓P), and RIOPS (*e* = 1.5)). In contrast to the ice strength, the ellipse ratio has rarely been considered as a tuning parameter in systematic parameter optimizations of sea ice models. Nevertheless, the few studies that include the ellipse ratio in their optimization to our knowledge (Sumata et al., 2019; Ungermann et al., 2017) show that a reduced ellipse ratio also improves the simulated sea ice extent, volume, and drift. Therefore, we recommend to also adjust the ellipse ratio to reduce the overall model-observation misfit.

We note here that our quantitative analysis of configuration parameters on deformation features does not include the effect of the oceanic drag, even though eddies (e.g., Cassianides et al., 2021), tides (e.g., Heil et al., 2008), and circulation patterns (e.g., Wang et al., 2021; Willmes & Heinemann, 2016) are known to impact sea ice dynamics. However, an analysis of the effect of these oceanic features requires a comprehensive analysis of the ocean fields for all simulations participating in SIREX, which has not been carried out by all groups given the different research interests and availability of ocean model output. We refrain from representing the complexity of the oceanic fields by their spatial and temporal resolution as done for the atmospheric boundary condition, because the horizontal resolutions of the coupled ocean and sea ice models are closely tied to each other, if not the same. For instance, the eddy-rich simulations in the comparison use the same grid in the ocean and sea ice model. Thus, it is not possible to unambiguously separate the effect of increased ocean complexity (e.g., eddies) and increased complexity in sea ice dynamics (e.g., heterogeneity of the ice strength) on the simulated sea ice deformation. The effect of the spatial resolution on the LKF statistics presented in Figure 9a therefore also includes any potential effect of more detailed ocean fields. We note however that the dynamical impacts of eddies on the sea ice is generally limited to the marginal ice zone, where the rheology is less important in the sea ice dynamics due to low ice concentration and a highly fractured ice pack (e.g., Manucharyan & Thompson, 2017). We argue that other aspects and parameterizations of the coupled ocean models will have only small effects on the deformation fields as most of the ocean parameterizations will affect thermodynamics more than dynamics. An exception may be the ice-ocean drag parameterization, 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. On short spatial and temporal scales, the stress divergence term may still be more important. We suggest a dedicated study of the interplay between ocean complexity and simulated sea ice deformation within a single model framework that allows to disentangle the effects of increased resolution in both ocean and sea ice models through dedicated sensitivity experiments.
## 6. Relationship Between LKF Statistics and the Scaling Analysis of Deformation

Scaling analysis of sea ice deformation has been the main tool to evaluate the realism of lead-permitting sea ice simulations. There are indications, however, that the computed scaling exponents are linked to the number of LKFs, but do not provide insights into other LKF statistics such as LKF lifetime or intersection angles (Hutter & Losch, 2020). In this section, we study if the deformation statistics obtained in SIREx Part I (Bouchat et al., 2022b) are linked to the LKF statistics in this study. To this end, we define skill metrics for the deformation rate PDFs (as the sum of the integrated area between the divergence and shear PDF of a SIREx simulation and RGPS) and for all scaling parameters (squared relative error of the parameter compared to RGPS, e.g., $((β_{model} - β_{RGPS}) / β_{RGPS})^2$). In this analysis, we group and average the multifractal parameters (degree of multifractality $μ$, fluctuation exponent $H$, and degree of heterogeneity $C_1$—Bouchat et al., 2022b) together to one skill metric each for spatial and temporal multifractal scaling. We test the correlation between all these deformation statistics metrics and all metrics defined for the LKF statistics (Figure 10). From all possible combinations, we only find a significant correlation (>50%) of the metrics for the PDF of deformation rates and LKF density and total length, as well as of the metrics for the spatial scaling exponent $β$ and spatial multifractal scaling with the LKF density and total length (Figure 10). All other combinations are not clearly correlated. This shows that good agreement of a simulation with observations in terms of the temporal (multi-fractal) scaling analysis or the more complicated spatio-temporal coupling of multifractal statistics does not guarantee realistic LKFs in the simulation. As the LKF density and total length are linked to the number of LKFs, our results generalize the findings of Hutter and Losch (2020) to a broader set of simulations. The link between the PDF of deformation rates and LKF density allows to use the PDF of deformation rates as an easy-to-compute metric to quickly assess a sea ice simulation before using more sophisticated analysis such as our feature-based comparison for a thorough evaluation.

## 7. Anomalies in Sea Ice Concentration and Thickness Along LKFs

So far we have defined LKFs as bands of deformation rates that exceed the local neighborhood. In the context of climate simulations, however, it makes sense to link LKFs to leads and pressure ridges, that is, reduced or increased sea ice concentration or thickness, because the interaction of sea ice with the atmosphere and ocean strongly depends on the anomalies in the concentration and thickness fields. For instance, the sea ice concentration and thickness in a lead determines the size of the heat fluxes from the ocean to the atmosphere. We compute estimates of concentration and thickness anomalies along LKFs for all SIREx simulations and visualize them in two dimensional PDFs (Figure 11). We define the anomalies as the difference between the average concentration along an LKF and the local mean concentration, which is computed as the average around the LKF weighted by a Gaussian kernel of 150 km. We use the same definition for thickness anomalies. We determine the kernel size by balancing the two effects: (a) To be large enough to average out the LKF information itself, but (b) still be small enough to take into account regional variations in the sea ice concentration and thickness. The anomaly analysis is restricted to pack ice regions by taking only LKFs into account that are at least 150 km away from the coast.

The majority of LKFs shows little to no variation in the concentration and thickness field (on average + 10.7 cm in pressure ridges and -12.8 cm in leads and +0.5% in closing and -1.0% in opening). Given the km-scale resolution of the models in this comparison, it makes sense that only some deformation features are associated with large-scale opening and closing, while the majority of LKFs represents smaller scale or pure shear deformation. The anomalies are caused by the divergent and convergent component of the deformation field along the line of failure of the ice that forms the LKF. The anomalies feed back into the ice strength favoring further deformation in sea ice models using the standard Hibler (1979) ice strength. This positive feedback cycle of deformation—Reduced concentration—Further deformation can amplify the initial anomalies, make them persistent, and thereby lead to more LKFs. Obviously, the magnitude of the anomalies depends on the model grid resolution, as a much larger area needs to be opened or ridged in a coarse resolution model to obtain the same variation in sea ice concentration and thickness. Therefore the feedback cycle is enhanced with increasing resolution and the number of LKFs' increases: The concentration anomalies in models with a resolution of ~5 km are as large as -10%, while the simulations MITgcm (2 km) and FESOM2 (1 km) with even higher resolution feature LKFs with concentration reduction of the order of 50% and higher. Also the average reduction over all openings in these simulations is higher than in the coarse resolution models.
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Relationship between skill metrics of the scaling analysis of deformation rates (Sea Ice Rheology Experiment Part I, Bouchat et al., 2022b) and the skill metrics of the Linear kinematic feature statistics. All skill metrics are normalized. For each combination of skill metrics a scatter plot is shown, where the correlation coefficient between both skill metrics is given in upper right in each subplot. Models are represented by their color codes given in Table 1 or Figure 1.</figcaption></figure>

In general, we find that thickness anomalies are more pronounced than concentration anomalies especially for coarse resolution models. In these models, the initial anomalies caused by deformation are smaller given the smaller effect of advection due to the large grid spacing. Therefore, the reduction in ice strength is less likely sufficient for further deformation and advection. Thermodynamic growth closes the ice cover, which slows down the deformation ice-strength feedback significantly and inhibits the deformation (the Hibler (1979) ice strength depends exponentially on the concentration and linearly on thickness). The thin ice in openings grows vertically only slowly by thermodynamic ice growth, so a long-term memory of deformation is retained in the ice. Under changed forcing conditions these long-term thickness anomalies are likely to be a seeding point for new deformation. To catch the effects of very short-lived concentration anomalies in climate simulations, short coupling time steps may be useful.

High-resolution simulations have the tendency to produce more negative thickness anomalies associated with leads than positive thickness anomalies that are found along pressure ridges, because the deformation ice-strength feedback accelerates opening. Only the MITgcm (2 km, ITD) that uses an active ITD together with the Rothrock (1975) ice strength stands out in this respect with a considerable amount of positive thickness anomalies with an average ridging of 27.6 cm. Simulations with reduced aspect ratios of the elliptical yield curve also produce higher mean ridging along LKFs (MITgcm, McGill and IFREMER), as with a “fat yield curve” convergent deformation is more likely for compressive stress states with confinement (Bouchat & Tremblay, 2017).
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Concentration and thickness anomalies along linear kinematic features (LKFs) for all models given in two dimensional probability distribution functions (PDFs). The anomaly is defined as the difference of average concentration and thickness along an LKF and the local mean of these fields that is determined by an average weighted by a Gaussian kernel of 150 km. Only LKFs in pack ice area and further than 150 km away from the coast are considered. The numbers summarize the mean concentration anomalies for opening and closing as well as the mean thickness anomalies for pressure ridges and leads given in centimeters. A typical lead shows both reduced ice thickness and concentration, whereas pressure ridges are characterized by increased ice thickness. Note that, the interpolation of the model output of the simulations on unstructured grids (FESOM, FESOM2, and neXtSIM) to regular grids leads to a slight underestimation of the anomalies for these simulations.</figcaption></figure>

8. Conclusions

Linear kinematic features emerge in the sea ice simulations with increasing resolution and/or with improved sea ice physics. In this comparison between different sea ice models and model configurations, only very few models reproduce some statistics of LKF properties, namely density, number, length, and growth rate, within an acceptable range as defined by the interannual variability of satellite-based RGPS deformation data. Most models, however, simulate unrealistic LKF distributions. This is particularly true for the intersection angle between pairs of LFKs, where none of the participation models reproduces the relative frequency distribution of observed angles. The LKF lifetime is also overestimated in nearly all simulations in the comparison. Even among the models with the highest skill scores, none simulates more than three LKF characteristics within the reference skill based on the interannual variability of satellite observations.

The models that have some skill use either high grid resolution, short time steps, different physics (brittle rheology or modified (E)VP yield curve), different numerical techniques (Lagrangian moving grid), high-resolution
atmospheric boundary conditions in space and time, or a combination of these factors. To advance sea ice model dynamics, one should carefully choose these parameters. But it also means that it is in principle possible to reproduce at least some LKF properties with currently available modeling techniques.

Model resolution, rheology, and solver accuracy have been reported before to affect the number of resolved LKFs (Bouchat & Tremblay, 2017; Girard et al., 2011; Hutter et al., 2018; Koldunov et al., 2019; Rampal et al., 2016, 2019; Sreen et al., 2017; Wang et al., 2016). Here, we find that the temporal resolution of the model (i.e., the time step) and the spatial resolution of the atmospheric boundary conditions also have a strong impact on the spatial and temporal characteristics of resolved deformation features. The temporal resolution of the atmospheric boundary conditions also affects the simulated LKF properties but with a lower impact. The dependence on the time step length may be related to the solver convergence, as shorter time steps imply smaller changes in external forcing between time levels and hence less work for the numerical solver.

From our analysis, it is not possible to unambiguously identify a single factor determining the accuracy of LKF representation in a model, because some factors only appear in combination. Still, our analysis suggests that a sea ice model produces more realistic LKFs if the model configuration meets at least three of the following requirements: (a) Time step smaller than 400 s, (b) grid spacing smaller than 5 km, (c) brittle visco-elastic rheology in combination with a Lagrangian grid, and (d) atmospheric boundary conditions at resolutions smaller than 60 km. The number of ITD classes and the temporal resolution of the atmospheric boundary conditions appear less important to explain the skill metrics scores. Thus, we recommend to choose the four parameters above with care when setting up sea ice simulations. Explicitly, models with a viscous-plastic rheology will always require very high grid resolution (see also Appendix A). The skill metrics defined in our study can be used as cost functions in systematic model parameter optimizations to quantify the observation-model misfit for LKF statistics. Based on our promising results for simulations with modified compressive and shear strength, we recommend to also consider the ellipse ratio *e* in such an optimization for (E)VP models.

SIREx is the first sea ice model intercomparison effort that includes the main three rheologies employed in state-of-the-art sea ice models. The (E)VP rheology, which is used in most, if not all, climate models, is well represented in our study with 30 out of 36 simulations. This emphasizes how important our findings are for improving the representation of sea ice conditions in climate projections. The MEB and EAP rheology, however, are each based on simulations of only one model code. For the MEB simulations, it is also not possible to cleanly disentangle whether it is the rheology alone or the combination with the Lagrangian grid that makes neXtSIM scoring higher than the average (E) VP model at the same resolution. A more diverse set of simulations or ideally a single modeling framework for all three rheologies would be required to unambiguously assess the effect of the rheology.

The combination of the deformation statistics from SIREx Part I and the LKF statistics from SIREx Part II shows that only the realism of deformation rate PDFs and spatial scaling analysis are correlated to LKF density skill metrics. This suggests that PDFs of deformation rates, which are relatively easy to compute, can be used to quickly assess sea ice simulations, for instance during model tuning, before applying more sophisticated metrics such as the LKF statistics of Part II or the scaling analysis of Part I. There is, however, no clear correlation between the temporal multifractal scaling parameters and the LKF statistics, so that the scaling metrics of sea ice deformation provide only some information on the number of deformation features, but that the feature-based evaluation of sea ice deformation used in this study is required to deduce other properties of deformation features.

Finally, our feature-based evaluation of LKFs uses sea ice deformation fields implicitly assuming that these deformation-field-based LKFs coincide with leads and pressure ridges. We find, however, that not all models develop the expected concentration and thickness anomalies along LKFs. In particular, the magnitude of these anomalies depends strongly on the model resolution. Only simulations with a resolution smaller than 2 km produce sea ice concentration down to 50% in grid cells tagged as “leads”. These differences due to model resolution could have a large impact on high-resolution coupled climate simulations to come, as ice-air-ocean interaction processes, such as heat flux or form drag and the extent to which they are resolved by an atmospheric model, depend fundamentally on the magnitude of the concentration and thickness along leads and pressure ridges. The evaluation of leads and pressure ridges as features of ice concentration and thickness fields will become feasible with increasing resolution of both simulations and satellite products. This will nicely complement the analysis of LKFs based exclusively on deformation data.
## Appendix A: Comparability of RGPS and Model Data

RGPS is a Lagrangian drift data set derived from SAR imagery. In the beginning of each winter, a network of virtual buoys is initialized on a regular grid with a grid spacing of 10 km (in coastal regions 25 km; Kwok, 1998). The position of each virtual buoy is updated every time a new SAR image is available, which covers the position of this virtual buoy. The timing of the individual position record is different for all buoys, as they were covered by different or multiple overflights. On average, the positions of all buoys are updated every 3 days. Deformation rates are computed from the virtual buoy positions using line integral approximations. The temporal and spatial scales associated with these deformation rates vary due to the irregular temporal sampling and distortion of the cells caused by the advection. The average temporal and spatial scale is 3 days and 10 km. The Lagrangian deformation data is interpolated onto a regular grid with a spacing of 12.5 km. This Eulerian deformation data set was used to derive the RGPS LKF data set used in our analysis (Hutter et al., 2019).

The principle of computing deformation rates from tracked displacements has interesting consequences with respect to the resolution of deformation features that are resolved in the data. All displacements originating from the fracture between two position records are recorded regardless of their duration, as long as they exceed the spatial resolution of the SAR images, which is roughly 100 m. Therefore, the deformation features do not have a specific temporal resolution. The spatial resolution of the SAR images limits only the magnitude of the deformation rates (Lindsay et al., 2003), but not the width of the deformation features causing the deformation. In the context of model-observation comparisons, it is reasonable to choose the spatial resolution of LKFs referring to their minimal width. As a minimal width of LKFs cannot be derived from RGPS, we refer to the RGPS data set having a 3-day temporal and 12.5 km spatial sampling rate, instead of using the term resolution.

The fact that the deformation derived from RGPS is not limited to a certain feature width nor a deformation duration time complicates the comparison with models that do not use a subgrid parameterization of small-scale deformation, but explicitly resolve deformation features. In our comparison, these are all (E)VP, and EAP models. In those models, 5–7 grid cells are needed to resolve an LKF as a discontinuity in the concentration and thickness fields (Williams & Tremblay, 2018). Thus, the effective resolution of these simulations is accordingly larger than model's grid spacing. This explains in part why the coarser resolution models in the comparison show fewer LKFs than RGPS (Figure 1). Since the effective resolution has not been quantified, adapting of the model-observation analysis to the effective resolution in our analysis is not possible.

The spatial sampling of RGPS does not allow to distinguish if the recorded displacement is caused by one or multiple deformation features. Very high resolution models could resolve multiple deformation features within a 12.5 × 12.5 RGPS cell, which would distort the model-observation comparison of LKF numbers. The highest resolution models in our comparison, however, have a resolution of 1–2 km. Thus, their effective resolution is in the range of the spatial sampling rate of RGPS and we can exclude that this affects our results.

The 3-day RGPS deformation rates can be seen as an “integral” of all deformation taking place within the 3-day sampling period of RGPS. We mimic this by flagging LKF pixels for each daily field and then combine three daily maps to a cumulative 3-day LKF map (for more details see Section 3.2). Assuming that the exponential tail of the LKF lifetime distribution (Figure 6) also holds for lifetimes smaller than 3 days, we potentially miss very short-lived LKFs in our model LKF data sets. However, LKFs are associated with strong velocity gradients and their signals are imprinted in the daily mean fields. In addition, forcing at a very short temporal scale is needed to create these short-lived deformations. Actually, we do not find differences in the deformation rate PDFs when using daily versus hourly snapshot or mean model diagnostics. Thus, we are confident that the different temporal sampling for models and RGPS does not have a significant effect on our analysis. Very high frequency model output could be used to quantify this effect in a dedicated study.

## Appendix B: Principles of Intersection Angle Selection

We study the intersection angle of LKFs because it is directly related to the fracture angle for conjugate faults. The fracture angle is the angle between the LKF and the direction of the main stress, therefore half of the intersection angle, and it provides insights into the fracture physics. To compute the fracture angle, the direction
<figure><img src="image_15.png" /><figcaption>**Figure B1.** Schematic of relative displacement, vorticity, stress direction, and intersection angle along two intersection linear kinematic features.</figcaption></figure>

of stress must be known. For RGPS and model simulations, however, only the drift and deformation data are provided. In our analysis, we use the vorticities derived from the relative displacement along two intersecting LKFs to determine the main direction of stress. The vorticity is computed from the sea ice drift fields. The two LKFs of a conjugate fault have vorticities of opposite sign (Figure B1). Thus, we label intersecting pairs of LKFs with vorticities of opposite sign as conjugate faults and compute the intersection angle relative to stress direction by measuring the angle from the LKF with positive vorticity to the LKF with negative vorticity. If both intersecting LKFs have vorticities of the same sign, it is not possible to determine the main stress direction causing the deformation. In this case, we take both possible intersection angles (angle measured from the LKF with positive vorticity to LKF with negative vorticity and vice versa) into account for the computation of the intersection angle PDF (Figure 5).

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

The Table C1 details the meaning of most acronyms relating to models, configurations, and reanalyzes used in the runs participating in this study.

<caption>**Table C1**</caption>

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

<table><thead><tr><th>Acronym</th><th>Meaning</th></tr></thead><tbody><tr><td><em>Models</em></td><td></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></tbody></table>
<table><caption>Table C1</caption><caption>Continued</caption><thead><tr><th>Acronym</th><th>Meaning</th></tr></thead><tbody><tr><td colspan="2"><em>Configurations</em></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"><em>Reanalyses</em></td></tr><tr><td>JRA55</td><td>Japanese 55-year Reanalysis</td></tr><tr><td>ERA-Interim</td><td>European Centre for Medium-Range Weather Forecasts (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>

Acknowledgments

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. T. Rasmussen was funded by the Danish State through the National Centre for Climate Research and by the SALIENSEAS project part of the ERA4CS program, which is co-funded by the Innovation Fund Denmark and the Horizon 2020 Framework Program of the European Union (Grant 690462). Q. Wang is supported by the German Helmholtz Climate Initiative REKLIM (Regional Climate Change). Data is available upon request. 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). 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. 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 Foundation 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 Modernization Program. N. Koldunov was supported by the project S1 (Diagnosis and Metrics in Climate Models) of the Collaborative Research Centre TRR 181 “Energy Transfer in Atmosphere and Ocean” funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Project 274762653. We thank the two anonymous reviewers whose constructive comments helped improve and clarify this manuscript. Open access funding enabled and organized by Projekt DEAL.

Data Availability Statement

All published code and data products originating from the Sea Ice Rheology Experiment is organized in the SIREx data collection: https://zenodo.org/communities/sirex: this includes the model output in netCDF format (Bouchat et al., 2022b), the model Lagrangian trajectories and deformation data (Bouchat et al., 2022c), and the RGPS Lagrangian composite deformation data (Bouchat & Hutter, 2022) as well as the LKF data (Hutter et al., 2022). The code of the LKF detection and tracking algorithm is available on github: https://github.com/ nhutter/lkf_tools.git (Hutter, 2019). The RGPS LKF data set is available on PANGAEA: https://doi.org/10.1594/ PANGAEA.898114 (Hutter et al., 2019b). 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/. For more details on the ANHA configuration, visit http://knossos.eas.ualberta.ca/anha/anhatable.php.

References

Blockley, E., Vancoppenolle, M., Hunke, E., Bitz, C., Feltham, D., Lemieux, J.-F., et al. (2020). The future of sea ice modeling: Where do we go from here? *Bulletin of the American Meteorological Society*, 101(8), E1304–E1311. https://doi.org/10.1175/bams-d-20-0073.1

Bouchat, A., Hutter, N., Chanut, J., Dupont, F., Dukhovskoy, D., Garric, G., et al. (2022a). 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

Bouchat, A., Hutter, N., Chanut, J., Dupont, F., Dukhovskoy, D., Garric, G., et al. (2022b). 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. (2022c). Model Lagrangian trajectories and deformation data analyzed in the Sea Ice Rheology Experiment—Part I (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]. 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(7), 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 spatiotemporal scaling statistics. *Journal of Geophysical Research: Oceans*, 125(8). https://doi.org/10.1029/2019JC015944

Bouillon, S., & Rampal, P. (2015). Presentation of the dynamical core of neXtSIM, a new sea ice model. *Ocean Modelling*, 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

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 Modelling*, **120**, 60–78. https://doi.org/10.1016/j.ocemod.2017.10.007

Cunningham, G. F., Kwok, R., & Banfield, J. (1994). Ice lead orientation characteristics in the Winter Beaufort Sea. In *Proceedings of IGARSS '94—1994 IEEE international geoscience and remote sensing symposium* (Vol. 3, pp. 1747–1749). https://doi.org/10.1109/IGARSS.1994.399553

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

Dempsey, J., Xie, Y., Adamson, R., & Farmer, D. (2012). Fracture of a ridged multi-year Arctic Sea ice floe. *Cold Regions Science and Technology*, **76–77**, 63–68. https://doi.org/10.1016/j.coldregions.2011.09.012

Dukhovskoy, D. S., Yashayaev, I., Proshutinsky, A., Bamber, J. L., Bashmachnikov, I. L., Chassignet, E. P., & Tedstone, A. J. (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

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

Erlingsson, B. (1988). Two-dimensional deformation patterns in sea ice. *Journal of Glaciology*, **34**(118), 301–308. https://doi.org/10.1017/s0022143000007061

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. *Annals of Glaciology*, **52**(57), 123–132. https://doi.org/10.3189/172756411795931499

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. *Deep Sea Research Part II: Topical Studies in Oceanography*, **55**(8), 943–962. (Ice Station POLarstern (ISPOL): Results of interdisciplinary studies on a drifting ice floe in the western Weddell Sea). https://doi.org/10.1016/j.dsr2.2007.12.026

Heorton, H. D. B. S., Feltham, D. L., & Tsamados, M. (2018). Stress and deformation characteristics of sea ice in a high-resolution, anisotropic sea ice model. *Philosophical Transactions of the Royal Society A: Mathematical, Physical & Engineering Sciences*, **376**(2129), 20170349. https://doi.org/10.1098/rsta.2017.0349

Hibler, W. D. (1979). A dynamic thermodynamic sea ice model. *Journal of Physical Oceanography*, **9**, 815–846. https://doi.org/10.1175/1520-0485(1979)009<0815:adtsim>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., Allard, R., Blain, P., Blockley, E., Feltham, D., Fichefet, T., et al. (2020). Should sea-ice modeling tools designed for climate research be used for short-term forecasting? *Current Climate Change Reports*, **1**, 1–16.

Hutchings, J. K., Heil, P., & Hibler, W. D. (2005). Modeling linear kinematic features in sea ice. *Monthly Weather Review*, **133**, 3481–3497. https://doi.org/10.1175/MWR3045.1

Hutter, N. (2015). *Viscous-plastic sea-ice models at very high resolution Master's Thesis*, University of Bremen, Alfred Wegener Institute, Helmholtz Centre for Polar and Marine research. https://doi.org/10.10013/epic.46129

Hutter, N. (2019). *Ikf_toolS: A code to detect and track Linear Kinematic Features (LKFs) in sea-ice deformation data (Version v1.0) [code]*.
Zenodo. https://doi.org/10.5281/zenodo.2560078

Hutter, N., Bouchat, A., Dupont, F., Dukhovskoy, D., Koldunov, N., Lee, Y., et al. (2022). Linear kinematic feature detected and tracked in sea-ice deformation simulation by all models participating in the Sea Ice Rheology Experiment and from RGPS (1.0) [Data set]. 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. (2019a). 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

Hutter, N., Zampieri, L., & Losch, M. (2019b). Linear Kinematic Features (leads & pressure ridges) detected and tracked in RADARSAT Geophysical Processor System (RGPS) sea-ice deformation data from 1997 to 2008 [Data set]. Alfred Wegener Institute, Helmholtz Centre for Polar and Marine Research, Bremerhaven, PANGAEA. https://doi.org/10.1594/PANGAEA.898114

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*. https://doi.org/10.1029/2018MS001485

Kwok, R. (1998). The RADARSAT geophysical processor system. In *Analysis of SAR data of the polar oceans* (pp. 235–257). Springer Berlin Heidelberg. 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. *Iutam symposium on scaling laws in ice mechanics and ice dynamics*, In J. Dempsey, & H. Shen (Eds.), (Vol. 94, pp. 315–322). Springer Netherlands. https://doi.org/10.1007/978-94-015-9735-7_26

Kwok, R., Hunke, E. C., Maslowski, W., Menemenlis, D., & Zhang, J. (2008). Variability of sea ice simulations assessed with RGPS kinematics. *Journal of Geophysical Research*, **113**(C11), 1–20. https://doi.org/10.1029/2008JC004783

Lemieux, J.-F., & Tremblay, B. (2009). Numerical convergence of viscous-plastic sea ice models. *Journal of Geophysical Research*, **114**(C5). https://doi.org/10.1029/2008jc005017

Lindsay, R. W., Zhang, J., & Rothrock, D. A. (2003). Sea-ice deformation rates from satellite measurements and in a model. *Atmosphere-Ocean*, **41**(1), 35–47. https://doi.org/10.3137/ao.410103

Linow, S., & Dierking, W. (2017). Object-based detection of linear kinematic features in sea ice. *Remote Sensing*, **9**(5). https://doi.org/10.3390/rs9050493

Losch, M., Fuchs, A., Lemieux, J.-F., & Vanselow, A. (2014). A parallel Jacobian-free Newton-Krylov solver for a coupled sea ice-ocean model. *Journal of Computational Physics Part A*, **257**, 901–911. https://doi.org/10.1016/j.jcp.2013.09.026

Lüpkes, C., Gryanik, V. M., Witha, B., Gryschka, M., Raasch, S., & Gollnik, T. (2008). Modeling convection over arctic leads with les and a non-eddy-resolving microscale model. *Journal of Geophysical Research*, **113**(C9), 1–17. https://doi.org/10.1029/2007JC004099

Madsen, K. S., Rasmussen, T. A. S., Ribergard, M. H., & Ringgaard, I. M. (2016). *High resolution sea ice modelling and validation of the Arctic with focus on South Greenlandic waters, 2004–2013*. Polarforschung. 85.
Mahoney, A. R., Eicken, H., Shapiro, L. H., Heinrichs, T., Meyer, F. J., & Gaylord, A. G. (2012). *Mapping and characterization of recurring spring leads and landfast ice in the Beaufort and Chukchi Seas* [Report]. University of Alaska, Coastal Marine Institute. (Final Report: OCS Study BOEM 2012-067).

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., Métaxian, J.-P., Grangeon, J., Roux, P.-F., & Haapala, J. (2011). Low-frequency bursts of horizontally polarized waves in the Arctic Sea-ice cover. *Journal of Glaciology*, 57(202), 231–237. https://doi.org/10.3189/002214311796405834

Martin, T., Tsamados, M., Schroeder, D., & Feltham, D. L. (2016). The impact of variable Sea ice roughness on changes in Arctic Ocean surface stress: A model study. *Journal of Geophysical Research: Oceans*, 121(3), 1931–1952. https://doi.org/10.1002/2015jc011186

McPhee, M. G., Kwok, R., Robins, R., & Coon, M. (2005). Upwelling of Arctic pycnocline associated with shear motion of sea ice. *Geophysical Research Letters*, 32(10), 2004g021819. https://doi.org/10.1029/2004g021819

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, https://doi.org/10.1038/s41598-018-24660-0

Muilwijk, M., Ilicak, M., Cornish, S. B., Danilov, S., Gelderloos, R., Gerdes, R., & Wang, Q. (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

Ólason, E., Rampal, P., & Dansereau, V. (2021). On the statistical properties of sea-ice lead fraction and heat fluxes in the Arctic. *The Cryosphere*, 15, 1053–1064. https://doi.org/10.5194/tc-15-1053-2021

Proshutinsky, A., Krishfield, R., & Timmermans, M.-L. (2020). Introduction to special collection on Arctic Ocean modeling and observational synthesis (FAMOS) 2: Beaufort gyre phenomenon. *Journal of Geophysical Research: Oceans*, 125(2), e2019JC015400. https://doi.org/10.1029/2019jc015400

Rampal, P., Bouillon, S., Ólason, 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., & Samaké, A. (2019). On the multi-fractal scaling properties of sea ice deformation. *The Cryosphere Discussions*, 1–45. https://doi.org/10.5194/tc-13-2457-2019

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. (2020). Non-normal flow rules affect fracture angles in sea ice viscous-plastic rheologies. *The Cryosphere Discussions*, 1–24. https://doi.org/10.5194/tc-2020-153

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

Schulson, E., Fortt, A., Iliescu, D., & Renshaw, C. (2006). On the role of frictional sliding in the compressive fracture of ice and granite: Terminal vs. post-terminal failure. *Acta Materialia*, 54(15), 3923–3932. https://doi.org/10.1016/j.actamat.2006.04.024

Schulson, E. M., & Hibler, W. D. (2004). Fracture of the Winter sea ice cover on the Arctic Ocean. *Comptes Rendus Physique*, 5(7), 753–767. https://doi.org/10.1016/j.crhy.2004.06.001

Spreen, 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

Stamoulis, C., & Dyer, I. (2000). Acoustically derived ice-fracture velocity in central arctic pack ice. *Journal of the Acoustical Society of America*, 108(1), 96–104. https://doi.org/10.1121/1.429448

Sumata, H., Kauker, F., Karcher, M., & Gerdes, R. (2019). Simultaneous parameter optimization of an Arctic Sea ice–ocean model by a genetic algorithm. *Monthly Weather Review*, 147(6), 1899–1926. https://doi.org/10.1175/mwr-d-18-0360.1

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*. https://doi.org/10.1002/2016JC012128

Walter, B. A., & Overland, J. E. (1993). The response of lead patterns in the Beaufort Sea to storm-scale wind forcing. *Annals of Glaciology*, 17, 219–226. https://doi.org/10.3189/S0260305500012878

Wang, K. (2007). Observing the yield curve of compacted pack ice. *Journal of Geophysical Research: Oceans*, 112(C5). https://doi.org/10.1029/2006JC003610

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, 4703–4725. https://doi.org/10.5194/tc-15-4703-2021

Wang, Q., Koldunov, N. V., Danilov, S., Sidorenko, D., Wekerle, C., Scholz, P., & Jung, T. (2020). Eddy kinetic energy in the Arctic Ocean from a global simulation with a 1-km Arctic. *Geophysical Research Letters*, 47(14), e2020GL088550. https://doi.org/10.1029/2020GL088550

Wernecke, A., & Kaleschke, L. (2015). Lead detection in Arctic Sea ice from Cryosat-2: Quality assessment, lead area fraction and width distribution. *The Cryosphere*, 9(5), 1955–1968. https://doi.org/10.5194/tc-9-1955-2015

Williams, J., & Tremblay, L. B. (2018). The dependence of energy dissipation on spatial resolution in a viscous-plastic sea-ice model. *Ocean Modelling*, 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). https://doi.org/10.3390/rs8010004