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

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

Table 1 Default model parameters.

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

$$\begin{array}{}\text{(1)}& {\mathit{\rho }}_{i}h{\displaystyle \frac{\partial \mathit{u}}{\partial t}}=\mathrm{\nabla }\cdot \mathit{\sigma }+\mathit{\tau },\end{array}$$

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:

$$\begin{array}{}\text{(2)}& {\displaystyle }\mathit{\tau }\approx {\mathit{\tau }}_{\mathrm{a}}-{\mathit{\rho }}_{\mathrm{w}}{C}_{\text{dw}}\left|\mathit{u}\right|\mathit{u},\end{array}$$

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

$$\begin{array}{}\text{(3)}& {\displaystyle }& {\displaystyle \frac{\partial h}{\partial t}}+\mathrm{\nabla }\cdot \left(h\mathit{u}\right)=\mathrm{0},\text{(4)}& {\displaystyle }& {\displaystyle \frac{\partial A}{\partial t}}+\mathrm{\nabla }\cdot \left(A\mathit{u}\right)=\mathrm{0},\end{array}$$

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)

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\partial \mathit{\sigma }}{\partial t}}+{\displaystyle \frac{\mathrm{1}}{\mathit{\lambda }}}\mathit{\sigma }=E\mathbf{C}:\dot{\mathit{ϵ}},\end{array}$$

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{ϵ}}$ is the (second-order) strain rate tensor. The tensors C and $\dot{\mathit{ϵ}}$ 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

$$\begin{array}{}\text{(6)}& {\displaystyle }& {\displaystyle }\mathbf{C}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}-{\mathit{\nu }}^{\mathrm{2}}}}\left(\begin{array}{ccc}\mathrm{1}& \mathit{\nu }& \mathrm{0}\\ \mathit{\nu }& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}-\mathit{\nu }\end{array}\right),\text{(7)}& {\displaystyle }& {\displaystyle }\left(\begin{array}{c}{\dot{\mathit{ϵ}}}_{xx}\\ {\dot{\mathit{ϵ}}}_{yy}\\ {\dot{\mathit{ϵ}}}_{xy}\end{array}\right)=\left(\begin{array}{c}\frac{\partial u}{\partial x}\\ \frac{\partial v}{\partial y}\\ \frac{\mathrm{1}}{\mathrm{2}}(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x})\end{array}\right),\end{array}$$

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

$$\begin{array}{}\text{(8)}& {\displaystyle }& {\displaystyle }E=Yh{e}^{-a(\mathrm{1}-A)}(\mathrm{1}-d),\text{(9)}& {\displaystyle }& {\displaystyle }\mathit{\lambda }={\mathit{\lambda }}_{\mathrm{0}}(\mathrm{1}-d{)}^{\mathit{\alpha }-\mathrm{1}},\end{array}$$

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), λ₀ (= 10⁵ 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,

$$\begin{array}{}\text{(10)}& F\left(\mathit{\sigma }\right)={\mathit{\sigma }}_{\text{II}}+\mathit{\mu }{\mathit{\sigma }}_{\mathrm{I}}-c<\mathrm{0},\end{array}$$

where

$$\begin{array}{}\text{(11)}& {\displaystyle }& {\displaystyle }{\mathit{\sigma }}_{\mathrm{I}}={\displaystyle \frac{{\mathit{\sigma }}_{xx}+{\mathit{\sigma }}_{yy}}{\mathrm{2}}},\text{(12)}& {\displaystyle }& {\displaystyle }{\mathit{\sigma }}_{\text{II}}=\sqrt{{\left({\displaystyle \frac{{\mathit{\sigma }}_{xx}-{\mathit{\sigma }}_{yy}}{\mathrm{2}}}\right)}^{\mathrm{2}}+{\mathit{\sigma }}_{xy}^{\mathrm{2}}},\end{array}$$

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

$$\begin{array}{}\text{(13)}& c=c_{\mathrm{0}}h{e}^{-a(\mathrm{1}-A)},\end{array}$$

where c₀ (= 10 kN m⁻²) 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

$$\begin{array}{}\text{(14)}& {\displaystyle \frac{\partial d}{\partial t}}={\displaystyle \frac{(\mathrm{1}-\mathrm{\Psi })(\mathrm{1}-d)}{T_{\mathrm{d}}}},\end{array}$$

where

$$\begin{array}{}\text{(15)}& \mathrm{\Psi }={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{c}}}{{\mathit{\sigma }}^{\prime }}}=min\left(\mathrm{1},{\displaystyle \frac{c}{{\mathit{\sigma }}_{\text{II}}^{\prime }+\mathit{\mu }{\mathit{\sigma }}_{\mathrm{I}}^{\prime }}}\right),\end{array}$$

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

