**Research article**
10 Dec 2021

**Research article** | 10 Dec 2021

# A generalized stress correction scheme for the Maxwell elasto-brittle rheology: impact on the fracture angles and deformations

Mathieu Plante and L. Bruno Tremblay

**Mathieu Plante and L. Bruno Tremblay**Mathieu Plante and L. Bruno Tremblay

- Department of Atmospheric and Oceanic Sciences, McGill University, Montréal, Québec, Canada

- Department of Atmospheric and Oceanic Sciences, McGill University, Montréal, Québec, Canada

**Correspondence**: Mathieu Plante (mathieu.plante@mail.mcgill.ca)

**Correspondence**: Mathieu Plante (mathieu.plante@mail.mcgill.ca)

Received: 06 Dec 2020 – Discussion started: 02 Feb 2021 – Revised: 16 Sep 2021 – Accepted: 01 Nov 2021 – Published: 10 Dec 2021

The Maxwell elasto-brittle (MEB) rheology uses a damage parameterization to represent the brittle fracture of sea ice without involving plastic laws
to constrain the sea ice deformations. The conventional MEB damage parameterization is based on a correction of super-critical stresses that binds
the simulated stress to the yield criterion but leads to a growth of errors in the stress field. A generalized damage parameterization is developed
to reduce this error growth and to investigate the influence of the super-critical stress correction scheme on the simulated sea ice fractures,
deformations and orientation of linear kinematic features (LKFs). A decohesive stress tensor is used to correct the super-critical stresses towards
different points on the yield curve. The sensitivity of the simulated sea ice fractures and deformations to the decohesive stress tensor is
investigated in uniaxial compression experiments. Results show that the decohesive stress tensor influences the growth of residual errors associated
with the correction of super-critical stresses, the orientation of the lines of fracture and the short-term deformation associated with the damage,
but it does not influence the long-term post-fracture sea ice deformations. We show that when ice fractures, divergence first occurs while the elastic
response is dominant, and convergence develops post-fracture in the long term when the viscous response dominates – contrary to laboratory
experiments of granular flow and satellite imagery in the Arctic. The post-fracture deformations are shown to be dissociated from the fracture
process itself, an important difference with classical viscous plastic (VP) models in which large deformations are governed by associative plastic
laws. Using the generalized damage parameterization together with a stress correction path normal to the yield curve reduces the growth of errors
sufficiently for the production of longer-term simulations, with the added benefit of bringing the simulated LKF intersection half-angles closer to
observations (from 40–50 to 35–45^{∘}, compared to 15–25^{∘} in observations).

Sea ice is a thin layer of solid material that insulates the polar oceans from the cold atmosphere. When sea ice fractures and a lead opens, large heat and moisture fluxes take place between the ocean and the atmosphere, significantly affecting the polar meteorology on short timescales and the climate system on long timescales (Maykut, 1982; Ledley, 1988; Lüpkes et al., 2008; Li et al., 2020). The refreezing of leads contributes to the sea ice mass balance (Wilchinsky et al., 2015; Itkin et al., 2018); the associated brine rejection drives the thermohaline ocean circulation in the Arctic and vertical eddies in the ocean mixed layer (Kozo, 1983; Matsumura and Hasumi, 2008). As such, the production of accurate seasonal-to-decadal projections using coupled models requires an accurate representation of sea ice deformations along linear kinematic features (LKFs).

As sea ice models are moving to higher spatial resolutions, they become increasingly capable of resolving LKFs (Hutter et al., 2018, 2021). The representation of smaller-scale fracture physics on the other hand yet remains a challenge, as most sea ice models are based on a continuum approach and rely on parameterizations to relate sea ice deformations to unresolved fractures. To this day, this is most commonly done using plastic rheologies or modifications thereof (Hibler, 1979; Hunke and Dukowicz, 1997), which have benefited from improved numerical schemes and efficiency to solve the highly non-linear momentum equation (Lemieux et al., 2008, 2014; Kimmritz et al., 2016; Koldunov et al., 2019). These models use plastic flow rules to represent the rate invariance of sea ice deformations at large spatio-temporal scales, in which the sea ice can be considered ductile, but neglect the influence of the smaller-scale physics associated with the brittle fractures. A number of other rheologies have been developed over the years to relate the sea ice deformations to the smaller-scale fracture physics (Tremblay and Mysak, 1997; Wilchinsky and Feltham, 2004; Schreyer et al., 2006; Sulsky and Peterson, 2011; Rampal et al., 2016; Dansereau et al., 2016; Damsgaard et al., 2018). This brings a diversity of sea ice rheologies, with different physical and numerical frameworks influencing the representation of sea ice deformations at different scales.

The Sea Ice Rheology Experiment (SIREx; Bouchat et al., 2021; Hutter et al., 2021), a coordinated effort between several ice–ocean modelling groups, assessed
the pan-Arctic sea ice deformation statistics simulated by different sea ice rheologies. SIREx included the classical viscous–plastic
(Hibler, 1979) and elastic–viscous–plastic (Hunke and Dukowicz, 1997) sea ice rheologies as well as the elastic–anisotropic (Wilchinsky and Feltham, 2004) and
Maxwell elasto-brittle (MEB; Dansereau et al., 2016) rheologies that include parameterizations of unresolved small-scale physics. All participating
sea ice models produced sea ice deformation characteristics that have previously been associated with brittle behaviour, such as the scaling and
spatio-temporal coupling of sea ice deformations (Bouchat et al., 2021), when run at sufficiently high resolution. The extent at which the inclusion of
smaller-scale fracture physics improves this brittle behaviour thus remains an open question. Additionally, all rheologies produce similar angles
between conjugate pairs of LKFs, a measure usually intimately related to the fracture mechanics and shear strength of a material
(Bardet, 1991; Wang, 2007), showing a peek probability at 90^{∘} while the observed angles are in the range of 30–50^{∘}
(Hutter et al., 2021). This calls for the improvement of sea ice rheological models, such as modifications of the mechanical strength parameters and
yield curve (Bouchat and Tremblay, 2017; Ringeisen et al., 2019; Dansereau et al., 2019), the use of non-associated flow rules (in the case of classical plastic models;
Ringeisen et al., 2021) or modifications of fine-scale fracture parameters (in the case of the elastic anisotropic plastic (EAP) and MEB rheologies).

In the Maxwell elasto-brittle (MEB) rheology (Dansereau et al., 2016), the smaller-scale fracture physics is represented by a damage parameterization that was derived for rock mechanics and seismic models (Amitrano et al., 1999; Amitrano and Helmstetter, 2006) and adapted for the large-scale modelling of sea ice (Girard et al., 2011; Bouillon and Rampal, 2015; Rampal et al., 2016). This parameterization aims at representing the brittle character of sea ice by using a damage parameter to represent the changes in material properties associated with fractures. This differs from parameterizations used in viscous plastic models in that the large-scale sea ice deformations are not governed by plastic or granular flow rules. Instead, the sea ice deformations in the MEB model are preconditioned by the presence of damage, and the development of LKFs is associated with the far-field stress concentration response to local damage, leading to the propagation of the damage (i.e. fractures) in space (Dansereau et al., 2019). While still based on the continuum assumption, it allows for brittle fractures to influence the sea ice dynamics over shorter timescales. It is currently used in the large-scale sea ice finite element model neXtSIM (Rampal et al., 2019) and a finite difference version was recently implemented in the McGill Sea Ice Model Version 5 (McGill SIM5) (Plante et al., 2020).

