Articles | Volume 15, issue 12
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

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

1 Introduction

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 (Maykut1982; Ledley1988; 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 (Kozo1983; Matsumura and Hasumi2008). 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 (Hibler1979; Hunke and Dukowicz1997), 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 Mysak1997; Wilchinsky and Feltham2004; Schreyer et al.2006; Sulsky and Peterson2011; 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 (Hibler1979) and elastic–viscous–plastic (Hunke and Dukowicz1997) sea ice rheologies as well as the elastic–anisotropic (Wilchinsky and Feltham2004) 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 (Bardet1991; Wang2007), 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 Tremblay2017; 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 Helmstetter2006) and adapted for the large-scale modelling of sea ice (Girard et al.2011; Bouillon and Rampal2015; 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 Helmstetter2006; Carrier et al.2015). In continuum damage mechanics, the damage parameter is derived instead from thermodynamic laws (Murakami2012) 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 Peterson2011).

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.

Table 1Default model parameters.

Download Print Version | Download XLSX

2 Model

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 Mysak1997; 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),

(1) ρ i h u t = σ + τ ,

where ρi is the ice density, h is the mean ice thickness, u (=ui^+vj^) 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:

(2) τ τ a - ρ w C dw | u | u ,

where ρw is the water density, Cdw 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)

(5) σ t + 1 λ σ = E C : ϵ ˙ ,

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 ϵ˙ is the (second-order) strain rate tensor. The tensors C and ϵ˙ 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  Rice2010) 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 (0<d<1), a (= 20) is the standard ice concentration parameter (Hibler1979; Rampal et al.2016), λ0 (= 105s,  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 Lensu2002).

2.3 Yield criterion

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

(10) F ( σ ) = σ II + μ σ I - c < 0 ,



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

(13) c = c 0 h e - a ( 1 - A ) ,

where c0 (= 10 kN m−2) is the cohesion of sea ice derived from observations (Sodhi1977; Tremblay and Hakakian2006; Plante et al.2020) or laboratory experiments (Timco and Weeks2010). 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 Td (= 1 s) as

(14) d t = ( 1 - Ψ ) ( 1 - d ) T d ,


(15) Ψ = σ c σ = min 1 , c σ II + μ σ I ,

is a damage factor (0<Ψ<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):

(16) σ t = - ( 1 - Ψ ) σ T d .

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

Figure 1(a) Mohr–Coulomb yield criterion (±σII+μσI=c, blue lines) in stress invariant space. σ is the uncorrected super-critical stress state, σc is the critical stress state for a given correction path angle γ (red dashed line) and c is the cohesion. The decohesive stress tensor σD is defined as the difference between σc and the scaled super-critical stress (Ψσ). (b) Proposed correction paths for various super-critical stresses σ that minimize the error amplification ratio (R), which consist of the standard parameterization for large tensile stresses (orange) and a correction path with γ= 45 for small tensile and compressive stresses (purple). The green line indicates the transition between the two formulations.


3 Generalized 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 Murakami2012).

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

(17) Ψ = σ IIc σ II ,

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

(18) σ II = ( 1 / tan ( γ ) ) σ I + B ,

where B (=σII-1/tan(γ)σI) 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),

(19) σ IIc = c + μ tan ( γ ) σ II - μ σ I 1 + μ tan ( γ ) .

The damage factor can then be written in terms of the super-critical stress state invariants (σI, σII), the correction path angle γ and the coefficient of internal friction μ as

(20) Ψ = c + μ tan ( γ ) σ II - μ σ I ( 1 + μ tan ( γ ) ) σ II .

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 ΨσI. 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

(21) σ t = - ( 1 - Ψ ) σ + σ D T d ,

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

(22)σID=σIc-ΨσI=c-Ψ(σII-μσI)μ,(23)σIID=0,(by definition).

When tanγ=σI/σII and σID=σIID=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 Ψ(σI,σII) can be written as (Plante et al.2020)

(24) δ Ψ = Ψ σ I 2 δ σ I 2 + Ψ σ II 2 δ σ II 2 ,

where δσI and δσII are the errors of the calculated stress invariants. Using Eq. (21) and re-writing δσI and δσII in terms of the relative error ϵ (i.e., δσI=ϵσI, δσII=ϵσII), we obtain


where R is the error amplification ratio.

Given that the uncorrected stress is close to the yield criterion (i.e. σII+μσI-c0), the error amplification ratio R tends to infinity for

(28) tan ( γ ) = - 1 / μ ,

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 arctan(-1/μ)<θ< 90.

Note that the error amplification ratio R is small for σI<0 but becomes infinitely large at the yield curve tip when σII 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 σI/σII=tanγ, green line in Fig. 1b). The damage factor is then defined as