$$\begin{array}{}\text{(16)}& {\displaystyle \frac{\partial {\mathit{\sigma }}^{\mathbf{\prime }}}{\partial t}}=-{\displaystyle \frac{(\mathrm{1}-\mathrm{\Psi }){\mathit{\sigma }}^{\mathbf{\prime }}}{T_{\mathrm{d}}}}.\end{array}$$

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 ($\pm {\mathit{\sigma }}_{\text{II}}+\mathit{\mu }{\mathit{\sigma }}_{\mathrm{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 Murakami, 2012).

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

$$\begin{array}{}\text{(17)}& \mathrm{\Psi }={\displaystyle \frac{{\mathit{\sigma }}_{\text{IIc}}}{{\mathit{\sigma }}_{\text{II}}^{\prime }}},\end{array}$$

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

$$\begin{array}{}\text{(18)}& {\mathit{\sigma }}_{\text{II}}=(\mathrm{1}/\tan (\mathit{\gamma }\left)\right){\mathit{\sigma }}_{\mathrm{I}}+B,\end{array}$$

where B (= ${\mathit{\sigma }}_{\text{II}}^{\prime }-\mathrm{1}/\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),

$$\begin{array}{}\text{(19)}& {\mathit{\sigma }}_{\text{IIc}}={\displaystyle \frac{c+\mathit{\mu }\tan \left(\mathit{\gamma }\right){\mathit{\sigma }}_{\text{II}}^{\prime }-\mathit{\mu }{\mathit{\sigma }}_{\mathrm{I}}^{\prime }}{\mathrm{1}+\mathit{\mu }\tan \left(\mathit{\gamma }\right)}}.\end{array}$$

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

$$\begin{array}{}\text{(20)}& \mathrm{\Psi }={\displaystyle \frac{c+\mathit{\mu }\tan \left(\mathit{\gamma }\right){\mathit{\sigma }}_{\text{II}}^{\prime }-\mathit{\mu }{\mathit{\sigma }}_{\mathrm{I}}^{\prime }}{(\mathrm{1}+\mathit{\mu }\tan (\mathit{\gamma }\left)\right){\mathit{\sigma }}_{\text{II}}^{\prime }}}.\end{array}$$

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⁻¹), 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⁻¹).

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

$$\begin{array}{}\text{(21)}& {\displaystyle \frac{\partial {\mathit{\sigma }}^{\prime }}{\partial t}}=-{\displaystyle \frac{(\mathrm{1}-\mathrm{\Psi }){\mathit{\sigma }}^{\prime }+{\mathit{\sigma }}_{\mathrm{D}}}{T_{\mathrm{d}}}},\end{array}$$

and the invariants of the decohesive stress tensor (σ_ID and σ_IID) are now defined as

$$\begin{array}{}\text{(22)}& {\displaystyle }& {\displaystyle }{\mathit{\sigma }}_{\text{ID}}={\mathit{\sigma }}_{\text{Ic}}-\mathrm{\Psi }{\mathit{\sigma }}_{\mathrm{I}}^{\prime }={\displaystyle \frac{c-\mathrm{\Psi }({\mathit{\sigma }}_{\text{II}}^{\prime }-\mathit{\mu }{\mathit{\sigma }}_{\mathrm{I}}^{\prime })}{\mathit{\mu }}},\text{(23)}& {\displaystyle }& {\displaystyle }{\mathit{\sigma }}_{\text{IID}}=\mathrm{0},\left(\text{by definition}\right).\end{array}$$

When $\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)