With the MEB rheology being relatively new, the extent to which the sea ice deformations are sensitive to the numerical and material strength parameters has not been thoroughly tested yet. Nonetheless, the orientation of the simulated faults in uniaxial compression experiments is known to be sensitive to the angle of internal friction and to the Poisson ratio (Dansereau et al., 2019). This sensitivity is attributed to the influence of these parameters on the far-field stress concentration response to local damage, which determines the direction of the damage propagation. This suggests that the simulated angle of fracture may be sensitive to the exact choice of damage parameterization, but has not yet been tested. Additionally, while the neXtSIM model performed well compared to other SIREx models, its different numerics (e.g. Lagrangian scheme with a triangular adaptive mesh) could also be responsible for the different scaling and localization statistics (Bouchat et al., 2021). The finite difference implementation of the MEB rheology in the McGill SIM5 model, on the other hand, shows fast growth of residual errors at the grid scale – in ideal experiments – that significantly affect the post-fracture sea ice deformations (Plante et al., 2020). These errors result from the stress correction scheme used in the MEB rheology to define the growth of damage and to bring super-critical stresses back to the yield curve. To our knowledge, defining the damage in terms of the super-critical stress correction is new and unique to the EB and MEB sea ice rheologies. For instance, many progressive damage models instead represent the damage parameter as a discrete function of the number of failure cycles (Amitrano and Helmstetter, 2006; Carrier et al., 2015). In continuum damage mechanics, the damage parameter is derived instead from thermodynamic laws (Murakami, 2012) to simulate material fatigue. In the elastic–decohesive (ED) rheology, material damage is not parameterized but a decohesive strain rate explicitly represents the material discontinuity associated with the ice fracture and reduces the material strength of sea ice, based on the orientation of the failure surface (Schreyer et al., 2006; Sulsky and Peterson, 2011).

In this paper, we present a generalization of the damage parameterization in which a decohesive stress tensor is introduced in the stress correction scheme such that the super-critical stresses can be brought back to the yield curve following different stress correction paths in the stress invariant space. The generalization is used to reduce the growth of the residual errors associated with the stress correction and tested in uniaxial loading experiments to examine the influence of the stress correction on the simulated sea ice fracture and deformations. The sensitivity of the simulated fracture angles to the decohesive stress tensor is also investigated to find the stress correction paths that present the added benefit of bringing the simulated fracture angles closer to observations.

This paper is organized as follows. In Sect. 2, we present the MEB rheology and governing equations. The generalized stress correction scheme is described in Sect. 3. The uniaxial loading experiment set-up is presented in Sect. 4 along with the definition of diagnostics used to quantify the growth of damage and of residual errors. Results are presented in Sect. 5, with a focus on the material behaviour in uniaxial compression experiments and its response to the changes in the damage parameterization. In Sect. 6, we provide a discussion on the generalized damage parameterization performance and other model sensitivities. Conclusions are summarized in Sect. 7.

## 2.1 Momentum and continuity equations

The simulations are run using the MEB model implemented on an Eulerian, finite difference Arakawa C grid in the McGill SIM5 (Tremblay and Mysak, 1997; Lemieux et al., 2008; Plante et al., 2020). The vertically integrated 2D momentum equation for sea ice can be written as (ignoring the sea surface tilt, the Coriolis term and the ice grounding terms),

where *ρ*_{i} is the ice density, *h* is the mean ice thickness, ** u** ($=u\widehat{i}+v\widehat{j}$) is the ice velocity
vector,

*σ*is the vertically integrated internal stress tensor and

*τ*is the net external surface stress from winds and ocean currents. This simplified formulation is appropriate for short-term uniaxial loading experiments but can result in small errors in ice velocity when using a realistic model domain and forcing (Turnbull et al., 2017). Following Plante et al. (2020), we define the uniaxial loading by a surface wind stress

*τ*

_{a}and prescribe an ocean at rest below the ice:

where *ρ*_{w} is the water density, *C*_{dw} is the water drag coefficient and ** u** is the sea ice velocity (see values in
Table 1).

The prognostic equations for the mean ice thickness *h* (volume per grid cell area) and concentration *A* are written as

where the thermodynamic source and sink terms are ignored.

## 2.2 Maxwell elasto-brittle rheology

The MEB model differs from classical sea ice models in that it represents the brittle character of sea ice using a damage parameter to represent the effect of local fracture on the large-scale sea ice material properties. The sea ice deformations in the MEB model thus occur post-fracture, rather than simultaneously as in most sea ice models using granular or plastic flow laws, and the formation of LKFs follows from the propagation of damage in space over short timescales during the fracture process.

In the MEB rheology, the ice behaves as a visco-elastic material with a fast elastic response to forcing and a slower viscous response that acts over a longer timescale. The governing equation for this visco-elastic material can be written as (Dansereau et al., 2016, 2017; Plante et al., 2020)

where *E* is the elastic stiffness defined as the vertically integrated Young modulus of sea ice, *λ* is the viscous relaxation timescale,
**C** is the (fourth-order) elastic tensor, : denotes the inner double tensor product and $\dot{\mathit{\u03f5}}$ is the (second-order) strain rate
tensor. The tensors **C** and $\dot{\mathit{\u03f5}}$ on the right-hand side of Eq. (5) can be written in matrix form by representing the
three independent components of the stress and strain tensors in a vector (see Rice, 2010) and the nine independent components of the elastic
modulus tensor in a 3 × 3 matrix as

where *ν* (= 0.33) is the Poisson ratio, which defines the relative amount of deformation on the plane parallel to the loading.

The relative importance of the elastic and viscous components (first and second terms on the left-hand side in Eq. 5) are
determined by the magnitude of the elastic modulus *E* and viscous relaxation timescale *λ*. *E* and *λ* are functions of the ice
thickness, concentration and damage, such that the elastic term dominates when the ice is undamaged while the viscous term dominates when the ice is
heavily fractured. The elastic modulus *E* and viscous relaxation timescale *λ* are written as

where *Y* (= 1 GPa) is the Young modulus of undeformed sea ice, *d* is the damage parameter ($\mathrm{0}<d<\mathrm{1}$), *a* (= 20) is the standard ice
concentration parameter (Hibler, 1979; Rampal et al., 2016), *λ*_{0} (= 10^{5} s, ≈ 1 d) is the viscous relaxation timescale for undamaged sea ice and *α* is a parameter defining the post-fracture transition to the viscous regime. This damage-based transition to
post-fracture viscosity represents a simplification of the observed plasticity (rate independence) of sea ice deformations
(Coon et al., 1974; Tuhkuri and Lensu, 2002).

## 2.3 Yield criterion

Damage (or fracture) occurs when the internal stress state exceeds the Mohr–Coulomb failure criterion,

where

where *σ*_{I} is the isotropic normal stress invariant (compression defined as negative), *σ*_{II} is the maximum shear
stress invariant, (*σ*_{xx},*σ*_{yy},*σ*_{xy}) are the components of the stress tensor, *μ* (= sin *ϕ*) is the coefficient of
internal friction of sea ice, *ϕ* (= 45^{∘}) is the angle of internal friction and *c* is the vertically integrated cohesion, defined as

where *c*_{0} (= 10 kN m^{−2}) is the cohesion of sea ice derived from observations (Sodhi, 1977; Tremblay and Hakakian, 2006; Plante et al., 2020) or laboratory
experiments (Timco and Weeks, 2010). No compressive or tensile strength cut-off is used in this analysis. The reader is referred to
Table 1 for a list of default model parameters.

## 2.4 Damage parameterization

The prognostic equation for the damage parameter *d* in the standard MEB rheology is parameterized using a relaxation term with timescale
*T*_{d} (= 1 s) as

where

is a damage factor ($\mathrm{0}<\mathrm{\Psi}<\mathrm{1}$), *σ*_{c} is the critical stress lying on the yield curve and *σ*^{′} is the uncorrected
stress state lying outside of the yield curve. Thermodynamic healing and the advection of damage are neglected as we are focusing on the ice fracture,
which occurs at a timescale (seconds) much shorter than the healing and advection timescales (hours). Adding these terms does not change the results
and conclusions presented in this paper but increases the localization of the ice fractures with higher damage values that in turn increases
ridging. These terms should be included in longer-term integration of the MEB model.