(29) Ψ = c + μ γ σ II - μ σ I ( 1 + μ γ ) σ II , if  σ I < σ II tan γ , c σ II + μ σ I , otherwise .

Figure 2Idealized domain for uniaxial compression simulations, with a solid boundary (Dirichlet conditions, u=v=0) at the bottom and open boundaries (Neumann conditions, u/n=0) on the sides and top. The initial conditions are h= 1 m and A= 100 % in a region of 250 km× 60 km in the center of the domain (white), with two 20 km wide bands of open water on each side (blue). The orientation of the LKFs (θ) is defined as half of the angle between conjugate pairs of fracture lines (orange lines).


4 Methods

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 (Herman2016) 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=0) is used at the bottom, and open-water Neumann boundary conditions (u/n=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 m2 (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 Mysak1997; 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−8N m2 (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

(30) ϵ asym = i = a b j = 1 n y | ( σ II ) i , j - ( σ II ) n x - i , j | i = a b j = 1 n y | ( σ II ) i , j ,

where i and j are the xy grid indices, respectively; nx and ny 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 (Rmax) 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):

(31) D = i = a b j = 1 n y d i , j t + Δ / 2 - d i , j t - Δ / 2 Δ .

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 Thomson1977; Schulson2004). The orientation of such compressive-shear fractures is often related to brittle fracture theories (e.g. to the development of wing cracks;  Schulson2004; Wachter et al.2009) and in terms of granular properties such as Coulombic friction or dilatancy (Erlingsson1988; Tremblay and Mysak1997; 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 (Coulomb1773; Mohr1900) relates the orientation of fractures to the angle of internal friction, as

(32) θ = π 4 - ϕ 2 .

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

(33) θ = π 4 - δ 2 .

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; Bardet1991).

In our experiment, the fracture angle is calculated graphically for each individual simulation. We define the uncertainty as ±tan(W/L)± 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.

Figure 3(a) Damage (unitless), (b) ice thickness (m, colour) and velocity vectors (m s−1), (c) mean normal strain rate invariant (ϵ˙I, d−1), and (d) maximum shear strain rate invariant (ϵ˙II, d−1), after 2 h of integration in the control simulation using the standard stress correction scheme.


5 Results

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 m2, 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 (θ= 43) and higher than seen in observations (θ= 15–25Marko and Thomson1977; Hibler III and Schulson2000; Schulson2004; 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; Vardoulakis1980; Balendran and Nemat-Nasser1993).

Figure 4Scatter plots of local stress invariants (σI vs. σII, in kN m−1, left column) and of the normal stresses and scaled strain rate invariants (σI vs. (1-d)3ϵ˙II, right column) in heavily damaged (d>0.9) grid cells, at t=57min (during the fracture development, top row), t=60min (a few minutes after the fracture, middle row) and t=90min ( 30 min after the fracture, bottom row). Colour indicates the local damage. The strain rates are normalized to account for the non-linear dependency of the viscosity η on the damage parameter. The gradual alignment of the points in the σI vs. (1-d)3ϵ˙II diagram indicates the development of a linear–viscous stress–strain relationship over time.


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.

Figure 5(a) Temporal evolution of the damage activity D; (b) the solution residual ϵres, asymmetry factor ϵasym and convergence criterion on ϵres; and (c) the maximum error amplification ratio Rmax, in the control simulation using the standard stress correction scheme.


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−8N m2). 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 Rmax, which indicates the level of amplification of residual errors in the simulations, at times by more than 1 order of magnitude locally (Rmax> 10).

Figure 6(a) Time evolution of the asymmetry factor ϵasym and (b) time series of the maximum error amplification ratio Rmax, in a sensitivity experiment on the stress correction path angle γ, using the generalized stress correction scheme.


5.2 Generalized stress correction

The generalized damage parameterization reduces the growth of residual errors, with decreasing asymmetry factor and maximum error amplification ratio Rmax for increasing path angle γ (Fig. 6). In particular, using γ> 0 stabilizes the damage parameterization and eliminates the large spikes in Rmax 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.

Figure 7Sensitivity of the LKF orientation θ on the stress correction path angle γ (degrees) in uniaxial loading experiments using the generalized stress correction schemes. The theoretical LKF angles from the Mohr–Coulomb and Roscoe theories are indicated by dashed–dotted and dashed lines, respectively, for reference.


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.

Figure 8Time evolution of the mean normal (a) and maximum shear (b) strain rate invariants integrated over the ice cover, in simulations using the generalized damage parameterization with a different stress correction path γ.


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-Nasser1993; Tremblay and Mysak1997), 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 Lensu2002) but do not represent the brittle component of the fractures or discontinuities in material properties.

