A generalized stress correction scheme for the MEB rheology: impact on the fracture angles and deformations

A generalized damage parameterization is developed for the Maxwell Elasto-Brittle (MEB) rheology that reduces the growth of residual errors associated with the correction of super-critical stresses. In the generalized parameterization, a decohesive stress tensor is used to bring the super-critical stresses back on the yield curve based on any correction path. The sensitivity of the simulated material behaviour to the magnitude of the decohesive stress tensor is investigated in uniaxial compression simulations. Results show that while the decohesive stress tensor influences the short-term fracture deformation 5 and orientation, the long-term post-fracture behaviour remains unchanged. Divergence first occurs when the elastic response is dominant followed by post-fracture shear and convergence when the viscous response dominates – contrary to laboratory experiment 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. Using the generalized damage parameterization together with a stress correction path normal to the yield curve brings the simulated fracture angles 10 closer to observations (from 40− 50◦ to 35− 45◦, compared to 20− 30◦ in observations) and reduces the growth of errors sufficiently for the production of longer-term simulations.

thus calls for the development of high-resolution sea ice forecast products that capture the finer-scale lead structures (Jung et al., 2016). 25 As sea ice models are moving to higher spatial resolutions, they become increasingly capable of resolving LKFs Bouchat and Tremblay, 2020). The simulation of the ice fractures yet represents a challenge. To this day, most sea ice models simulate the motion of sea ice using plastic rheologies or modifications thereof (Hibler, 1979;Hunke, 2001).
While several improvements were made on the numerics and efficiency of the methods used to solve the highly non-linear momentum equation (Hunke, 2001;Lemieux et al., 2008Lemieux et al., , 2014Kimmritz et al., 2016;Koldunov et al., 2019), the physics 30 governing the ice fracture remains mostly the same. A number of rheologies have however been developed over the years in an attempt to simulate the observed sea-ice deformations (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). Among these new approaches, a damage parameterization derived for rock mechanics and seismology models (Amitrano et al., 1999;Amitrano and Helmstetter, 2006) was adapted for the large scale modelling of sea ice (Girard et al., 2011;Bouillon and Rampal, 2015). 35 This parameterization uses a damage parameter to represent the changes in material properties associated with fractures. While still based on the continuum assumption, it allows for fractures to propagate on short time-scales in the sea-ice cover. It is used in the Elasto-Brittle (EB Bouillon and Rampal, 2015;Rampal et al., 2016) and Maxwell Elasto-Brittle (MEB Dansereau et al., 2016) rheologies, implemented in the large scale sea-ice Finite Element model neXtSIM (Rampal et al., 2019) and, recently, in the Finite Difference McGill sea ice model . 40 The damage parameterization is relatively new, and it remains unclear to what extent differences in material behaviour are associated with the damage or to other rheological parameters. One known difference is the fracture development associated with local damage, stress concentration and damage propagation, rather than prescribed by an associative normal flow rule as in the standard VP models. The fracture angle simulated by the MEB and standard VP models are nonetheless in the same range (θ = 35 − 55 • , Dansereau et al., 2019;Hutter et al., 2020), which is larger than those derived from high-resolution satellite 45 observations (θ = 20 − 45  and in-situ observations (θ = 20 − 30 • Marko and Thomson, 1977;Schulson, 2004). In the standard VP model, modifications of the mechanical strength parameters (compressive and shear) and the use of non-associated flow rules lead to smaller fracture angles that are more in line with observations (Ringeisen et al., 2019(Ringeisen et al., , 2020. In the MEB rheology, the fracture angles can be reduced by increasing the angle of internal friction or the Poisson ratio (Dansereau et al., 2019). These sensitivities suggest that modifications to the damage parameterization could be used to bring 50 the simulated fracture angles closer to observations, but has not yet been tested.
The MEB rheology also presents some numerical challenges associated with the growth of residual errors associated with the damage parameterization at the grid scale . These errors can be attributed to the stress correction scheme, a numerical tool used to define the growth of damage and to bring the super-critical stresses back to the yield curve. Other progressive damage models instead represent the damage parameter as a discrete function of the number of failure cycles 55 (Main, 2000;Amitrano and Helmstetter, 2006;Carrier et al., 2015). In continuum damage mechanics, a damage potential derived from thermodynamic laws (Murakami, 2012) is used to simulate the 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 (Schreyer et al., 2006;Sulsky and Peterson, 2011).
In this paper, we present a generalization of the damage parameterization that reduces the growth of the residual errors asso-60 ciated with the stress correction and brings the simulated fracture angle of sea ice in simple uniaxial loading experiments closer to observations. Inspired by the work of Schreyer et al. (2006) and (Sulsky and Peterson, 2011), we introduce a decohesive stress associated with the fracture of sea ice and test its influence on the simulated sea-ice fracture and deformations in uniaxial loading experiments.
The paper is organised as follows. In section 2, we present the MEB rheology and governing equations. The generalized 65 stress correction scheme is described in section 3. The uniaxial loading experimental set-up is presented in section 4 along with the definition of diagnostics used to quantify the growth of damage and the growth of residual errors. Results are presented in section 5, with a focus on the material behaviour in uniaxial compression experiments and its response to the changes in the damage parameterization. In section 6, we discuss the influence of the stress correction and seeded heterogeneity. Conclusions are summarized in section 7.  Tremblay and Mysak, 1997;Lemieux et al., 2008;Plante et al., 2020). The vertically integrated 2D momentum equation for sea ice, forced with surface friction only (i.e. ignoring the sea surface tilt, the coriolis and the ice 75 grounding terms), can be written as: where ρ i is the ice density, h is the mean ice thickness, u (= uî + vĵ) 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 80 model domain and forcing (Turnbull et al., 2017). Following , we define the uniaxial loading by a surface wind stress τ a and prescribe an ocean at rest below the ice: where ρ w is the water density, C dw is the water drag coefficient and u is the sea ice velocity (see values in Table 1).
The prognostic equations for the mean ice thickness h (volume per grid cell area) and concentration A are written as: where the thermodynamic source an sink terms are ignored.

90
In the MEB rheology, the ice behaves as a visco-elastic material with a fast elastic response and a viscous response over a longer-time scale. The governing equation for this visco-elastic material can be written as (Dansereau et al., 2016(Dansereau et al., , 2017Plante et al., 2020): where E is the elastic stiffness defined as the vertically integrated Young Modulus of sea ice, λ is the viscous relaxation time-95 scale, C is the elastic tensor (fourth order), ":" denotes the inner double tensor product and˙ is the strain rate tensor. The elastic tensor C and strain rate tensor˙ can be written is matrix form 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 parameter ruling the dependency of the material strength properties on sea-ice concentration (Hibler, 1979;Rampal 110 et al., 2016), λ 0 (= 10 5 s, ≈1 day) is the viscous relaxation time scale for undamaged sea ice and α is a parameter ruling the post-fracture transition to the viscous regime.

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

115
where σ I is the isotropic normal stress invariant (compression defined as negative), σ II is the maximum shear stress invariant, is the angle of internal friction, and c is the vertically integrated cohesion, defined as: where c 0 (= 10 kN m −2 ) is the cohesion of sea ice derived from observations (Sodhi, 1997;Tremblay and Hakakian, 2006; 120 Plante et al., 2020) or laboratory experiments (Timco and Weeks, 2010). No compressive or tensile strength cut-off are used in this analysis. The reader is referred to Table 1 for a list of default model parameters.

Damage parameterization
The prognostic equation for the damage parameter d in the standard MEB rheology is parameterized using a relaxation term with time scale T d (= 1 s) as: where 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 advection are neglected as we are focusing on the ice fracture.
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. 12: 3 Generalized stress correction 135 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). To this end, we chose to define the damage factor in terms of the shear stress invariant only, as: where the critical shear stress invariant σ IIc is defined by the intersection point between the correction path and the yield curve (see Fig 1). After some algebra, we obtain: 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: 145 In this manner, the correction of super-critical stresses can follow any line in the stress invariant space provided that the damage increases when ice fractures (Ψ < 1, or γ < 90 • ). The generalized formulation now allows for the use of a yield curve without cohesion (c = 0 kN m −1 ), something that is not possible in the standard parameterization otherwise Ψ is identically equal to 0 (see Eq. 13).
Note that using a stress correction path other than the standard path to the origin means that the corrected normal stress 150 differs from the scaled super-critical stress Ψσ I . We define this difference as the decohesive stress tensor needed to for the corrected stress to follow the stress correction path γ (see Fig. 1). The stress correction equation (Eq. 14) then becomes: The invariants of the decohesive stress tensor (σ ID ,σ IID ) are therefore written as: When tan γ = σ I /σ II and σ ID = σ IID = 0, we obtain the standard damage parameterization of Dansereau et al. (2016) as a special case where the stress correction path depends on the super-critical stress state. 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 that they both determine the change in stress state associated with the development of a fracture. In the present 160 scheme, σ D is derived from the stress correction path, while the decohesive strain rate in Schreyer et al. (2006) is derived from the opening of a lead based on granular theory.

Projected error
The error δΨ on the damage factor Ψ(σ I , σ II ) can be written as: where (δσ I ,δσ II ) are the errors on the calculated stress invariants. Expanding the derivative terms (using Eq. 18) 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.

170
Assuming that the uncorrected stress is close to the yield criterion (i.e. σ II + µσ I − c ∼ 0), this relation indicates that the error amplification ratio R goes to infinity if: 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 175 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, correspond 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. 22). This behaviour is opposite to that of the standard stress correction scheme, which has small R values 180 in tension and large values in compression . To minimize the errors for all stress states, we blend the two schemes (i.e. Eq. 17 in compression and Eq. 13 in tension, see Fig. 1b). We set the transition between the two schemes at the points where they are both equal (i.e., at σ I /σ II = tan γ, see green line in Fig 1b). The damage factor is then defined as:  (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 .
A solution to the non-linear momentum and constitutive equations (Eqs. 1 and 5) is found using a Picard solver. The Picard   We monitor the growth of the residual error in the simulations using a normalised domain-integrated asymmetry factor ( asym ) in the maximum shear stress invariant field (σ II ), defined as:

Damage activity
We define the damage activity D as the total damage integrated over the original ice domain in a 1 minute interval:

215
This parameter is analog to the damage rate in (Dansereau et al., 2016(Dansereau et al., , 2017. Note that this definition of damage activity (or damage rate) emphasizes activity in undamaged ice and is not sensitive to activity in already heavily damaged ice.

Fracture angle
When loaded in uniaxial compression, a granular material fails in diamond-shaped shear fractures (e.g. see Marko and Thomson, 1977;Ringeisen et al., 2019). We define the fracture angle θ as the angle between the y-axis and the fracture lines (see  Fig. 2). The orientation of these fracture lines have been measured in laboratory using in uniaxial loading experiments. Several theories were developed to relate the fracture angle in terms of material parameters. The most common is the Mohr-Coulomb theory (Coulomb, 1773;Mohr, 1900), where the fracture angle is related to the angle of internal friction as: This theory tends to underestimate the fracture angle of granular materials in laboratory experiments (Bardet, 1991). In the 225 Roscoe (1970) theory, the fracture angle is defined instead in terms of the angle of dilatancy (δ) of the granular material: If δ = φ, the two theories give the same fracture angle θ. In general, the fracture angle 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 230 ± tan(W/L) ∼ ±2 • , where W is the fracture width (typically a few grid cells wide in our results, 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 well defined.

Control simulation: standard damage parameterization
In the control simulation, a pair of conjugate fracture lines first appear when the surface forcing τ a = 0.29 N/m, along with 235 secondary fracture 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 fracture lines 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 in the high range seen in observations (θ = ∼20-40 • Marko and Thomson, 1977;Hibler III and Schulson, 2000;Schulson, 2004;Hutter et al., 2020).
When the ice fractures, the initial response is mostly elastic with divergence along the fracture line. The resulting stress concentration influences the propagation of the fracture in space over short time-scales (seconds) governed by the elastic waves speed. The sea-ice deformation continues to occur post-fracture in the damaged ice and, over time, the response transitions 245 from elastic to viscous-dominated as the Maxwell viscosity dissipates the elastic stresses and creates permanent viscous deformations. 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 for instance 4 c,f,i). The simulation reaches steady state with deformations that are fully viscous and localized in the heaviest damage areas (Fig. 4g-i). This causes a predominance of shear and convergence deformation along the fracture line throughout the simulation. The asymmetries in the solution are very small at the beginning of the simulation (t ≤ 57min), and do not grow until fractures occur (Fig. 5a-b). As the fractures 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 −6 N m −2 ). The growth in asym are associated with large values of damage error amplification ratio R (reaching ∼20, Fig. 5b). Since asym is a domain-integrated quantity, it increases in time following large local error growths R. This 255 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 indicate the level of amplification of residual errors in the simulations, at times by more than one order of magnitude locally (R max > 10).

Generalized stress correction 260
The generalized damage parameterization reduces the growth of residual errors, with decreasing error amplification ratio R max for increasing path angle γ (Fig. 6a). This results in an overall reduction of the asymmetry factor asym (Fig. 6b), allowing for the production of longer-term simulations that include post-fracture deformations. This improvement is only significant when using γ > 0. For γ < 0, the maximum error amplification ratio R max remains important with periods when the residual error increases by up to two orders of magnitude locally.

265
Results show that the fracture angle is sensitive to the decohesive stress tensor, with decreasing fracture angle θ for increasing stress correction path angle γ (Fig. 7). This finding is in line with results from Dansereau et al. (2019), where the fracture angle was related to the far-field stress associated with the collective damage. In the MEB model, the far-field stresses directly depends on corrected stress state, including σ D in the generalized damage parameterization. Increasing the correction path angle γ reduces the fracture angles in better agreement to observations.

270
Along the fracture lines, the correction path angle γ influences the time-integration required to reach the same damage and deformation rates (Fig. 8). This 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 followed by a transition to viscous deformations where shear and convergence defor-275 mations are predominant (Fig. 8a). In contrast with results from the VP model and from typical granular material behaviour, divergent post-fracture deformation is only present when tensile stresses develop, e.g. at the intersection between conjugate lines of fracture.

Sensitivity to φ and ν
Repeating the experiment using different angles of internal friction (φ) shows that the fracture angle decreases with increasing 280 φ. The simulated fracture 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 fracture angle 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 approaches the horizontal, fracture yields large stress corrections but small damage increases (i.e., Ψ = 1), such that the angle of fracture is 285 governed by the stress correction and is weakly sensitive other model parameters. Based on these results, we suggest the use of a correction path that is normal to the yield criterion (γ = arctan µ, see black points in Fig. 9).
Decreasing the angle of internal friction reduces the shear strength of sea ice for a given normal stress, such that the fracture develops earlier in the simulation (i.e. under smaller surface forcing, Fig. 10). It also reduces the divergence associated with the elastic response when ice fractures and increase the convergence in the post-fracture viscous regime. This result is typical 290 for granular material, with smaller fracture angles associated with larger angles of dilatancy and divergence during the fracture development.
The fracture angle 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 also in Dansereau et al., 2019). Note that the 295 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 fracture angle is not affected by the changes in Poisson ratio thus indicates that the stress concentration and propagation of the fracture in space is mainly controlled by the stress correction rather than by the relaxation of material properties with damage. We speculate that the sensitivity of the fracture angle to the Poisson ratio in the standard stress correction scheme stems from the dependency of the stress correction path angle to the 300 super-critical stress state (i.e. γ = tan −1 (σ I /σ II )).

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 accu-305 mulate and influence the temporal evolution of the solution. In regions of heavily damaged ice, the accumulated 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 fracture angles, without significantly impacting the 310 material deformations. Using a large correction path angle γ (> 45 • ), however, significantly slows the damage production and reduces the simulated sensitivity of the fracture angle 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 fracture 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. Nemat -Nasser, 1993). This indicates that the decohesive stress tensor cannot be used to influence the deformations associated to the fracture of ice in the MEB rheology unless other parameterizations, such as including a decohesive strain tensor during the fractures (e.g., see Schreyer et al., 2006;Sulsky and Peterson, 2011), are added to the rheology.
The viscous dissipation timescale (λ) in our model is set based on observations (∼ 10 5 , Tabata, 1955;Hata and Tremblay, 2015), and is one order of magnitude smaller than in other MEB implementations (Dansereau et al., 2016;  Note that the results presented above neglect heterogeneity in the ice cover, a factor that is responsible for much of the brittle material behaviour in progressive damage models (Amitrano and Helmstetter, 2006). Heterogeneity was neglected in 340 the analysis above to isolate the growth of the residual errors. While including heterogeneities does not change the overall physics and sensitivity to the damage parameterization, it creates irregular sliding planes instead of the linear diamond shape fractures (Fig. 13a), naturally creating contact points where ridging occurs with lead opening elsewhere along the fracture lines.
This results in a form of granular dilatancy typical of granular materials.

345
We propose a generalized stress correction scheme for the damage parameterization to reduce the growth of residual errors in the MEB sea ice model. 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.
The sensitivity of the fracture angles and sea-ice deformations to these changes are investigated in the context of the uniaxial compression experiment similar to those presented in Ringeisen et al. (2019).

350
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 fractures, contrary to laboratory experiments of granular materials and satellite observations of sea ice. The use of a decohesive stress tensor influences the fracture angle of sea ice, but does not influence the type of deformation rates (convergence and shear), nor the simulated dilatancy. Future work will involve the modification of the non-linear relationship between the Maxwell viscosity and the 355 damage. We also show that the sensitivity of the fracture angle to the Poisson ratio, seen when using the standard damage parameterization, disappears when using the generalizes stress correction scheme with a fixed stress correction path. This suggests that in the MEB model, the stress concentration and fracture propagation is 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 360 reduces the growth of residual errors and allows for the production of longer term simulations with post-fracture deformations.
Using this stress correction path also reduces the fracture angles by ∼5 • , bringing them in the range of 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.