When the ice fractures, the damage factor Ψ is used to scale the super-critical stresses back towards the yield curve. The prognostic equation
for the temporal evolution of the super-critical stress tensor *σ*^{′} is written as a relaxation equation of the same form as in
Eq. (14):

This stress correction scheme corresponds to scaling all the individual stress components by the factor Ψ, such that the stress state is corrected back onto the yield curve in the stress invariant space by following a line passing through the origin. This results in a dependency of the stress correction magnitude and of the damage on the super-critical stress state; i.e., the stress correction path becomes increasingly parallel to the yield curve for increasing compressive super-critical stresses, which also increases the numerical errors (Plante et al., 2020). We hereafter refer to this scheme as the “standard stress correction”.

We propose a generalized damage parameterization where the super-critical stresses are corrected back to the yield curve along a line oriented at any
angle *γ* from the *y* axis in the stress invariant space (see Fig. 1). This generalization is developed with the goal of reducing the
growth rate of the numerical errors in the MEB model by removing the dependency of the stress correction path on the super-critical stress state,
while keeping the changes in the damage parameterization to a minimum so that it can be easily added to other MEB model implementations (and other
damage-based models). In the MEB model, the exact path along which the super-critical stresses is returned to the yield curve is not known a priori,
as the stress state never exceeds the yield criterion in reality. The proposed generalization allows us to investigate the influence of the
super-critical stress correction path angle on the simulated fractures and deformations. Other physically meaningful modifications of the stress
correction that are based on thermodynamics principles are left for future work (see for instance Murakami, 2012).

We define the damage factor in the generalized damage parameterization in terms of the shear stress invariant only as

where *σ*_{IIc} is the critical shear stress invariant. The equation defining the stress correction path with angle *γ* (see
Fig. 1) can be written as

where *B* (= ${\mathit{\sigma}}_{\text{II}}^{\prime}-\mathrm{1}/\mathrm{tan}\left(\mathit{\gamma}\right){\mathit{\sigma}}_{\mathrm{I}}^{\prime}$) is defined from the super-critical stress state
(*σ*^{′}). The critical shear stress invariant (*σ*_{IIc}) is then defined as the intersection point between the yield curve
(Eq. 10) and the stress correction path (Eq. 18),

The damage factor can then be written in terms of the super-critical stress state invariants (${\mathit{\sigma}}_{\mathrm{I}}^{\prime}$,
${\mathit{\sigma}}_{\text{II}}^{\prime}$), the correction path angle *γ* and the coefficient of internal friction *μ* as

In this manner, the correction of super-critical stresses can follow any path in the stress invariant space provided that the damage increases when
ice fractures (Ψ<1 or *γ* < 90^{∘}). This formulation can also be used with a yield curve with zero isotropic tensile strength
(i.e. *c* = 0 kN m^{−1}), as opposed to the standard parameterization in which case any super-critical stress state is returned to the
origin (see Eq. 15 when *c* = 0 N m^{−1}).

Note that using a stress correction path other than the standard path to the origin means that the corrected normal stress differs from the scaled super-critical stress $\mathrm{\Psi}{\mathit{\sigma}}_{\mathrm{I}}^{\prime}$. We define this difference as the decohesive stress tensor (see Fig. 1), which is added to the damage parameterization to keep the corrected stress state on a given stress correction path. This effectively changes the stress correction while keeping the scalar definition of the damage parameter. The stress correction equation (Eq. 16) in the generalized damage parameterization then becomes

and the invariants of the decohesive stress tensor (*σ*_{ID} and *σ*_{IID}) are now defined as

When $\mathrm{tan}\mathit{\gamma}={\mathit{\sigma}}_{\mathrm{I}}^{\prime}/{\mathit{\sigma}}_{\text{II}}^{\prime}$ and ${\mathit{\sigma}}_{\text{ID}}={\mathit{\sigma}}_{\text{IID}}=\mathrm{0}$, we obtain the standard damage parameterization of Dansereau et al. (2016).

Note that the decohesive stress tensor used in this parameterization has a similar role as the decohesive strain rates used in the elastic–decohesive
model (Schreyer et al., 2006). In Schreyer et al. (2006), the decohesive strain represents the discontinuity in sea ice displacement associated with a
fracture and relaxes the effective stress rates. It is derived from a decohesion function that depends on the mode of failure. Here, we do not define
the strain discontinuity associated with the fractures but use the decohesive stress tensor *σ*_{D} to prescribe the orientation at
which the stress state is relaxed back onto the yield curve. This only indirectly influences the local strain rate via the constitutive equation.

## 3.1 Projected error

The error *δ*Ψ on the damage factor $\mathrm{\Psi}({\mathit{\sigma}}_{\mathrm{I}}^{\prime},{\mathit{\sigma}}_{\text{II}}^{\prime})$ can be written as (Plante et al., 2020)

where $\mathit{\delta}{\mathit{\sigma}}_{\mathrm{I}}^{\prime}$ and $\mathit{\delta}{\mathit{\sigma}}_{\text{II}}^{\prime}$ are the errors of the calculated stress invariants. Using
Eq. (21) and re-writing $\mathit{\delta}{\mathit{\sigma}}_{\mathrm{I}}^{\prime}$ and $\mathit{\delta}{\mathit{\sigma}}_{\text{II}}^{\prime}$ in terms of the
relative error *ϵ* (i.e., $\mathit{\delta}{\mathit{\sigma}}_{\mathrm{I}}^{\prime}=\mathit{\u03f5}{\mathit{\sigma}}_{\mathrm{I}}^{\prime}$, $\mathit{\delta}{\mathit{\sigma}}_{\text{II}}^{\prime}=\mathit{\u03f5}{\mathit{\sigma}}_{\text{II}}^{\prime}$), we obtain

where *R* is the error amplification ratio.

Given that the uncorrected stress is close to the yield criterion (i.e. ${\mathit{\sigma}}_{\text{II}}^{\prime}+\mathit{\mu}{\mathit{\sigma}}_{\mathrm{I}}^{\prime}-c\sim \mathrm{0}$),
the error amplification ratio *R* tends to infinity for

which corresponds to a path that runs parallel to the yield curve. This result is consistent with the instabilities in the standard stress correction
scheme during ridging reported in Plante et al. (2020), given that a line passing through the origin is nearly parallel to the Mohr–Coulomb yield curve
for large compressive stresses. In contrast, the path that maximizes the denominator (smallest error growth) has *γ* = 90^{∘}. This
path, however, corresponds to Ψ=1 and does not create damage. The possible stress correction path angles *γ* thus lie in the range
$\mathrm{arctan}(-\mathrm{1}/\mathit{\mu})<\mathit{\theta}$ < 90^{∘}.

Note that the error amplification ratio *R* is small for ${\mathit{\sigma}}_{\mathrm{I}}^{\prime}<\mathrm{0}$ but becomes infinitely large at the yield curve tip when
${\mathit{\sigma}}_{\text{II}}^{\prime}$ approaches 0 (see Eq. 25). This behaviour is opposite to that of the standard stress correction
scheme, which has small *R* values in tension and large values in compression (Plante et al., 2020). For this reason, we use both schemes
(i.e. Eq. 20 in compression and Eq. 15 in tension; see Fig. 1b) and set the transition between
the two schemes at the points where their paths are the same (i.e., at ${\mathit{\sigma}}_{\mathrm{I}}^{\prime}/{\mathit{\sigma}}_{\text{II}}^{\prime}=\mathrm{tan}\mathit{\gamma}$,
green line in Fig. 1b). The damage factor is then defined as