Figure 9Sensitivity of the LKF orientation (θ, degrees) on the angle of internal friction (ϕ, degrees) in uniaxial loading experiments using different correction path angles (γ). The correction path angle γ=atan (μ) implies that the stress correction path is perpendicular to the yield curve. The theoretical LKF orientation from the Mohr–Coulomb and Roscoe theories is indicated by dashed–dotted and dashed lines, respectively, for reference.


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

Figure 10Time evolution of (a) the mean normal strain rate invariant integrated over the ice cover (d−1) and (b) the maximum shear strain rate invariant integrated over the ice cover (d−1), when using different angles of internal friction ϕ, with a stress correction path normal to the yield curve (γ=arctan (μ)).


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  Bolton1986).

Figure 11Sensitivity of the LKF orientation (θ, degrees) to the Poisson ratio (ν, unitless), in uniaxial loading experiments using different correction path angles (γ). The theoretical orientations from the Mohr–Coulomb and Roscoe theories are indicated by dashed–dotted and dashed lines, respectively, for reference.


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 Barrat2018). 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. γ=tan-1(σI/σII)).

Figure 12Time evolution of the mean normal strain rate invariant integrated over the ice cover (d−1) using a stress correction path normal to the yield curve (γ=arctan (μ)) with α=3 (blue), α=1 and a longer viscous dissipation timescale (λ=108 s).


6 Discussion

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-Nasser1993). 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 Peterson2011), are added to the rheology.

The viscous dissipation timescale (λ) in our model is set based on observations ( 105;  Tabata1955; Hata and Tremblay2015) 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 105107, 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.

Figure 13(a) Damage (unitless), (b) ice thickness (m, colour) and velocity vectors (m s−1), (c) mean normal strain rate invariant (ϵ˙I, d−1), and (d) maximum shear strain rate invariant (ϵ˙II, d−1) after 2 h of integration in using the generalized stress correction scheme with γ= 45 and including heterogeneity in the initial material cohesion field. The heterogeneous cohesion (c0) field is defined locally at each grid cell by picking a random number between 7.0 and 13.0 kN m−2. The remaining initial conditions are the same as all other simulations.


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.

7 Conclusion

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 (σII) 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.

Code availability

Our sea ice model code is available upon request.

Data availability

Our model outputs are available upon request.

Author contributions

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.

Competing interests

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.

Financial support

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.

Review statement

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,, 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,, 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,, 1993. a, b, c

Bardet, J.: Orientation of shear bands in frictional soils, J. Eng. Mech.-ASCE, 117, 1466–1484,, 1991. a, b

Bolton, M. D.: The strength and dilatancy of sands, Geotechnique, 36, 65–78,, 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,, 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,, 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,, 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,, 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,, 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,, 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,, 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,, 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,, 2019. a, b, c, d, e, f

Erlingsson, B.: Two-Dimensional Deformation Patterns in Sea Ice, J. Glaciol., 34, 301–308,, 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,, 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,, 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,, 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,, 2000. a

Hunke, E. C. and Dukowicz, J. K.: An Elastic–Viscous–Plastic Model for Sea Ice Dynamics, J. Phys. Oceanogr., 27, 1849–1867,<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,, 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,, 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,, 2018. a

Karimi, K. and Barrat, J.-L.: Correlation and shear bands in a plastically deformed granular medium, Sci. Rep., 8, 4021,, 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,, 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,, 2019. a

Kozo, T. L.: Initial model results for Arctic mixed layer circulation under a refreezing lead, J. Geophys. Res.-Oceans, 88, 2926–2934,, 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,, 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,, 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,, 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,, 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,, 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,, 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,, 2008. a

Maykut, G. A.: Large-scale heat exchange and ice production in the central Arctic, J. Geophys. Res.-Oceans, 87, 7971–7984,, 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,, 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,, 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,, 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,, 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,, 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,, 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,, 2021. a, b

Roscoe, K. H.: The Influence of Strains in Soil Mechanics, Géotechnique, 20, 129–170,, 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,, 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,, 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,, 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,, 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,, 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,, 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,<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,, 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,, 2017.  a

Vardoulakis, I.: Shear band inclination and shear modulus of sand in biaxial tests, Int. J. Numer. Anal. Met., 4, 103–119,, 1980. a

Wachter, L., Renshaw, C., and Schulson, E.: Transition in brittle failure mode in ice under low confinement, Acta Mater., 57, 345–355,, 2009. a

Wang, K.: Observing the yield curve of compacted pack ice, J. Geophys. Res.-Oceans, 112, C05015,, 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,, 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,, 2015. a

Short summary
We propose a generalized form for the damage parameterization such that super-critical stresses can return to the yield with different final sub-critical stress states. In uniaxial compression simulations, the generalization improves the orientation of sea ice fractures and reduces the growth of numerical errors. Shear and convergence deformations however remain predominant along the fractures, contrary to observations, and this calls for modification of the post-fracture viscosity formulation.