$$\begin{array}{}\text{(24)}& \mathit{\delta }\mathrm{\Psi }=\sqrt{{\left({\displaystyle \frac{\partial \mathrm{\Psi }}{\partial {\mathit{\sigma }}_{\mathrm{I}}^{\prime }}}\right)}^{\mathrm{2}}\mathit{\delta }{{\mathit{\sigma }}_{\mathrm{I}}^{\prime }}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial \mathrm{\Psi }}{\partial {\mathit{\sigma }}_{\text{II}}^{\prime }}}\right)}^{\mathrm{2}}\mathit{\delta }{{\mathit{\sigma }}_{\text{II}}^{\prime }}^{\mathrm{2}}},\end{array}$$

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{ϵ}{\mathit{\sigma }}_{\mathrm{I}}^{\prime }$, $\mathit{\delta }{\mathit{\sigma }}_{\text{II}}^{\prime }=\mathit{ϵ}{\mathit{\sigma }}_{\text{II}}^{\prime }$), we obtain

$$\begin{array}{}\text{(25)}& {\displaystyle }\mathit{\delta }\mathrm{\Psi }& {\displaystyle }=\sqrt{\begin{array}{c}\frac{{\mathit{\mu }}^{\mathrm{2}}}{(\mathrm{1}+\mathit{\mu }\tan (\mathit{\gamma }){)}^{\mathrm{2}}{\mathit{\sigma }}_{\text{II}}^{{}^{\prime }\mathrm{2}}}{\mathit{ϵ}}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{I}}^{{}^{\prime }\mathrm{2}}+\frac{(c-\mathit{\mu }{\mathit{\sigma }}_{\mathrm{I}}^{\prime }{)}^{\mathrm{2}}}{(\mathrm{1}+\mathit{\mu }\tan (\mathit{\gamma }){)}^{\mathrm{2}}{\mathit{\sigma }}_{\text{II}}^{{}^{\prime }\mathrm{4}}}{\mathit{ϵ}}^{\mathrm{2}}{\mathit{\sigma }}_{\text{II}}^{{}^{\prime }\mathrm{2}}\end{array}},\text{(26)}& {\displaystyle }& {\displaystyle }=\mathrm{\Psi }\mathit{ϵ}\sqrt{{\displaystyle \frac{{\mathit{\mu }}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{I}}^{{}^{\prime }\mathrm{2}}+(c-\mathit{\mu }{\mathit{\sigma }}_{\mathrm{I}}^{\prime }{)}^{\mathrm{2}}}{(c+\mathit{\mu }\tan (\mathit{\gamma }){\mathit{\sigma }}_{\text{II}}^{\prime }-\mathit{\mu }{\mathit{\sigma }}_{\mathrm{I}}^{\prime }{)}^{\mathrm{2}}}}},\text{(27)}& {\displaystyle }& {\displaystyle }=\mathrm{\Psi }\mathit{ϵ}R,\end{array}$$

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

$$\begin{array}{}\text{(28)}& \tan \left(\mathit{\gamma }\right)=-\mathrm{1}/\mathit{\mu },\end{array}$$

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 (-\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 }=\tan \mathit{\gamma }$, green line in Fig. 1b). The damage factor is then defined as