## 4.1 Experiment setup

We test the numerical and material behaviour of the MEB model and the generalized damage parameterization in uniaxial compression
experiments. Uniaxial experiments are designed to present conditions similar to those in laboratory experiments and have been used with MEB
(Dansereau et al., 2016), VP (Ringeisen et al., 2019) and discrete element (Herman, 2016) models to assess ice fracture characteristics, LKF angles and
intermittency. In this analysis, we use the experiment designed by Ringeisen et al. (2019) to test the sensitivity of the residual error growth, sea ice
deformation and LKF orientation on the correction path angle *γ* in the generalized stress correction scheme. The model domain is
250 km × 100 km with 1 km spatial resolution. The initial conditions are 1 m ice thickness and 100 %
concentration in the middle 60 km of the domain with two narrow bands of open water (20 km width) on each side
(Fig. 2). A solid-wall Dirichlet boundary condition ($u=v=\mathrm{0}$) is used at the bottom, and open-water Neumann boundary conditions
($\partial u/\partial n=\mathrm{0}$) are used on the top and sides. In all experiments, the forcing is specified by a downward surface stress
*τ*_{a} (see Eq. 2) over the entire domain. This differs from Ringeisen et al. (2019) and Dansereau et al. (2016) where the
upper boundary is represented by a moving wall acting as external forcing. The magnitude of *τ*_{a} is ramped up from 0 to
0.60 N m^{2} (corresponding to ∼ 20 m s^{−1} winds or ∼ 0.33 m s^{−1} surface currents) in a 2 h
period and then remains constant.

Note that all simulations are performed without including heterogeneity in order to clearly identify the model performance (both numerics and physics), unless specified otherwise. This allows us to quantify the growth of residual numerical errors in a problem with full symmetry and their impact on the simulated LKF orientation and post-fracture sea ice deformations.

## 4.2 Numerical approaches

The MEB model is implemented in the McGill Sea Ice Model Version 5 (McGill SIM5) using an Eulerian, second-order finite difference numerical scheme
(Tremblay and Mysak, 1997; Lemieux et al., 2014; Plante et al., 2020). The equations are discretized in space using an Arakawa C grid and in time using a semi-implicit backward
Euler scheme (Plante et al., 2020). A solution to the non-linear momentum and constitutive equations (Eqs. 1 and 5) is
found using a Picard solver. The Picard solver uses an outer loop in which the equations are linearized and solved at each iteration using a
preconditioned flexible general minimum residual method (FGMRES, Lemieux et al., 2008). The non-linear terms are then updated and the linear problem
solved again until the residual error *ϵ*_{res}, defined as the L2 norm of the solution residual vector, is lower than
10^{−8} N m^{2} (Lemieux et al., 2014, for details). The prognostic equations for the tracers (Eqs. 3, 4
and 14) are updated within the outer loop iteration using an implicit–explicit (IMEX) approach (Lemieux et al., 2014). The reader is
referred to Plante et al. (2020) for more details.

## 4.3 Diagnostics

### 4.3.1 Field asymmetry

We monitor the influence of the residual errors on the model solution in the simulations using a normalized domain-integrated asymmetry factor
(*ϵ*_{asym}) in the maximum shear stress invariant field (*σ*_{II}). This diagnostic measures the asymmetry in the model
solution about the *y* axis (the vertical center line) and represents a measure of the numerical accuracy given that the model equations, initial
conditions and boundary conditions are all fully symmetric. The asymmetry factor is defined as

where *i* and *j* are the *x*–*y* grid indices, respectively; *n*_{x} and *n*_{y} are the number of grid cells in the *x* and *y* directions; and *a* and *b* are the
indices of the first and last ice-covered grid cells on the *x* axis.

Note that the field asymmetry measures the degradation of the originally fully symmetric problem as numerical errors are integrated and includes the
physical response to the integrated errors. This is in contrast with the residual error amplification ratio *R*, which is a measure of the local
amplification of the residual error by the damage parameterization at a given time step. The maximum *R* values in the domain at each time step
(*R*_{max}) are also shown below to visualize the contribution of the damage parameterization to the growth of the residual errors.

### 4.3.2 Damage activity

We quantify the development of fractures in the experiments using the damage activity *D*, defined as the total damage integrated over the original
ice domain in a given time interval Δ (= 60 s):

This parameter is analogous to the damage rate in Dansereau et al. (2016, 2017) and is used to identify the time at which the ice fractures. Note that this definition of damage activity (or damage rate) emphasizes activity in undamaged ice (i.e. new fractures) and is not sensitive to activity in already heavily damaged ice.

### 4.3.3 Fracture angle

The angles between conjugate LKFs in the Arctic are often discussed in relation with the orientation of the smaller-scale brittle fractures observed in the laboratory under uniaxial compression loads (i.e., Marko and Thomson, 1977; Schulson, 2004). The orientation of such compressive-shear fractures is often related to brittle fracture theories (e.g. to the development of wing cracks; Schulson, 2004; Wachter et al., 2009) and in terms of granular properties such as Coulombic friction or dilatancy (Erlingsson, 1988; Tremblay and Mysak, 1997; Overland et al., 1998).

Here, we define the fracture angle *θ* as the angle between the *y* axis and the fracture lines (see Fig. 2), and we compare the
simulated fracture angles in our experiments to two theories that are often used to describe the orientation of fractures: the Mohr–Coulomb fracture
theory and the Roscoe theory of dilatancy. Widely used in geoscience and engineering, the Mohr–Coulomb theory (Coulomb, 1773; Mohr, 1900) relates the
orientation of fractures to the angle of internal friction, as

In the Roscoe theory (Roscoe, 1970), the fracture angle is defined instead in terms of the angle of dilatancy (*δ*) of the granular material:

If *δ*=*ϕ*, the two theories give the same fracture angle *θ*. In general, the fracture angle in geomaterial and soils falls between
values predicted by the Mohr–Coulomb and Roscoe theories with zero dilatancy (*δ*=0) (Arthur et al., 1977; Bardet, 1991).

In our experiment, the fracture angle is calculated graphically for each individual simulation. We define the uncertainty as
± $\mathrm{tan}(W/L)\sim $ ± 2^{∘}, where *W* is the fracture width (typically a few grid cells wide, or ∼ 2–5 km) and *L* is
the fracture length (∼ 45 km). This error increases to ± 6^{∘} for the few cases where the fracture is not as localized.

## 5.1 Control simulation: standard damage parameterization

In the control simulation, a pair of conjugate LKFs first appear when the surface forcing *τ*_{a} = 0.29 N m^{2}, along
with secondary lines that are the results of interactions between the ice floe and the solid boundary that extends across the full width of the domain
at the base (Fig. 3). All LKFs are oriented at 39^{∘} from the *y* axis, smaller than reported by Dansereau et al. (2019) using a
finite element implementation of the same model ($\mathit{\theta}=\sim $ 43^{∘}) and higher than seen in observations (*θ* = ∼ 15–25^{∘}; Marko and Thomson, 1977; Hibler III and Schulson, 2000; Schulson, 2004; Hutter et al., 2021). This orientation also falls in between that predicted
by the Mohr–Coulomb (*θ* = 22.5^{∘}) and Roscoe theories (*θ* = 4^{∘} when *δ*=0), in accord with the common
observation that both the angle of internal friction and the dilatancy (*δ*) are important in defining the fault orientation
(Arthur et al., 1977; Vardoulakis, 1980; Balendran and Nemat-Nasser, 1993).

The deformation along the fully developed LKFs in our experiment is mostly shear and convergent (i.e. ridging, Fig. 3c and d). This
contrasts with the early stage of the LKF development during which the material response to the new damage is elastic and shows mostly divergent
deformations (see the positive strain rates in Fig. 4b). This elastic response to damage influences the propagation of the
fractures in space at short timescales (seconds) governed by the elastic wave speed. The convergent deformations only develop over a longer
timescale as the sea ice deformation continues post-fracture in the damaged ice, and the deformation transitions from the elastic- to the
viscous-dominated regime. This transition is clearly seen in the development of a linear dependence between stress and strain rate invariants (scaled
by (1−*d*)^{3}), where the slope corresponds to the viscosity (see the transition from Fig. 4b and d–f). The simulation
reaches steady state with deformations that are fully viscous and localized in the heaviest damage areas (Fig. 4e
and f). This causes a predominance of shear and convergence deformation along the LKFs throughout the simulation.

The asymmetries in the solution are very small at the beginning of the simulation (*t*≤ 57 min) and do not grow until fractures occur
(Fig. 5a and b). As the LKFs develop, small errors grow rapidly, with *ϵ*_{asym} increasing in large steps crossing
multiple orders of magnitude. Note that the model is always iterated to convergence with a strict residual error tolerance
(*ϵ*_{res} = 10^{−8} N m^{2}). The steep growth in *ϵ*_{asym} is associated with large (>1) values of the
error amplification ratio *R* (see Eq. 27), which reach ∼ 20 in the control simulation (Fig. 5b). Since
*ϵ*_{asym} is a domain-integrated quantity, it increases in time following large local error growths *R*. This illustrates the
long-range and long-term influence of residual errors, which act on the development of the future fractures. Note that *ϵ*_{asym}
saturates when the *σ*_{II} field is no longer symmetric and becomes insensitive to additional error growth. We assess the precision of
the solution using the maximum error amplification ratio *R*_{max}, which indicates the level of amplification of residual errors in the simulations,
at times by more than 1 order of magnitude locally (*R*_{max} > 10).

## 5.2 Generalized stress correction

The generalized damage parameterization reduces the growth of residual errors, with decreasing asymmetry factor and maximum error amplification ratio
*R*_{max} for increasing path angle *γ* (Fig. 6). In particular, using *γ* > 0^{∘} stabilizes the
damage parameterization and eliminates the large spikes in *R*_{max} seen in the control simulation or when using *γ* < 0^{∘}, where the
amplification ratio *R* increases by up to 2 orders of magnitude locally (Fig. 6b). The increased stability results in an
overall smaller and smoother growth of the asymmetry factor *ϵ*_{asym} (Fig. 6a), allowing for longer-term
symmetrical simulations that include post-fracture deformations. Note that despite this improvement, the asymmetry factor *ϵ*_{asym}
still grows over time as the simulations remain sensitive to the residual errors in heavily damaged ice, due to the non-linear relationship between
the sea ice deformation and the damage. This effect is less important when using large correction path angles (*γ* > 45^{∘}) due to a
slower LKF development, as discussed below.

Results show that the LKF orientation is sensitive to the decohesive stress tensor, with a decreasing angle *θ* for increasing stress correction
path angle *γ* (Fig. 7). This finding is in line with results from Dansereau et al. (2019), where the orientation of faults was
related to the far-field stress associated with the collective damage. In the MEB model, the far-field stresses directly depend on the corrected
stress state, which includes *σ*_{D} in the generalized damage parameterization. Increasing the correction path angle *γ* reduces
the LKF angles, in better agreement with observations.

The correction path angle *γ* influences the time integration required to reach the same damage and deformation rates
(Fig. 8) along the LKFs. This is due to the fact that increasing the angle *γ* reduces the amount of damage for the same
super-critical stress state because the stress correction path approaches the horizontal and Ψ is closer to 1. The simulated ice deformations are
otherwise mostly insensitive to the correction path angle; i.e. all simulations have divergence during the initial elastic response when the ice
fractures are followed by a transition to viscous deformations where shear and convergence deformations are predominant (Fig. 8a). In
contrast with plastic flow (Ringeisen et al., 2019, 2021) or typical granular material behaviour (Balendran and Nemat-Nasser, 1993; Tremblay and Mysak, 1997), divergent
post-fracture deformation is only present when tensile stresses develop, e.g. at the intersection between conjugate LKFs. This behaviour stems from
the use of post-fracture viscosity to represent the large-scale sea ice deformations and differs from classical VP models, which represent the
observed plasticity of sea ice deformations at the macro-scale (Coon et al., 1974; Tuhkuri and Lensu, 2002) but do not represent the brittle component of the
fractures or discontinuities in material properties.

## 5.3 Sensitivity to *ϕ* and *ν*

Repeating the experiment using different angles of internal friction (*ϕ*) shows that the LKF orientations decrease with increasing *ϕ*. The
simulated angles *θ* fall within the envelope from the Mohr–Coulomb and Roscoe theories, except for small angles of internal friction
(*ϕ* < 20^{∘}), a value that is rarely observed for granular materials (Fig. 9). Note that the sensitivity of
the LKF orientation to the coefficient of internal friction also disappears for small angles of internal friction (*ϕ* < 20^{∘}) when
using a large correction path angle (*γ* = 60^{∘} in Fig. 7). When both the stress correction path and the yield
criterion approach horizontal, fracture yields large stress corrections but small damage increases (i.e., Ψ=1), such that the LKF orientation
is mostly governed by the stress correction and weakly sensitive to other model parameters. Based on these results, we suggest the use of a correction
path that is normal to the yield criterion (*γ*=arctan *μ*; see black points in Fig. 9).

Decreasing the angle of internal friction reduces the shear strength of sea ice for a given normal stress, such that the fracture develops earlier in the simulation (i.e. under smaller surface forcing, Fig. 10). It also reduces the divergence associated with the elastic response when ice fractures and increases the convergence in the post-fracture viscous regime. This result is typical for granular material, with smaller fault orientations (larger angles of internal friction) associated with larger angles of dilatancy (e.g. the sawtooth model of Bolton, 1986).

The orientation of LKFs is not sensitive to the Poisson ratio when the generalized stress correction scheme is used with a fixed stress correction
path angle *γ* (Fig. 11). This is in contrast with simulations using the standard stress correction scheme, where the
fracture angle decreases with increasing *ν* (see blue points in Fig. 11 and Dansereau et al., 2019). Note that the
Poisson ratio also affects the amount of shear and normal stress concentration associated with a local discontinuity in material properties
(Karimi and Barrat, 2018). The fact that the LKF orientation is not affected by the changes in Poisson ratio thus indicates that the stress concentration and
propagation of the fracture in space are mainly controlled by the stress correction rather than by the relaxation of material properties with
damage. We speculate that the sensitivity of the LKF orientation to the Poisson ratio in the standard stress correction scheme stems from the
dependency of the stress correction path angle to the super-critical stress state
(i.e. *γ* = ${\mathrm{tan}}^{-\mathrm{1}}({\mathit{\sigma}}_{\mathrm{I}}^{\prime}/{\mathit{\sigma}}_{\text{II}}^{\prime})$).

The results presented above show that the generalized stress correction scheme reduces the growth of the residual error associated with the damage
parameterization. Despite the improvement, some asymmetries are still present in the simulations (*ϵ*_{asym} < 10^{−2}). This is
due to the memory in the damage parameter (i.e. an integrated quantity) where residual errors accumulate and influence the temporal evolution of the
solution. In regions of heavily damaged ice, the integrated errors in the damage parameter result in large errors in the stress state due to the cubic
dependence of the Maxwell viscosity *η* on *d* (Eq. 9). Future work includes replacing this formulation with a function that
decreases the sensitivity of the Maxwell viscosity *η* for small changes in *d* around *d*=1.