$$\begin{array}{}\text{(29)}& \mathrm{\Psi }=\left\{\begin{array}{ll}\frac{c+\mathit{\mu }\mathit{\gamma }{\mathit{\sigma }}_{\text{II}}^{\prime }-\mathit{\mu }{\mathit{\sigma }}_{\mathrm{I}}^{\prime }}{(\mathrm{1}+\mathit{\mu }\mathit{\gamma }){\mathit{\sigma }}_{\text{II}}^{\prime }},& \text{if }{\mathit{\sigma }}_{\mathrm{I}}^{\prime }<{\mathit{\sigma }}_{\text{II}}^{\prime }\tan \mathit{\gamma },\\ \frac{c}{{\mathit{\sigma }}_{\text{II}}^{\prime }+\mathit{\mu }{\mathit{\sigma }}_{\mathrm{I}}^{\prime }},& \text{otherwise}.\end{array}\right.\end{array}$$

Figure 2 Idealized domain for uniaxial compression simulations, with a solid boundary (Dirichlet conditions, $u=v=\mathrm{0}$) at the bottom and open boundaries (Neumann conditions, $\partial u/\partial n=\mathrm{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 (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² (corresponding to ∼ 20 m s⁻¹ winds or ∼ 0.33 m s⁻¹ 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⁻⁸ N m² (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

$$\begin{array}{}\text{(30)}& {\mathit{ϵ}}_{\text{asym}}={\displaystyle \frac{{\sum }_{i=a}^b{\sum }_{j=\mathrm{1}}^{n_y}\mathrm{|}({\mathit{\sigma }}_{\text{II}}{)}_{i,j}-({\mathit{\sigma }}_{\text{II}}{)}_{n_x-i,j}\mathrm{|}}{{\sum }_{i=a}^b{\sum }_{j=\mathrm{1}}^{n_y}\mathrm{|}({\mathit{\sigma }}_{\text{II}}{)}_{i,j}}},\end{array}$$

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

$$\begin{array}{}\text{(31)}& D=\sum _{i=a}^b\sum _{j=\mathrm{1}}^{n_y}{\displaystyle \frac{d_{i,j}^{t+\mathrm{\Delta }/\mathrm{2}}-d_{i,j}^{t-\mathrm{\Delta }/\mathrm{2}}}{\mathrm{\Delta }}}.\end{array}$$

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

$$\begin{array}{}\text{(32)}& \mathit{\theta }={\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}-{\displaystyle \frac{\mathit{\varphi }}{\mathrm{2}}}.\end{array}$$

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

$$\begin{array}{}\text{(33)}& \mathit{\theta }={\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}-{\displaystyle \frac{\mathit{\delta }}{\mathrm{2}}}.\end{array}$$

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 ± $\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.

Figure 3 (a) Damage (unitless), (b) ice thickness (m, colour) and velocity vectors (m s⁻¹), (c) mean normal strain rate invariant (${\dot{\mathit{ϵ}}}_{\mathrm{I}}$, d⁻¹), and (d) maximum shear strain rate invariant (${\dot{\mathit{ϵ}}}_{\text{II}}$, d⁻¹), 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 m², 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).

Figure 4 Scatter plots of local stress invariants (σ_I vs. σ_II, in kN m⁻¹, left column) and of the normal stresses and scaled strain rate invariants (σ_I vs. $(\mathrm{1}-d{)}^{\mathrm{3}}{\dot{\mathit{ϵ}}}_{\text{II}}$, right column) in heavily damaged (d>0.9) grid cells, at t=57 min (during the fracture development, top row), t=60 min (a few minutes after the fracture, middle row) and t=90 min (∼ 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. $(\mathrm{1}-d{)}^{\mathrm{3}}{\dot{\mathit{ϵ}}}_{\text{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)³), 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 R_max, 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⁻⁸ N m²). 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).

Figure 6 (a) Time evolution of the asymmetry factor ϵ_asym and (b) time series of the maximum error amplification ratio R_max, 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 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.

Figure 7 Sensitivity 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 8 Time 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-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.

Figure 9 Sensitivity 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 10 Time evolution of (a) the mean normal strain rate invariant integrated over the ice cover (d⁻¹) and (b) the maximum shear strain rate invariant integrated over the ice cover (d⁻¹), 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  Bolton, 1986).

Figure 11 Sensitivity 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 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. γ = ${\tan }^{-\mathrm{1}}({\mathit{\sigma }}_{\mathrm{I}}^{\prime }/{\mathit{\sigma }}_{\text{II}}^{\prime })$).

Figure 12 Time evolution of the mean normal strain rate invariant integrated over the ice cover (d⁻¹) using a stress correction path normal to the yield curve (γ=arctan (μ)) with α=3 (blue), α=1 and a longer viscous dissipation timescale (λ=10⁸ 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⁻²). 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⁵;  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⁵–10⁷, 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⁻¹), (c) mean normal strain rate invariant (${\dot{\mathit{ϵ}}}_{\mathrm{I}}$, d⁻¹), and (d) maximum shear strain rate invariant (${\dot{\mathit{ϵ}}}_{\text{II}}$, d⁻¹) 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 (c₀) field is defined locally at each grid cell by picking a random number between 7.0 and 13.0 kN m⁻². 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 (${\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.