Overall, the use of a decohesive stress tensor yields smaller simulated LKF angles, without significantly impacting the material deformations. Using a
large correction path angle *γ* (> 45^{∘}), however, significantly slows the damage production and reduces the simulated sensitivity of
the LKF orientation to the mechanical strength parameters. Based on these results, we suggest using a correction path that is normal to the yield
criterion (*γ*=arctan *μ*). This value brings the simulated LKF angles closer to observations (see black points in
Fig. 9) and reduces the amplification of residual errors, while correcting the super-critical stresses towards the closest
point on the yield curve. Our implementation thus represents a generalization of the damage parameterization that can be easily implemented
numerically and used to improve the performance of MEB models. Whether these improvements are also seen in the context of pan-Arctic simulations
remains to be tested and is the subject of future work.

The simulation results show that in the MEB model, the damage develops at short timescales during which the elastic component of the rheology is
important, while most of the deformations occur post-fracture over a longer timescale in the heavily damaged ice. This is in contrast with plastic
models, in which a flow rule simultaneously dictates both the LKF development and the relative amount of shear and normal deformations occurring along
the LKFs. The decoupling between the development of damage and the post-fracture deformations in the MEB model explains that the type of deformations
in the LKFs remains similar (uniaxial convergence, i.e. ridging, contrary to observation; Stern et al., 1995) despite the use of a different stress
correction path *γ*. This behaviour stems from the dominance of the viscous regime post-fracture: lead opening cannot occur when the stress state
is compressive and remains limited to locations where tensile stresses are present, such as at the intersection of the LKFs. This is contrary to
granular theories, in which the distribution of contact normals determines the amount of ridging or lead opening (i.e. dilatancy) that is occurring
when forced in uniaxial compression (Balendran and Nemat-Nasser, 1993). This indicates that the decohesive stress tensor cannot be used to influence the
deformations associated with the fracture of ice in the MEB rheology unless other parameterizations, such as including a decohesive strain tensor during
the fractures (e.g., see Schreyer et al., 2006; Sulsky and Peterson, 2011), are added to the rheology.

The viscous dissipation timescale (*λ*) in our model is set based on observations (∼ 10^{5}; Tabata, 1955; Hata and Tremblay, 2015) and is 1
order of magnitude smaller than in other MEB implementations (Dansereau et al., 2016; Rampal et al., 2019). The results from the model are robust with respect to
the exact value of *λ* for a range 10^{5}–10^{7}, with the increased *λ* being compensated for by larger damage values along the LKFs. For even
larger *λ* values, divergent deformations persist longer in the simulation, and the transition from an elastic- to viscous-dominated regime occurs
later in the simulation (see Fig. 12), decreasing the overall convergence along the LKFs. If the transition to the viscous regime is removed
(e.g. by setting *α*=1), divergence dominates throughout the simulations and reaches large values as the leads open. The elastic waves,
however, are no longer dissipated in the LKFs, leading to large and noisy deformation fields (divergence or convergence). These findings call for a
different viscosity dependence on damage, leading to both dissipation of elastic waves and a more realistic post-fracture deformation field.

Note that the results presented above were presented using a single space and time resolution and ice sample aspect ratio and without using heterogeneity. While the exact localization of the LKFs in the simulations is affected by these parameters, the overall physics and sensitivity to the damage parameterization are robust to these changes. For instance, repeating the experiment by doubling the space resolution or the width of the ice sample does not change the LKF position and orientation (not shown). On the other hand, adding heterogeneity changes the LKF development by forming irregular sliding planes instead of the linear diamond shapes (Fig. 13a), naturally creating contact points where ridging occurs with lead opening elsewhere along the LKFs. This effectively creates a form of dilatancy typical of granular materials (see alternating divergence and convergence in Fig. 13c) and leads to the formation of many secondary fractures, but the overall LKF orientations and their sensitivities otherwise remain the same as presented in this paper. Heterogeneity was also documented to be responsible for the localization and intermittency of the sea ice fractures, properties that are not investigated in our paper. These properties and their sensitivity to the decohesive stress tensor and other physical or numerical parameters require more investigation and are the subject of future work.

We propose a generalized stress correction scheme for the damage parameterization to reduce the growth of residual errors in the MEB sea ice model documented in Plante et al. (2020). To this end, we scale the damage factor Ψ based on the super-critical maximum shear stress invariant (${\mathit{\sigma}}_{\text{II}}^{\prime}$) only, together with a decohesive stress tensor defining the path from the super-critical stress state to the yield curve. With this added flexibility to the choice of stress correction path, we determine the influence of the super-critical stress correction on the simulated sea ice deformations and LKF orientation in the context of uniaxial compression experiments similar to those presented in Ringeisen et al. (2019). This knowledge will serve as a basis for the development of other components to the damage parameterization to improve the simulated sea ice deformations.

Our results show that in the MEB rheology, most of the deformations occur post-fracture in heavily damaged ice, where the viscous term is dominant. This causes a predominance of convergence (ridging) in the LKFs, contrary to laboratory experiments of granular materials and satellite observations of sea ice. The use of a decohesive stress tensor influences the LKF orientation in the sea ice cover but does not influence the type of deformation rates (convergence and shear) or the simulated dilatancy. Future work will involve the modification of the non-linear relationship between the Maxwell viscosity and the damage. We also show that the sensitivity of the LKF orientation to the Poisson ratio, seen when using the standard damage parameterization, disappears when using the generalized stress correction scheme with a fixed stress correction path. This suggests that in the MEB model the stress concentration and fracture propagation are governed by the stress correction rather than by the relaxation of the mechanical properties associated with the damage.

Based on our results, using the generalized damage parameterization with a stress correction path normal to the yield curve reduces the growth of
residual errors and allows longer-term simulations with post-fracture deformations. Using this stress correction path also reduces the orientation of
LKFs by ∼ 5^{∘}, bringing them closer to observations. Despite these improvements, some error growth remains inherent to the formulation of
the damage parameterization. Whether this might be improved by removing the dependency of the damage parameters on the damage factor (and on the
super-critical stress state) will be explored in future work.

Our sea ice model code is available upon request.

Our model outputs are available upon request.

MP coded the model, ran all the simulations, analyzed results and led the writing of the manuscript. BT participated in regular discussions during the course of the work and edited the manuscript.

The contact author has declared that neither they nor their co-author have any competing interests.

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

This work is a contribution to the research program of Québec-Océan and to the ArcTrain International Training Program. We thank the three anonymous reviewers for their useful comments and suggestions during the open discussion process. We also thank Amélie Bouchat, Damien Ringeisen, Martin Losch and Jean-François Lemieux for useful discussions during the implementation of the MEB model and the generalized stress correction.

We are grateful to the Fonds de recherche du Québec – Nature et technologies (FRQNT) for financial support to Mathieu Plante during the course of this work as well as to the Natural Science and Engineering and Research Council (NSERC) Discovery Program and the Environment and Climate Change Canada Grant & Contribution for grants awarded to Bruno Tremblay.

This paper was edited by Yevgeny Aksenov and reviewed by three anonymous referees.

Amitrano, D. and Helmstetter, A.: Brittle creep, damage and time to failure in rocks, J. Geophys. Res.-Sol. Ea., 111, B11201, https://doi.org/10.1029/2005JB004252, 2006. a, b

Amitrano, D., Grasso, J.-R., and Hantz, D.: From diffuse to localised damage through elastic interaction, Geophys. Res. Lett., 26, 2109–2112, 1999. a

Arthur, J. R. F., Dunstan, T., Al-Ani, Q. A. J. L., and Assadi, A.: Plastic deformation and failure in granular media, Géotechnique, 27, 53–74, https://doi.org/10.1680/geot.1977.27.1.53, 1977. a, b

Balendran, B. and Nemat-Nasser, S.: Double sliding model for cyclic deformation of granular materials, including dilatancy effects, J. Mech. Phys. Solids, 41, 573–612, https://doi.org/10.1016/0022-5096(93)90049-L, 1993. a, b, c

Bardet, J.: Orientation of shear bands in frictional soils, J. Eng. Mech.-ASCE, 117, 1466–1484, https://doi.org/10.1061/(ASCE)0733-9399(1991)117:7(1466), 1991. a, b

Bolton, M. D.: The strength and dilatancy of sands, Geotechnique, 36, 65–78, https://doi.org/10.1680/geot.1986.36.1.65, 1986. a

Bouchat, A. and Tremblay, B.: Using sea-ice deformation fields to constrain the mechanical strength parameters of geophysical sea ice, J. Geophys. Res.-Oceans, 122, 5802–5825, https://doi.org/10.1002/2017JC013020, 2017. a

Bouchat, A., Hutter, N. C., Chanut, J., Dupont, F., Dukhovskoy, D. S., Garric, G., Lee, Y. J., Lemieux, J.-F., Lique, C., Losch, M., Maslowski, W., Myers, P. G., Òlason, E. Ö., Rampal, P., Rasmussen, T. A. S., Talandier, C., Tremblay, B., and Wang, Q.: Sea Ice Rheology Experiment (SIREx), Part I: Scaling and statistical properties of sea-ice deformation fields, Earth and Space Science Open Archive, p. 36, https://doi.org/10.1002/essoar.10507397.1, 2021. a, b, c

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

Carrier, A., Got, J.-L., Peltier, A., Ferrazzini, V., Staudacher, T., Kowalski, P., and Boissier, P.: A damage model for volcanic edifices: Implications for edifice strength, magma pressure, and eruptive processes, J. Geophys. Res.-Sol. Ea., 120, 567–583, https://doi.org/10.1002/2014JB011485, 2015. a

Coon, M. D., Maykut, G. A., Pritchard, R. S., Rothrock, D. A., and Thorndike, A. S.: Modeling the pack ice as an elastic-plastic material, AIDJEX bulletin, 24, 1–105, https://doi.org/10.1017/CBO9781107415324.004, 1974. a, b

Coulomb, C.: Test on the applications of the rules of maxima and minima to some problems of statics related to architecture, Mem. Math. Phys., 7, 343–382, 1773. a

Damsgaard, A., Adcroft, A., and Sergienko, O.: Application of Discrete Element Methods to Approximate Sea Ice Dynamics, J. Adv. Model. Earth Sy., 10, 2228–2244, https://doi.org/10.1029/2018MS001299, 2018. a

Dansereau, V., Weiss, J., Saramito, P., and Lattes, P.: A Maxwell elasto-brittle rheology for sea ice modelling, The Cryosphere, 10, 1339–1359, https://doi.org/10.5194/tc-10-1339-2016, 2016. a, b, c, d, e, f, g, h, i

Dansereau, V., Weiss, J., Saramito, P., Lattes, P., and Coche, E.: Ice bridges and ridges in the Maxwell-EB sea ice rheology, The Cryosphere, 11, 2033–2058, https://doi.org/10.5194/tc-11-2033-2017, 2017. a, b

Dansereau, V., Démery, V., Berthier, E., Weiss, J., and Ponson, L.: Collective Damage Growth Controls Fault Orientation in Quasibrittle Compressive Failure, Phys. Rev. Lett., 122, 085501, https://doi.org/10.1103/PhysRevLett.122.085501, 2019. a, b, c, d, e, f

Erlingsson, B.: Two-Dimensional Deformation Patterns in Sea Ice, J. Glaciol., 34, 301–308, https://doi.org/10.3189/S0022143000007061, 1988. a

Girard, L., Bouillon, S., Weiss, J., Amitrano, D., Fichefet, T., and Legat, V.: A new modeling framework for sea-ice mechanics based on elasto-brittle rheology, Ann. Glaciol., 52, 123–132, https://doi.org/10.3189/172756411795931499, 2011. a

Hata, Y. and Tremblay, L. B.: Anisotropic internal thermal stress in sea ice from the Canadian Arctic Archipelago, J. Geophys. Res.-Oceans, 120, 5457–5472, https://doi.org/10.1002/2015JC010819, 2015. a

Herman, A.: Discrete-Element bonded-particle Sea Ice model DESIgn, version 1.3a – model description and implementation, Geosci. Model Dev., 9, 1219–1241, https://doi.org/10.5194/gmd-9-1219-2016, 2016. a

Hibler, W. D.: A dynamic thermodynamic sea ice model, J. Phys. Oceanogr., 9, 815–846, 1979. a, b, c

Hibler III, W. D. and Schulson, E. M.: On modeling the anisotropic failure and flow of flawed sea ice, J. Geophys. Res.-Oceans, 105, 17105–17120, https://doi.org/10.1029/2000JC900045, 2000. a

Hunke, E. C. and Dukowicz, J. K.: An Elastic–Viscous–Plastic Model for Sea Ice Dynamics, J. Phys. Oceanogr., 27, 1849–1867, https://doi.org/10.1175/1520-0485(1997)027<1849:AEVPMF>2.0.CO;2, 1997. a, b

Hutter, N., Losch, M., and Menemenlis, D.: Scaling Properties of Arctic Sea Ice Deformation in a High-Resolution Viscous-Plastic Sea Ice Model and in Satellite Observations, J. Geophys. Res.-Oceans, 123, 672–687, https://doi.org/10.1002/2017JC013119, 2018. a

Hutter, N. C., Bouchat, A., Dupont, F., Dukhovskoy, D. S., Koldunov, N. V., Lee, Y. J., Lemieux, J.-F., Lique, C., Losch, M., Maslowski, W., Myers, P. G., Òlason, E. Ö., Rampal, P., Rasmussen, T. A. S., Talandier, C., Tremblay, B., and Wang, Q.: Sea Ice Rheology Experiment (SIREx), Part II: Evaluating simulated linear kinematic features in high-resolution sea-ice simulations, Earth and Space Science Open Archive, p. 35, https://doi.org/10.1002/essoar.10507396.1, 2021. a, b, c, d

Itkin, P., Spreen, G., Hvidegaard, S. M., Skourup, H., Wilkinson, J., Gerland, S., and Granskog, M. A.: Contribution of Deformation to Sea Ice Mass Balance: A Case Study From an N-ICE2015 Storm, Geophys. Res. Lett., 45, 789–796, https://doi.org/10.1002/2017GL076056, 2018. a

Karimi, K. and Barrat, J.-L.: Correlation and shear bands in a plastically deformed granular medium, Sci. Rep., 8, 4021, https://doi.org/10.1038/s41598-018-22310-z, 2018. a

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

Koldunov, N. V., Danilov, S., Sidorenko, D., Hutter, N., Losch, M., Goessling, H., Rakowsky, N., Scholz, P., Sein, D., Wang, Q., and Jung, T.: Fast EVP Solutions in a High-Resolution Sea Ice Model, J. Adv. Model. Earth Sy., 11, 1269–1284, https://doi.org/10.1029/2018MS001485, 2019. a

Kozo, T. L.: Initial model results for Arctic mixed layer circulation under a refreezing lead, J. Geophys. Res.-Oceans, 88, 2926–2934, https://doi.org/10.1029/JC088iC05p02926, 1983. a

Ledley, T. S.: A coupled energy balance climate-sea ice model: Impact of sea ice and leads on climate, J. Geophys. Res.-Atmos., 93, 15919–15932, https://doi.org/10.1029/JD093iD12p15919, 1988. a

Lemieux, J.-F., Tremblay, B., Thomas, S., Sedláček, J., and Mysak, L. A.: Using the preconditioned Generalized Minimum RESidual (GMRES) method to solve the sea-ice momentum equation, J. Geophys. Res.-Oceans, 113, C10004, https://doi.org/10.1029/2007JC004680, 2008. a, b, c

Lemieux, J.-F., Knoll, D. A., Losch, M., and Girard, C.: A second-order accurate in time IMplicit–EXplicit (IMEX) integration scheme for sea ice dynamics, J. Comput. Phys., 263, 375–392, https://doi.org/10.1016/j.jcp.2014.01.010, 2014. a, b, c, d

Li, X., Krueger, S. K., Strong, C., Mace, G. G., and Benson, S.: Midwinter Arctic leads form and dissipate low clouds, Nat. Commun., 11, 206, https://doi.org/10.1038/s41467-019-14074-5, 2020. a

Lüpkes, C., Vihma, T., Birnbaum, G., and Wacker, U.: Influence of leads in sea ice on the temperature of the atmospheric boundary layer during polar night, Geophys. Res. Lett., 35, L03805, https://doi.org/10.1029/2007GL032461, 2008. a

Marko, J. R. and Thomson, R. E.: Rectilinear leads and internal motions in the ice pack of the western Arctic Ocean, J. Geophys. Res., 82, 979–987, https://doi.org/10.1029/JC082i006p00979, 1977. a, b

Matsumura, Y. and Hasumi, H.: Brine-Driven Eddies under Sea Ice Leads and Their Impact on the Arctic Ocean Mixed Layer, J. Phys. Oceanogr., 38, 146–163, https://doi.org/10.1175/2007JPO3620.1, 2008. a

Maykut, G. A.: Large-scale heat exchange and ice production in the central Arctic, J. Geophys. Res.-Oceans, 87, 7971–7984, https://doi.org/10.1029/JC087iC10p07971, 1982. a

Mohr, O.: Welche Umstände bedingen die Elastizitätsgrenze und den Bruch eines Materials, Z. Ver. Dtsch. Ing., 46, 1572–1577, 1900. a

Murakami, S.: Continuum damage mechanics: a continuum mechanics approach to the analysis of damage and fracture, https://doi.org/10.1007/978-94-007-2666-6, 2012. a, b

Overland, J. E., McNutt, S. L., Salo, S., Groves, J., and Li, S.: Arctic sea ice as a granular plastic, J. Geophys. Res.-Oceans, 103, 21845–21867, https://doi.org/10.1029/98JC01263, 1998. a

Plante, M., Tremblay, B., Losch, M., and Lemieux, J.-F.: Landfast sea ice material properties derived from ice bridge simulations using the Maxwell elasto-brittle rheology, The Cryosphere, 14, 2137–2157, https://doi.org/10.5194/tc-14-2137-2020, 2020. a, b, c, d, e, f, g, h, i, j, k, l, m, n

Rampal, P., Bouillon, S., Ólason, E., and Morlighem, M.: neXtSIM: a new Lagrangian sea ice model, The Cryosphere, 10, 1055–1073, https://doi.org/10.5194/tc-10-1055-2016, 2016. a, b, c

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

Rice, J. R.: Solid Mechanics, Harvard University, Cambridge, MA, 2010. a

Ringeisen, D., Losch, M., Tremblay, L. B., and Hutter, N.: Simulating intersection angles between conjugate faults in sea ice with different viscous–plastic rheologies, The Cryosphere, 13, 1167–1186, https://doi.org/10.5194/tc-13-1167-2019, 2019. a, b, c, d, e, f

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

Roscoe, K. H.: The Influence of Strains in Soil Mechanics, Géotechnique, 20, 129–170, https://doi.org/10.1680/geot.1970.20.2.129, 1970. a

Schreyer, H. L., Sulsky, D. L., Munday, L. B., Coon, M. D., and Kwok, R.: Elastic-decohesive constitutive model for sea ice, J. Geophys. Res.-Oceans, 111, C11S26, https://doi.org/10.1029/2005JC003334, 2006. a, b, c, d, e

Schulson, E. M.: Compressive shear faults within arctic sea ice: Fracture on scales large and small, J. Geophys. Res.-Oceans, 109, C07016, https://doi.org/10.1029/2003JC002108, 2004. a, b, c

Sodhi, D. S.: Ice arching and the drift of pack ice through restricted channels, Crrel report 77-18, Cold Regions Research And Engineering Laboratory, Hanover, New Hampshire, 1977. a

Stern, H. L., Rothrock, D. A., and Kwok, R.: Open water production in Arctic sea ice: Satellite measurements and model parameterizations, J. Geophys. Res.-Oceans, 100, 20601–20612, https://doi.org/10.1029/95JC02306, 1995. a

Sulsky, D. and Peterson, K.: Toward a new elastic–decohesive model of Arctic sea ice, Physica D, special Issue: Fluid Dynamics: From Theory to Experiment, 240, 1674–1683, https://doi.org/10.1016/j.physd.2011.07.005, 2011. a, b, c

Tabata, T.: A Measurement of Visco-Elastic Constants of Sea Ice, Journal of the Oceanographical Society of Japan, 11, 185–189, 1955. a

Timco, G. and Weeks, W.: A review of the engineering properties of sea ice, Cold Reg. Sci. Technol., 60, 107–129, https://doi.org/10.1016/j.coldregions.2009.10.003, 2010. a

Tremblay, L.-B. and Hakakian, M.: Estimating the Sea Ice Compressive Strength from Satellite-Derived Sea Ice Drift and NCEP Reanalysis Data, J. Phys. Oceanogr., 36, 2165–2172, https://doi.org/10.1175/JPO2954.1, 2006. a

Tremblay, L.-B. and Mysak, L. A.: Modeling Sea Ice as a Granular Material, Including the Dilatancy Effect, J. Phys. Oceanogr., 27, 2342–2360, https://doi.org/10.1175/1520-0485(1997)027<2342:MSIAAG>2.0.CO;2, 1997. a, b, c, d, e

Tuhkuri, J. and Lensu, M.: Laboratory tests on ridging and rafting of ice sheets, J. Geophys. Res.-Oceans, 107, 8-1–8-14, https://doi.org/10.1029/2001JC000848, 2002. a, b

Turnbull, I. D., Torbati, R. Z., and Taylor, R. S.: Relative influences of the metocean forcings on the drifting ice pack and estimation of internal ice stress gradients in the Labrador Sea, J. Geophys. Res.-Oceans, 122, 5970–5997, https://doi.org/10.1002/2017JC012805, 2017. a

Vardoulakis, I.: Shear band inclination and shear modulus of sand in biaxial tests, Int. J. Numer. Anal. Met., 4, 103–119, https://doi.org/10.1002/nag.1610040202, 1980. a

Wachter, L., Renshaw, C., and Schulson, E.: Transition in brittle failure mode in ice under low confinement, Acta Mater., 57, 345–355, https://doi.org/10.1016/j.actamat.2008.09.021, 2009. a

Wang, K.: Observing the yield curve of compacted pack ice, J. Geophys. Res.-Oceans, 112, C05015, https://doi.org/10.1029/2006JC003610, 2007. a

Wilchinsky, A. V. and Feltham, D. L.: A continuum anisotropic model of sea-ice dynamics, P. Roy. Soc. Lond. A Math., 460, 2105–2140, https://doi.org/10.1098/rspa.2004.1282, 2004. a, b

Wilchinsky, A. V., Heorton, H. D. B. S., Feltham, D. L., and Holland, P. R.: Study of the Impact of Ice Formation in Leads upon the Sea Ice Pack Mass Balance Using a New Frazil and Grease Ice Parameterization, J. Phys. Oceanogr., 45, 2025–2047, https://doi.org/10.1175/JPO-D-14-0184.1, 2015. a