the Creative Commons Attribution 4.0 License.

the Creative Commons Attribution 4.0 License.

# Landfast sea ice material properties derived from ice bridge simulations using the Maxwell elasto-brittle rheology

### Mathieu Plante

### Bruno Tremblay

### Martin Losch

### Jean-François Lemieux

The Maxwell elasto-brittle (MEB) rheology is implemented in the Eulerian finite-difference (FD) modeling framework commonly used in classical viscous-plastic (VP) models. The role of the damage parameterization, the cornerstone of the MEB rheology, in the formation and collapse of ice arches and ice bridges in a narrow channel is investigated. Ice bridge simulations are compared with observations to derive constraints on the mechanical properties of landfast sea ice. Results show that the overall dynamical behavior documented in previous MEB models is reproduced in the FD implementation, such as the localization of the damage in space and time and the propagation of ice fractures in space at very short timescales. In the simulations, an ice arch is easily formed downstream of the channel, sustaining an ice bridge upstream. The ice bridge collapses under a critical surface forcing that depends on the material cohesion. Typical ice arch conditions observed in the Arctic are best simulated using a material cohesion in the range of 5–10 kN m^{−2}. Upstream of the channel, fracture lines along which convergence (ridging) takes place are oriented at an angle that depends on the angle of internal friction. Their orientation, however, deviates from the Mohr–Coulomb theory. The damage parameterization is found to cause instabilities at large compressive stresses, which prevents the production of longer-term simulations required for the formation of stable ice arches upstream of the channel between these lines of fracture. Based on these results, we propose that the stress correction scheme used in the damage parameterization be modified to remove numerical instabilities.

The term landfast ice designates sea ice that is attached to the coastlines, acting as an immobile and seasonal extension of the land. It starts to form in shallow water in the early stages of the Arctic freeze-up (Barry et al., 1979; Reimnitz et al., 1978) and grows throughout the Arctic winter, usually reaching its maximum extent in early spring (Yu et al., 2014). Typically, large landfast ice areas can form and remain stable due to the presence of islands or via the grounding of ice keels on the ocean floor (Reimnitz et al., 1978; Mahoney et al., 2007; Selyuzhenok et al., 2017). Where the water is too deep for grounding, landfast ice can also form where ice floes are jammed in narrow passages between islands or pieces of grounded ice. In the Canadian Arctic Archipelago (CAA), this type of ice is referred to as landlocked. The resulting ice bridges, also called ice arches for their characteristic arching edges (Fig. 1), can have a profound influence on sea ice circulation, via the closure of gateways (Melling, 2002; Kwok, 2005), and on regional hydrography, via the formation of winter polynyas downstream of the arches (Barber and Massom, 2007; Dumont et al., 2010; Shroyer et al., 2015). Most studies about ice arches focus on the Nares Strait (Fig. 1) and Lincoln Sea ice bridges (Kozo, 1991; Dumont et al., 2008; Dansereau et al., 2017; Moore and McNeil, 2018; Vincent, 2019), which affect the export of thick multi-year ice into Baffin Bay (Kwok and Cunningham, 2010; Ryan and Münchow, 2017). Ice arches, however, are a seasonal feature in several locations of the Canadian Arctic Archipelago (Marko and Thomson, 1977; Sodhi, 1997; Melling, 2002) and are also present in the Kara and Laptev seas (Divine et al., 2004; Selyuzhenok et al., 2015; Olason, 2016), where they play a role in the formation of extensive landfast ice covers.

Despite decades of observations (Melling, 2002; Kwok, 2005; Moore and McNeil, 2018; Ryan and Münchow, 2017), the formation, persistence and breakup of ice arches remain difficult to predict. It is, however, clear from modeling studies that the ability of sea ice to form arches relates to the material properties of sea ice. A number of studies showed that ice arches are produced if the rheology includes sufficient material cohesion (Ip, 1993; Hibler et al., 2006; Dumont et al., 2008). Using the ellipse yield curve of Hibler (1979), this can be achieved either by decreasing the yield curve ellipse aspect ratio (Kubat et al., 2006; Dumont et al., 2008) or by extending the ellipse towards larger isotropic tensile strength (Beatty and Holland, 2010; Olason, 2016; Lemieux et al., 2016). The range of parameter values that are appropriate for the production of ice bridges varies between different numerical studies, suggesting that different forcing or model implementations may influence the ice arch formation (Olason, 2016; Lemieux et al., 2016, 2018).

In recent years, new rheologies were proposed to reproduce the observed characteristics of ice failure, such as the preferred orientation of the lines of fracture (Wilchinsky and Feltham, 2004; Schreyer et al., 2006) or the brittle behavior of sea ice at small scales (Girard et al., 2011; Dansereau et al., 2016). Among this effort, a brittle damage parameterization (Amitrano et al., 1999) was implemented in the neXtSIM model (Rampal et al., 2016), as part of the elasto-brittle (Girard et al., 2011, EB) and Maxwell elasto-brittle (Dansereau et al., 2016, MEB) rheologies. The MEB rheology was shown to produce ice arches in the Nares Strait region that remain stable for several days, and arch fractures that are part of the landfast ice breakup process (Dansereau et al., 2017). The simulated stable ice arches in Dansereau et al. (2017) are located downstream of either Smith Sound or Kennedy channel (see the orange curve in Fig. 1). These locations differ from the observed ice arch positions in Nares Strait upstream of these channels (see, e.g., Fig. 1) or in the Lincoln Sea (Vincent, 2019), which are well reproduced by standard VP or elastic viscous-plastic (EVP) models (e.g., Dumont et al., 2008; Rasmussen et al., 2010). Whether this difference in behavior stems from the different physics of MEB and VP rheologies or it is just due to the different numerical framework used in both models remains an open question.

The EB and MEB models have so far been implemented using Lagrangian advection schemes and/or finite-element methods (e.g., Rampal et al., 2016; Dansereau et al., 2017). These numerical features, however, make it difficult to compare the different MEB and EB physics with that of the standard VP or EVP rheologies of the modeling community, as these are usually implemented on Eulerian finite-difference (FD) numerical frameworks. In this paper, we present our implementation of the MEB rheology on the FD numerical framework of the McGill VP sea ice model (Tremblay and Mysak, 1997; Lemieux et al., 2008, 2014). To our knowledge, it is the first time the MEB rheology has been implemented on the numerical platform of a VP model such that its different physics can be assessed independently from the numerical implementation. With this model, we investigate the role of the damage parameterization and the material strength parameters in the formation of ice arches, using an idealized model domain capturing the basic features of real-life geometries where ice arches are observed. We also identify a numerical issue associated with the damage parameterization, which significantly impacts long simulations.

The paper is organized as follows. In Sect. 2, we present the implementation of the Maxwell elasto-brittle rheology in our FD numerical framework. A detailed analysis of the breakup of the ice bridge simulated by the MEB rheology is presented in Sect. 3, along with a sensitivity analysis of the results with respect to the material parameters. The MEB model performance in simulating compressive fractures is discussed in Sect. 4, with summarized conclusions in Sect. 5.

## 2.1 Momentum and continuity equations

The 2D momentum equation describing the motion of sea ice is written as follows:

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

**is the vertically integrated internal stress tensor, and**

*σ***($={\mathit{\tau}}_{\mathrm{a}}+{\mathit{\tau}}_{\mathrm{w}}$) is the total external surface forcings from winds and ocean currents. Note that we write the momentum equation in terms of the vertically integrated internal sea ice stresses (i.e., ∇⋅**

*τ**σ*) as standard in VP models (e.g., Hibler, 1979; Hunke and Dukowicz, 1997; Wilchinsky and Feltham, 2004), as opposed to the mean internal sea ice stresses (i.e., ∇⋅

**) used in previous implementations of the MEB rheology (Dansereau et al., 2016; Rampal et al., 2016). We assume no grounding of ice on the ocean floor and neglect the Coriolis term. This omission is appropriate for landfast ice but can result in small errors in drifting ice (Turnbull et al., 2017). The advection of momentum (which scales as**

*σ**ρ*

_{i}

*H*[

*U*]

^{2}∕

*L*, where

*H*, [

*U*] and

*L*are the characteristic ice thickness, velocity, and length scales) is 3 orders of magnitude smaller than characteristic air or ocean surface stresses (Zhang and Hibler, 1997; Hunke and Dukowicz, 1997). At the edge of an ice arch where a discontinuity in sea ice drift is present at the grid scale (2 km in our case), it remains 2 orders of magnitude smaller than other terms in the momentum equation.

The total surface stress is defined in terms of an effective stress (*τ*_{LFI}) that represents the combined wind and ocean forces acting on the landfast ice and an additional water drag term that only acts on the drifting ice. That is, using the standard bulk formula (with air and water turning angles set to zero, McPhee, 1979), we have the following equations:

where *ρ*_{a} and *ρ*_{w} are the air and water densities, *C*_{da} and *C*_{dw} are the air and water drag coefficients (see values in Table 2), and *u*_{a} and *u*_{w} are the surface air and water velocities. Note that the cross terms *u*_{w}** u** have been neglected. This equation is therefore exact for landfast ice, the focus of this study, and constitutes an approximation only for ice drifting over an ocean current. Below, we specify

*τ*_{LFI}and give the characteristic wind speed and ocean current equivalent to this forcing for reference.

The continuity equations used for the temporal evolution of the mean ice thickness *h* (volume per grid cell area) and concentration *A* ($\mathrm{0}<A<\mathrm{1}$) are written as follows:

where *S*_{h} and *S*_{A} are thermodynamic sink and source terms for ice thickness and compactness, respectively. As we are only interested in the dynamical behavior of the sea ice model, all thermodynamics are turned off so that *S*_{h}=0 and *S*_{A}=0. Mechanical redistribution (i.e., ridging) is taken into account by capping the ice concentration at 1 % (or 100 %) in convergence. As the mean ice thickness *h* is allowed to grow, the capping increases the actual ice thickness (Schulkes, 1995).

## 2.2 Rheology

### 2.2.1 Visco-elastic regime

Following Dansereau et al. (2016), we consider the ice as a visco-elastic brittle material behaving like a stiff spring and strong dashpot in series if the stresses are relatively small. The corresponding stress–strain relation is that of a Maxwell visco-elastic material:

where *λ* is the viscous time relaxation ($\mathit{\lambda}=\frac{\mathit{\eta}}{E}$, *η* being the vertically integrated viscosity), *E* is the vertically integrated elastic stiffness, **C** is the elastic modulus tensor and “:” denotes the double dot product of tensors. In generalized matrix form, the tensors **C** and $\dot{\mathit{\u03f5}}$ are written as follows:

where *ν* (=0.33) is the Poisson ratio. The components of the elastic modulus tensor **C** are derived using the plane stress approximation (i.e., following the original assumption that the vertical stress components are negligible; see Rice, 2010). Note that we neglect the advection of stress in the time derivative of Eq. (7) as we focus on landfast ice.

The visco-elastic regime of the MEB model (before fracture) is dominated by a fast and reversible elastic response (first term on the left-hand side of Eq. 7), with a slow viscous dissipation acting over longer timescales (second term on the left-hand side). The reversibility of the elastic deformations implies that the elastic strains return to zero if all loads are removed. This results from a memory of the previous elastic stress and strain states given by the time derivative in Eq. (7). The Maxwell viscosity term, although orders of magnitude lower than the other terms in the visco-elastic regime, leads to a slow viscous dissipation of this elastic stress memory over long timescales determined by *λ* (days in our case).

While Eq. (7) is similar in form to the stress–strain relationship of the EVP model (Hunke, 2001), the elastic component in the EVP model was introduced to improve the computational efficiency of the VP model by allowing for an explicit numerical scheme and efficient parallelization (Hunke and Dukowicz, 1997). In the MEB model, the elastic component represents the elastic behavior of sea ice, while the viscous relaxation component is introduced to dissipate the elastic strains into permanent deformations. The use of a viscous component is consistent with the observation of viscous creep (Tabata, 1955; Weeks and Assur, 1967) and viscous relaxation in field experiments (Tucker and Perovich, 1992; Sukhorukov, 1996; Hata and Tremblay, 2015b). The viscous relaxation term is also analogous to the viscous term in the thermal stress models of Lewis (1993) and Hata and Tremblay (2015a).

### 2.2.2 Damage parameterization

In the MEB model, the brittle fracture is simulated using a damage parameterization, which is based on progressive damage models originally developed in the field of rock mechanics to reproduce the nonlinear (brittle) behavior in rock deformation and seismicity (Cowie et al., 1993; Tang, 1997; Amitrano and Helmstetter, 2006). In these models, the material damage associated with microcracking is simulated by altering the material properties (e.g., the Young modulus or the material strength) at the model element (or local) scale. If heterogeneity is present in the material, the damage parameterization simulates the self-organization of the microcracks in a macroscopic line of fracture, as observed in laboratory experiments. It was first used for large-scale sea ice modeling by Girard et al. (2011) and is now implemented in the Lagrangian dynamic–thermodynamic sea ice model neXtSIM (Rampal et al., 2019).

The sea ice deformations associated with the brittle fractures are parameterized by a gradual decrease in the elastic stiffness *E* and viscosity *η* at the local scale, and consequently as a local increase in the magnitude of the deformation associated with a given stress state. The local increase in deformations results in the concentration of internal stresses in adjacent grid cells, leading to the propagation of the fractures in space. The decrease in elastic stiffness and viscosity is set by a damage parameter *d* representing the weakening of the ice upon fracturing (Bouillon and Rampal, 2015). The damage parameter has a value of 0 for undamaged sea ice and 1 for fully damaged ice.

The damage increases when the stress state exceeds a critical stress, defined by the Mohr–Coulomb criterion. This yield criterion is standard for granular materials and in agreement with laboratory experiments (Schulson et al., 2006) and field observations (Weiss et al., 2007). We also investigate the use of a compressive cutoff to limit the uniaxial compression (${\mathit{\sigma}}_{\mathrm{2}}={\mathit{\sigma}}_{\mathrm{I}}-{\mathit{\sigma}}_{\mathrm{II}}$; see Fig. 2). In terms of the stress invariants *σ*_{I} and *σ*_{II}, this can be written as follows:

where

where *σ*_{I} is the isotropic normal stress (defined as negative in compression), *σ*_{II} is the maximum shear stress, *c* is the vertically integrated cohesion, *μ* (=sin *ϕ*) is the coefficient of internal friction of ice, *ϕ* is the angle of internal friction and *σ*_{c} is the vertically integrated uniaxial compressive strength. The parameterization of *c* and *σ*_{c} follows the form of the internal sea ice pressure in the standard VP model with the ice concentration parameter *a* set to 20 (Hibler, 1979). The cohesion *c*_{0} and compressive strength ${\mathit{\sigma}}_{{\mathrm{c}}_{\mathrm{0}}}$ are the material properties derived from in situ observations (see Table 1 for values and references) and laboratory experiments (Timco and Weeks, 2010). Model parameters used in this and other studies are listed in Tables 2 and 3.

Following Rampal et al. (2016), the introduction of damage upon failure is proportional to the local stress in excess of the yield criterion. A damage factor Ψ ($\mathrm{0}<\mathrm{\Psi}<\mathrm{1}$) is used to scale the stress back on the yield curve. It is defined as follows (see Appendix A for the derivation Ψ):

where *σ*_{f} is the corrected stress lying on the yield curve and *σ*^{′} is the prior stress state that exceeds the yield criterion. Note that the stress components are all scaled by the same damage factor, such that the path of the stress correction in stress invariant space follows a line from the uncorrected stress state to the origin (see Fig. 2). The stress correction path does not correspond to a flow rule: the magnitude of the excess stress is only used to increase the damage parameter. It determines the magnitude of the strain associated with a stress state but otherwise does not change the visco-elastic relationship in Eq. (7).

The temporal evolution of the damage parameter follows a simple relaxation with a damage timescale *T*_{d} (Dansereau et al., 2016):

where *T*_{d} is set to the advective timescale associated with the propagation of elastic waves in undamaged ice (i.e., ${T}_{\mathrm{d}}=\mathrm{\Delta}x/{c}_{\mathrm{e}}$, Δ*x* being the spatial resolution of the model and *c*_{e} the elastic wave speed). Consequently, the damage at any given time is a function of the previously accumulated damage. No damage healing process was included in this study as we focus on the breakup of ice bridges at small timescales. For the same reason, the advection of damage is neglected. The relaxation timescale (${T}_{\mathrm{d}}/(\mathrm{\Psi}-\mathrm{1})$) in Eq. (14) is time step dependent via its dependency on the damage factor Ψ. That is, a larger time step yields larger stress increments and larger excess stresses at each time level, decreasing the timescale for the damage relaxation. The sensitivity of the damage parameterization on the model time step led Dansereau et al. (2016) to suggest that the model time step be set to exactly *T*_{d}, otherwise the damage could travel faster than the elastic waves. We argue that while this point is true when using a fixed damage reduction parameter (as in Amitrano et al., 1999; Girard et al., 2011), the use of a damage factor Ψ relates the damage parameter to the rate of changes in the stress state, which is associated with the propagation of elastic waves. The propagation of damage in space is thus bounded by the elastic wave speed, and a smaller time step (0.5 s in this study) should be used to respect the Courant–Friedrichs–Lewy (CFL) criterion associated to the elastic waves.

The elastic stiffness *E* and Maxwell viscosity *η* are written as a nonlinear function of *d*, with a dependency on the ice thickness and sea ice concentration inspired by the ice strength parameterization of Hibler (1979):

where *Y* (=1 GN m^{−2}, smaller than in Bouillon and Rampal, 2015, and similar to Dansereau et al., 2016; see Table 3) is the Young modulus of undeformed sea ice, *η*_{0} is the viscosity of undeformed sea ice and *α* is an integer set to 4 that determines the smoothness of the transition from visco-elastic behavior to the post-fracture viscous behavior (Dansereau et al., 2016). Note that *E* and *η* are defined as in previous implementations except for the linear dependence in ice thickness required because of the use of vertically integrated stress *σ*.

The relaxation time constant *λ* in Eq. (7) is then written as follows:

where *λ*_{0} ($={\mathit{\eta}}_{\mathrm{0}}/Y={\mathrm{10}}^{\mathrm{5}}$ s, smaller than in Dansereau et al., 2016, but in agreement with observations; see
Table 1) is a parameter that corresponds to the viscous relaxation timescale in undamaged sea ice. In the limit when *λ*_{0} tends to infinity, the MEB rheology tends to the elasto-brittle rheology (Girard et al., 2011).

Note that when a fracture is developing, the stress state is constantly brought back to the yield curve while the damage and the deformation increase. This is comparable to the plastic regime of the standard VP model of Hibler (1979): in the VP model, the nonlinear bulk and shear viscous coefficients reduce with increasing strain rates, such that the stress state (the product of the two) remains on the yield curve while the deformation increases. However, the plastic deformations in the VP model are defined by a normal flow rule, which also determines the orientation of the strain rate tensor (Bouchat and Tremblay, 2017; Ringeisen et al., 2019). In the MEB model, the large deformation associated to the damage is governed by the visco-elastic relationship of Eq. (7) and the yield curve does not directly determine the orientation of the strain rate tensor. The two models also differ post-fracturing: the VP model does not have a memory of past deformations other than via the continuity equation and its impact on the ice thickness and concentration. In the MEB rheology, the damage corresponds to a material memory of past deformations even if the thickness and concentration remain unchanged.

The nonlinear relationship of the viscous relaxation timescale on *d* and *A* ensures that the viscous term is very small in undamaged ice, and dominant in heavily damaged ice (see Eq. 7, where *λ* appears in the denominator). In this case, the deformations are large, irreversible and viscous. This is different from the standard VP and EVP models in which there is no change in the constitutive equation before or after the ice fracture. The dependency of *λ* on the ice concentration also ensures that the total stress tends toward zero for low concentration (i.e., in free drift) but not in a continuous (*A*∼1) but heavily damaged ice.

## 2.3 Numerical approaches

This model was coded using an Eulerian FD implicit numerical scheme and is the first implementation of the MEB model on the same numerical framework as the standard VP model. This implementation was motivated by the need for a direct comparison between the VP and the MEB rheologies independent of the different numerical approaches. It presents a significant change from previous implementations that use finite-element methods with a triangular mesh (Rampal et al., 2016; Dansereau et al., 2016) and/or Lagrangian advection schemes (Rampal et al., 2016). In the standard VP numerical framework, the stress components do not appear explicitly in the momentum equation. Instead they are written in terms of the nonlinear viscous coefficients and strain rates. For the MEB model, this is accomplished by treating the stress memory term from the time derivation of Eq. (7) as an additional forcing term. The damage parameterization is therefore the only new module to be coded.

### 2.3.1 Time discretization

The model equations are discretized in time using a semi-implicit backward Euler scheme. The uncorrected stress at time level *n* can then be written using Eq. (7) as follows:

where *n*−1 is the previous time level and where

Note that ${{\mathit{\sigma}}^{\mathbf{\prime}}}^{n}$ is a function of *σ*^{n−1}, which we refer to as the stress memory. Equation (18) is then substituted in the stress divergence term of Eq. (1), so that the *x* and *y* components of the momentum equation can be expanded as follows:

where *C*_{1}, *C*_{2} and *C*_{3} are the components of the tensor *C* (Eq. 8) and where the stress memory terms have been included in the forcing, i.e.,

The MEB rheology equations can then be implemented in a VP model by setting the VP bulk and shear viscosity to ${\mathit{\zeta}}_{\mathrm{VP}}=\mathit{\xi}\frac{{C}_{\mathrm{1}}+C\mathrm{2}}{\mathrm{2}}$ and *η*_{VP}=*ξ**C*_{3}, respectively, setting the pressure term *P*=0, and adding the stress memory terms.

The variables *E*^{n} and *λ*^{n} in Eq. (18) to (21) are discretized explicitly as follows:

using

### 2.3.2 Space discretization

The model equations are discretized in space using a centered finite-difference scheme on an Arakawa C grid. In this grid, the diagonal terms of the ** σ** and $\dot{\mathit{\u03f5}}$ tensors are naturally computed at the cell centers and the off-diagonal terms at the grid nodes. The

*x*component of the momentum equation is written as follows:

where

The shear terms in Eqs. (32) and (36) (${\dot{\mathit{\u03f5}}}_{xy}$, *ξ*_{z} and *γ*_{z}) are thus defined at the lower-left grid node rather than at the grid center. The staggering of the stress components is unavoidable when using the C grid and requires node approximations for the scalar values *h*, *A* and *d* (Losch et al., 2010). This is treated on our Cartesian grid with square cells by approximating the scalar prognostic variables at the nodes (*h*_{z}, *A*_{z} and *d*_{z}) using a simple average of the neighboring cell centers, i.e.,

and similarly for *A*_{z} and *d*_{z}. The stress–strain coefficients *ξ*_{z} and *γ*_{z} are then computed using (*h*_{z}, *A*_{z} and *d*_{z}) in Eqs. (15), (17) and (19).

The shear stress at the cell center must also be approximated when computing the stress invariants in the stress correction scheme (Eq. 13). Averaging the shear stress components from the neighboring nodes (as in Eq. 37 for the scalars) causes a checkerboard instability in the solution because of the staggered shear stress corrections and memories. To avoid this, the mean shear stress at the cell center is defined using an average of the neighboring shear stress increments (${\mathit{\xi}}_{z}^{n}{\dot{\mathit{\u03f5}}}_{xy}^{n}$), which are integrated in another shear stress memory term, defined at the grid center, i.e.,

where ${{\mathit{\sigma}}_{xy\phantom{\rule{0.25em}{0ex}}i,j}^{\prime}}^{n}{|}_{C}$ is the uncorrected shear stress at the grid center, ${\stackrel{\mathrm{\u203e}}{\left({\mathit{\xi}}_{z}^{n}{\dot{\mathit{\u03f5}}}_{xy}^{n}\right)}}_{i,j}$ is the shear stress increment averaged as in Eq. (37) and ${\mathit{\sigma}}_{xy\phantom{\rule{0.25em}{0ex}}i,j}^{n-\mathrm{1}}{|}_{C}$ is the corrected shear stress at the grid center from the previous time step. Note that the approximations in Eqs. (37) and (38) are required due to the use of a FD scheme, a notable difference with the other MEB implementations using finite-element methods (Dansereau et al., 2016; Rampal et al., 2019).

### 2.3.3 Numerical solution

With *n**x* tracer points in the *x* direction and *n**y* in the *y* direction, the spatial discretization on our C grid leads to a system of $N=\left(ny\right(nx+\mathrm{1})+nx(ny+\mathrm{1}\left)\right)$ nonlinear equations for the velocity components. By stacking all the *u* components followed by the *v* ones, we form the vector ** u** of size

*N*. The nonlinear system of equations (momentum) for

*u*^{n}and the other discretized equations (Eqs. 24–31) are solved simultaneously using an IMplicit–EXplicit (IMEX) approach (Lemieux et al., 2014). As described in the algorithm below, this procedure is based on a Picard solver (Lemieux et al., 2008), which involves an outer loop (OL) iteration. At each OL iteration

*k*, the nonlinear system of equations is linearized and solved using a preconditioned Flexible General Minimum RESidual method (FGMRES). The latest iterate

*u*^{k}is used to solve explicitly the damage and continuity equations. This iterative process is conducted until the L2 norm of the solution residual falls below a set tolerance of ${\mathit{\u03f5}}_{\mathrm{res}}={\mathrm{10}}^{-\mathrm{10}}$ N m

^{−2}. The uncorrected stresses ${{\mathit{\sigma}}^{\mathbf{\prime}}}^{n}$ is then scaled by the damage factor Ψ

^{n}and stored as the stress memory

*σ*^{n}for the following time step. This numerical scheme differs from that of Dansereau et al. (2017), who solve the equations using tracers (

*h*,

*A*,

*d*) from the previous time level.

- 1.
Start with initial iterate

*u*^{0},

do*k*=1,*k*_{max};- 2.
Linearize the nonlinear system of equations by using ${\mathit{u}}^{n,k-\mathrm{1}}$, ${h}^{n,k-\mathrm{1}}$, ${A}^{n,k-\mathrm{1}}$, and ${d}^{n,k-\mathrm{1}}$;

- 3.
Calculate

*u*^{n,k}by solving the linear system of equations with FGMRES; - 4.
Calculate ${\mathrm{\Psi}}^{n,k}=f\left({{\mathit{\sigma}}^{\prime}}^{n,k}\right)$;

- 5.
Calculate ${h}^{k,n}=f({h}^{n,k-\mathrm{1}},{\mathit{u}}^{n,k})$,

${A}^{n,k}=f({A}^{n,k-\mathrm{1}},{\mathit{u}}^{n,k})$,

${d}^{n,k}=f({d}^{n,k-\mathrm{1}},{\mathit{u}}^{n,k},{\mathrm{\Psi}}^{n,k})$; - 6.
Calculate ${E}^{k,n}=f({d}^{n,k},{h}^{n,k},{A}^{n,k})$,

${\mathit{\lambda}}^{n,k}=f({d}^{n,k},{h}^{n,k},{A}^{n,k})$; - 7.
If the Picard solver converged to a residual <

*ϵ*_{res}, stop;

enddo

- 2.
- 8.
Update the stress memory ${\mathit{\sigma}}^{n}={\mathrm{\Psi}}^{n}{{\mathit{\sigma}}^{\prime}}^{n}$.

A simple upstream advection scheme is used for *h*^{k,n} and *A*^{k,n} in step 5. Note that steps 4, 5, 6 and 8 are performed for all the grid points.

In the following, we present a series of idealized simulations to document the formation and break-up of ice arches with the MEB rheology and their sensitivity to the choice of mechanical strength parameters. Results from these simulations and observations are used to constrain the material parameters used in sea ice models. Here, we define an ice arch as the location of the discontinuity in the sea ice velocity (and later in the ice thickness and concentration fields) and the ice bridge as the landfast ice upstream of the ice arch.

Our model domain is 800 km ×200 km with a spatial resolution of 2 km (Fig. 3). The boundary conditions are periodic on the left and right, closed on the top, and open on the bottom. Two islands, separated by a narrow channel 200 km long and 60 km wide, are located 300 km away from the top and bottom boundaries. The initial conditions for sea ice are zero ice velocity, uniform 1 m ice thickness, 100 % concentration and zero damage. A southward forcing *τ*_{LFI} (see Eq. 4) is imposed on the ice surface and ramped up from 0 to 0.625 N m^{−2} (corresponding to 20 m s^{−1} winds or 0.33 m s^{−1} surface currents) in a 10 h period, a rate well below the adjustment timescale associated with elastic waves. The solution can therefore be considered as steady state at all times, which allows us to determine the critical forcing associated with a fracture event.

## 3.1 Control run

The breakup of landfast ice in our simulation proceeds through a series of fracture events that are highly localized in time (see Fig. 4) and space (see Figs. 5 and 6), separated by periods of elastic stress buildup (low brittle failure activity). Two major fracture events are seen in the simulation (Stage B and Stage D in Fig. 4). The first corresponds to the failure of ice in tension with the development of an ice arch on the downstream side of the channel (Fig. 5). The damage occurs on very short timescales (within minutes) and preconditions the formation of an arching flaw lead downstream of the ice bridge over longer timescales (Fig. 5b), in accordance with results from Dansereau et al. (2017). The second event corresponds to the collapse of the landfast ice bridge with the breakup of ice within and upstream of the channel (Fig. 6). As for the downstream ice arch, the lines of fractures are formed on short timescales and precondition the location of ridging on the advection timescale (Fig. 6b). The three remaining periods during which few new brittle fractures occur correspond to an elastic landfast ice regime (Stage A); a stable downstream ice arch regime (Stage C); and a drift ice regime when ice flows within, downstream, and upstream of the channel (Stage E).

In the first stage of the simulation, elastic stress builds up but remains inside the yield curve in the entire domain, such that there is no brittle failure activity (Fig. 4, Stage A). The sea ice in the elastic regime behaves as an elastic plate and deformations are linearly related to the internal stresses. The elastic stresses are determined by the orientation of the surface forcing with respect to the coastlines: there are large tensile stresses on the downstream coastlines, compressive stresses on the upstream coastlines and shear stresses on the four corners of the channel (Fig. 7). At the vertical line of symmetry (away from channel openings, Fig. 7a, dashed blue line), the simulated stress field is in good agreement with the analytical solutions from a 1D version of the momentum equation, giving us confidence in the numerical implementation of the model (see Appendix B and Fig. 8). Upstream and downstream of the channel, both stress invariants are important, reaching a maximum in magnitude at the channel corners and decreasing to a local minimum at the center of the channel. In this configuration, the second principal stress alignment (Fig. 7c) is along the *x* direction downstream of the coastlines (where the ice is in uniaxial tension) and along the *y* direction upstream of the coastlines (where the ice is in uniaxial compression). In the downstream end of channel, the second principal stress alignment follows the shape of an arch, transitioning to a vertical alignment towards the upstream channel entrance.

### 3.1.1 Downstream ice arch

The formation of the downstream ice arch is initiated at a surface forcing of ∼0.02 N m^{−2}. The initial fractures are located at the downstream corners of the channel where the stress state reaches the critical shear strength for positive (tensile) normal stresses. The fractures then propagate from these locations and form an arch (see Fig. 5a). The progression of the fracture into an ice arch is helped by the concentration of stresses at the channel corners and around the subsequent damage. That is, the damage permanently decreases the elastic stiffness, which leads to locally larger elastic deformations and increases the load in the surrounding areas, leading to the propagation of the fractures in space through regions where the internal stress state was originally subcritical. In the first order, the arching progression of the fracture from the channel corners follows the second principal stress direction (i.e., a failure in uniaxial tension on the plane perpendicular to the maximum tensile stress; see Fig. 7c). This differs from the expected angle of fracture in a Coulombic material of $\mathit{\theta}=\pm (\mathit{\pi}/\mathrm{4}-\mathit{\varphi}/\mathrm{2})$ from the second principal stress orientation (Ringeisen et al., 2019), as reported in Dansereau et al. (2019).

A second period of low brittle fracture activity follows the formation of the ice arch (period C in Fig. 4). In this stage, the ice downstream of the ice arch is detached from the land boundaries and starts to drift. The nonzero brittle fracture activity in this stage is due to the increased damage in regions of already damaged ice; since the local stress state lies on the yield curve, the increasing forcing constantly increases the stress states beyond the yield criterion, leading to further damage. Upstream of the ice arch, the elastic stresses show little change from Stage A, except for their increase in magnitude due to higher forcing (Fig. 9). As the yield parameters (*c*, *σ*_{c}) are not a function of the damage, tensile fracturing does not reduce the critical stress. This results in large tensile and shear stresses persisting along and north of the ice arch after its formation. The formation of a stress-free surface could be obtained by modifying the formulations of *c* and *σ*_{c0} such that they depend on the damage.

### 3.1.2 Ice bridge collapse

The second break-up event (Stage D in Fig. 4) corresponds to the fracture of ice upstream of the channel and the collapse of the ice bridge. The fractures are initiated at a surface forcing of 0.13 N m^{−2} on the upstream corners of the islands where the internal stress reaches the critical shear strength for negative (compressive) normal stresses. The propagation of damage from these locations is composed of two separate fractures (see Fig. 6a). First, a shear fracture progresses downstream along the channel walls, resulting in the decohesion of the landfast ice in the channel from the channel walls. The decohesion of the ice bridge increases the load on the downstream ice arch and on the landfast ice upstream of the channel. Second, a shear fracture propagates upstream from the channel corners at an angle 58^{∘} from the coastline. The shear fracture orientation corresponds to an angle *θ*=32^{∘} from the second principal stress orientation (Fig. 7c), which also deviates from the theoretical 22.5^{∘} in a granular material with *ϕ*=45^{∘} (Ringeisen et al., 2019).

Once the lines of fracture are completed, the ice bridge collapses and the ice in the channel starts to drift (Stage E). In this stage, landfast ice only remains in two wedges of undeformed ice upstream from the islands in which high compressive stress remains present (see Fig. 10a). The remaining continuous areas of undamaged ice drift downward into the funnel as a solid body with uniform velocity, with ridges building at the fracture lines. The ridge building is highly localized (see Fig. 6b) but slowly expands in the direction perpendicular to the lines of fracture. This follows from the increase in material strength with ice thickness, resulting in larger compressive stresses along the ridge such that the ice fracture occurs in the neighboring thinner ice, in a succession of fracture events that are localized in time (see peaks in Stage E in Fig. 4).

## 3.2 Sensitivity to mechanical strength parameters

The Mohr–Coulomb yield criterion defines the shear strength of sea ice as a linear function of the normal stress on the fracture plane. In stress invariant coordinates (*σ*_{I}, *σ*_{II}), this can be written in terms of two material parameters: the cohesion *c* and the coefficient of internal friction *μ*=sin *ϕ* (Fig. 2). The isotropic tensile strength (i.e., the tip of the yield curve) is then a linear function of the two (${\mathit{\sigma}}_{\mathrm{t}}=c/\mathit{\mu}$). In this section, we investigate the influence of these material parameters and of the use of a uniaxial compressive strength criterion on the simulated ice bridge.

### 3.2.1 Cohesion

Changing only the cohesion *c*_{0} (with a fixed internal angle of friction *ϕ*) moves the entire yield curve along the first stress invariant (*σ*_{I}) axis. For example, a higher cohesion increases the isotropic tensile strength ${\mathit{\sigma}}_{{\mathrm{t}}_{\mathrm{0}}}={c}_{\mathrm{0}}/\mathrm{sin}\mathit{\varphi}$ and also increases the shear strength uniformly for all normal stress conditions. In the ice bridge simulations, the choice of cohesion influences the critical surface forcing associated with the different stages of the simulations but does not change the series of events described in Sect. 3.1 or the orientation of the ice fractures. This is in agreement with results from Dansereau et al. (2017).

The critical surface forcing associated with the ice bridge breakup can be related to the cohesion using the 1D steady-state momentum equation (see Appendix B for details). Assuming an infinite channel running in the *y* direction, the shear stress along the channel walls (*σ*_{xy}) is given by

where *W* is the channel width (see Fig. 3). Using the yield criterion (Eq. 10) with *σ*_{I}=0 (i.e., *σ*_{II}=*c*), the maximum sustainable surface forcing *τ*_{LFIc} can be related to the cohesion as follows:

In the simulations, the critical forcing for the complete decohesion of ice bridges (point 5 in Figs. 4 and 6) with different widths follows the simple 1D model (Fig. 11). This indicates that although the fracture is initiated at a weaker forcing due to the concentration of stress at the channel corners, the ice arch sustains the increasing load such that the ice bridge remains stable.

Given that ice bridges and arches with a width of ∼60 km are frequent in the CAA (e.g., Nares Strait, Lancaster Sound or Prince Regent Inlet) and that the surface stresses regularly exceed 0.15 N m^{−2} (e.g., corresponding to a wind speed of 10 m s^{−1} or a tidal current of ∼0.15 m s^{−1}), this suggests a lower bound on the cohesion of sea ice of at least 5 kN m^{−1} (see yellow curve in Fig. 11). Similarly, the fact that the ice bridges are rarely larger than 100 km (some are seen intermittently in the Kara Sea, Divine et al., 2004) indicates that the cohesion of sea ice should be smaller than 10 kN m^{−1} (see red curve in Fig. 11). This range (5–10 kN m^{−1}) is lower than records from ice stress buoys measurements, which measure both thermal and mechanical internal stresses at smaller scales (40 kN m^{−2}, Weiss et al., 2007) but agree with estimates from ice arch observations (Sodhi, 1997). Note that higher forcing may be frequent in areas associated with strong tides, although these locations correspond to unstable landfast ice areas and recurrent polynyas (Hannah et al., 2009). Our estimates therefore provide a meaningful bound to be used in sea ice models.

### 3.2.2 Angle of internal friction

The angle of internal friction *ϕ*, analogous to the static friction between two solids, determines the constant of proportionality between the shear strength and the normal stress (*μ*=sin *ϕ*; see Eq. 10 and Fig. 2). Varying the angle of internal friction changes the shear strength of ice under tensile and compressive stresses in opposite ways: when increasing the angle of internal friction, the shear strength of ice in tension is reduced while that of ice in compression is increased (and vice versa). This affects the critical forcing associated with the downstream and upstream ice fractures. When decreasing *ϕ*, the downstream ice arch (Stage B) forms under a stronger forcing, and a weaker forcing is required for the development of the upstream lines of fracture. As such, while the cohesion determines the stability of the landfast ice in the channel, the collapse of the ice bridge also requires the uniaxial fracture of ice upstream of the channel, which is sensitive to the angle of internal friction. The angle of internal friction also determines the shape of the ice fractures: decreasing *ϕ* leads to an increase in the curvature of the downstream ice arch and intensifies the departure of the upstream lines of fracture from the *y* axis (see Fig. 12). The simulated orientations of the fracture lines (32 and 45^{∘} for *ϕ*=20 and 45^{∘}) differ from the orientations of 35 and 22.5^{∘} predicted
by the Mohr–Coulomb theory and do not vary linearly with the internal angle of friction.

### 3.2.3 Tensile strength

The yield curve modifications discussed above (varying *c*_{0} and *ϕ*) also change the tensile strength (both uniaxial and isotropic) of ice. The tensile strength determines the magnitude of the critical surface forcing necessary for the formation of the downstream ice arch (Stage B). The tensile stresses downstream from the islands can be approximated using the 1D version of the momentum equation as a function of the fetch distance *F*_{down} (see Fig. 3) from the islands to the bottom boundary of the domain (derivation in Appendix B):

This can be written as a function of the material parameters using a simplified Mohr–Coulomb criterion (Eq. 10) for the 1D case (Appendix B):

where $\mathit{\nu}=\mathrm{1}/\mathrm{3}$ was used. Substituting *σ*_{yy} from Eq. (41) into Eq. (42), the yield criterion can be written in terms of the surface forcing and the material parameters:

Using our cohesion estimates ($\mathrm{5}<c<\mathrm{10}$ kN m^{−1}), angles of internal friction in the range of observations (30 and 45^{∘}) and a typical surface forcing of 0.15 N m^{−2}, this would suggest stable bands of landfast ice of extent *F*_{down}∼6–13 km to be sustainable. This is similar to observations in the Arctic, where bands of landfast ice rarely exceed tens of kilometers unless anchor points are provided by stamukhi (Mahoney et al., 2014).

### 3.2.4 Compressive strength criterion

While not used in other MEB implementations (Dansereau et al., 2016, 2017), the compressive cutoff offers a limit on the simulated uniaxial compression, which can reach unrealistically large values and cause numerical instabilities (see Sect. 4). Including a compressive strength criterion (${\mathit{\sigma}}_{\mathrm{I}}-{\mathit{\sigma}}_{\mathrm{II}}>{\mathit{\sigma}}_{\mathrm{c}}$) can modify the upstream fracture event (Stage D) via the development of uniaxial compression fractures along the upstream coast of the islands if the uniaxial compressive stress upstream of the islands exceeds the ice strength typically observed in the field (∼40 kN m^{−2}; see Table 1). The critical surface forcing for the development of a compressive fracture can be approximated using the 1D version of the momentum equation. The maximum normal stress at the upstream coast of the islands is as follows:

where *F*_{up} is the distance between the top boundary of the domain and the upstream coasts of the islands (see Fig. 3). In the ideal case, the compression strength criterion is as follows:

The compression criterion can thus be written as a function of the surface forcing as follows:

Whether the ice will fail in shear (Mohr–Coulomb criterion) or in compression can be evaluated by substituting *τ*_{LFI} from Eq. (39) into Eq. 46, yielding the following criterion:

If this condition is met, the compression strength criterion does not influence the simulation, and the upstream shear fracture lines develop as in the control simulation (Fig. 13a). If the left-hand side of Eq. (47) is much smaller that *σ*_{c}, a compression fracture occurs before the ice bridge breakup and a ridge forms along the upstream coastlines, propagating in the channel entrance, while the ice in the channel remains landfast (Fig. 13b). If the terms are of a similar order, the decohesion of the ice bridge and the compression fractures are initiated simultaneously, such that the compression fracture occurs along the upstream coastlines but not in the channel entrance, as the ice starts to drift in and upstream of the channel (Fig. 13c).

In the Arctic, ice arches are commonly observed upstream of narrow channels, where granular floes jam when forced into the narrowing passage. This requires the ice to not to be landfast in the channel itself (Vincent, 2019), as opposed to the simulations presented above where the ice is initially landfast in the model domain. Contrary to results presented in Dansereau et al. (2017) where the presence of floes is simulated by a random seeding of weaknesses in the initial ice field, unstable ice arches upstream of the channel are not present in our simulations. Instead, our experiment simulates the propagation of ice fractures through the landfast ice upstream of a channel, which are akin to a failure in uniaxial compression (Dansereau et al., 2016; Ringeisen et al., 2019).

In theory, the angle of internal friction governs the intersection angle between lines of fracture (Marko and Thomson, 1977; Pritchard, 1988; Wang, 2007; Ringeisen et al., 2019). That is, the lines of fracture are oriented at an angle $\mathit{\theta}(=\mathit{\pi}/\mathrm{2}-\mathit{\varphi}/\mathrm{4}$) with the second principal stress direction, where the ratio of shear to normal stress is largest. In our simulations, the angles of fracture, although sensitive to the angle of internal friction, do not follow this theory. The fact that different angles of internal friction yield the same fracture orientation (e.g., for *ϕ*=20^{∘} and *ϕ*=30^{∘}; see Fig. 12) indicates that the orientation is not directly associated to the yield criterion in the MEB rheology (there is no flow rule in the MEB rheology). However, the orientation of the lines of fracture do have a sensitivity to the angle of internal friction, which suggests that the deformations are at least indirectly influenced by the yield criterion. This is in accord with previous results showing that the fracture orientation is determined by the concentration of stress along lines damage instability (Dansereau et al., 2019). This raises the question whether the lines of fracture may be influenced by the stress correction path used in the damage parameterization, which determines the stress state associated with the fractures. These questions are left for future work and will be addressed using a simple uniaxial loading numerical experiments (e.g., Ringeisen et al., 2019).

We speculate that in a longer simulation, ice would eventually jam between the upstream lines of fracture, resulting in the formation of a stable ice arch upstream of the channel. This is suggested by the orientation of the second principal stress component upstream of the channel (Fig. 10c). Longer-term simulations, however, are prevented by the presence of numerical instabilities associated with the current damage parameterization. As the integration progresses, the simulated fields lose their longitudinal symmetry at about the center line of the domain. This loss of symmetry occurs more rapidly as the residual norm increases (Fig. 14) and is not due to a difficulty in solving the equations: the nonlinear solver converges rapidly, within six iterations, given the small time step required by the CFL criterion to resolve the elastic waves. The errors are instead related to the integration of the residual norms in the model memory terms in the constitutive equation. The integrated error is only dissipated over a large number of time steps, such that the error in the solution is orders of magnitude larger than the set residual norm tolerance. This limits the current analysis to short-term simulations in which this issue remains negligible.

An error propagation analysis shows that the instabilities are largely attributed to the stress correction scheme and the computation of the damage factor Ψ (Eq. 13). Assuming that the model is iterated to convergence such that the uncorrected stress state has a relative error *ϵ*, the error on the corrected stress is as follows (see derivation in Appendix C):

where

If ${\mathit{\sigma}}_{\mathrm{I}}^{\prime}>\mathrm{0}$ (tensile stress state), $\mathrm{0}<R<\mathrm{1}$ (triangle inequality) and the error on the memory terms (*ϵ*_{M}) is of the same order as that of the uncorrected stress state ($\mathit{\u03f5}\le {\mathit{\u03f5}}_{\mathrm{M}}\le \sqrt{\mathrm{2}}\mathit{\u03f5}$). If ${\mathit{\sigma}}_{\mathrm{I}}^{\prime}<\mathrm{0}$ (compressive stress state), we have *R*≥1, the error on the stress memory can become orders of magnitude larger than that of the uncorrected stress state, and the model accuracy and convergence properties are greatly reduced. These errors are stored in the memory terms and accumulate at each fracture event. Note that as the elastic stress memory is dissipated over the viscous relaxation timescale, and this issue could be mitigated by decreasing the viscous coefficients *η*_{0}. Using a compressive strength cutoff capping also offers a limit to the uniaxial compression and reduces this instability. Another solution could be using a nonlinear yield curve which converges to the Tresca criterion (*σ*_{II}=const) for large compressive stresses (e.g., the yield criterion of Schreyer et al., 2006). We, however, argue that this issue in the damage parameterization should be treated by bringing the stress back onto the yield curve along a different path (e.g., following a line perpendicular to the curve). It might also be possible to use a different stress correction path to constrain the orientation of the lines of fractures to the yield criterion. This will be assessed in future work.

The MEB rheology is implemented in the Eulerian FD numerical framework of the McGill sea ice model. We show that the discretized Maxwell stress–strain relationship can be written in a form that resembles that of the VP model, with an additional memory term. The MEB rheology is then simply implemented by redefining the VP viscous coefficients in terms of the MEB parameters and by adding the damage parameterization in a separate module. To our knowledge, it is the first time the MEB rheology has been implemented in the same framework as a VP or EVP model. This will allow direct comparison of these models using the same numerical platform in future work.

In idealized ice bridge simulations, we show that the damage parameterization allows the ice fractures in the MEB model to propagate over large distances at short timescales. This process relies on the memory of the past deformations included in the model, which causes a concentration of stresses close to the preexisting damage. We also show that while the choice of yield curve influences the localization and orientation of the ice fractures, the angles of fracture propagation differ from those expected from the Mohr–Coulomb theory. This is consistent with results from Dansereau et al. (2019) that show that the fracture orientation is determined by the planes of damage instability. Preliminary results suggest that the orientation of the fracture lines are influenced by the stress correction scheme. This will be the subject of future work.

The stress correction scheme in the damage parameterization (Rampal et al., 2016) is also found to cause a problematic increase in the numerical errors in the stress memory terms. The growth of errors depends on the magnitude of the compressive stress associated with the ice failure. These errors accumulate in the memory term at each fracture event, creating numerical artifacts that dominate the solutions over time. We argue that this weakness of the damage parameterization should be treated as a numerical issue. In previous MEB implementations, asymmetries are expected due to either the asymmetric coastlines and forcing (Rampal et al., 2016) or the material heterogeneity used to initialize the model (Dansereau et al., 2016), such that this instability difficult to detect. A possible solution to this problem would be to use a nonlinear yield curve that converges to the Tresca criterion for large compressive stresses (e.g., the yield criterion of Schreyer et al., 2006). It may also be possible to eliminate this numerical noise by using a different stress correction scheme that does not follow a path to the origin. This will be assessed in future work.

The simulated breakup of the landfast ice bridge occurs with two main fracture events. First, an ice arch develops at the downstream end of the channel, shaping the edge of the ice bridge in the channel. This ice arch forms in all simulations and is stable in shape as long as the ice bridge remains in place, with a curvature that increases for smaller angles of internal friction. Second, shear fractures are formed at the upstream end of the channel, resulting in the decohesion of the channel ice bridge and the formation of landfast wedges upstream of the islands. Based on the simulation results, we determined that the parameterized cohesion most consistent to the observed ice bridges in the Arctic is in the range of 5–10 kN m^{−2}, lower than stress buoys that measure both dynamical and thermal stresses at smaller scales but in the range of values previously associated to ice arch observations.

Let ${\mathit{\sigma}}_{\mathrm{I}}^{\prime}$ and ${\mathit{\sigma}}_{\mathrm{II}}^{\prime}$ be the stress invariant at time level *n* before the correction is applied and *σ*_{If} and *σ*_{IIf} be the corrected stress invariant lying on the yield curve. Following Bouillon and Rampal (2015), we use a damage factor Ψ ($\mathrm{0}<\mathrm{\Psi}<\mathrm{1}$) to reduce the elastic stiffness and bring the stress state onto the yield curve, i.e.,

Substituting these relations into the Mohr–Coulomb criterion (${\mathit{\sigma}}_{\mathrm{IIf}}+\mathit{\mu}{\mathit{\sigma}}_{\mathrm{If}}=c$), we solve for Ψ:

Note that this relation implies that the stress correction is done following a line from the stress state (${\mathit{\sigma}}_{\mathrm{I}}^{\prime}$, ${\mathit{\sigma}}_{\mathrm{II}}^{\prime}$) to the origin (see Fig. 2). This scheme stems from applying the damage factor to each individual stress component. Other paths could be used for the correction (e.g., following a vertical or horizontal line) but would require the use of a different stress factor for the different components of the stress tensor. This could be used to cure the error propagation problem when large compressive stresses are present (see Appendix C).

Considering an infinite channel of landfast ice (** u**=0) along the

*y*direction with forcing

*τ*

_{LFY}=

*τ*

_{y}, we write the 1D steady-state momentum equation as follows:

where we have neglected the $\partial /\partial y$ terms. In this case, the normal stress is zero in the entire channel and the stress invariants are *σ*_{I}=0, *σ*_{II}=*σ*_{xy}. The shear stress at any arbitrary point *x* across the channel can be determined by integrating Eq. (B1) from the channel center (*x*=0) to *x*:

By symmetry, the maximum shear stresses in the channel are located at the channel walls at $x=\pm \frac{W}{\mathrm{2}}$, where *W* is the width of the channel. The maximum shear stress invariant on the channel walls is then

Similarly, we find the analytical solution for the normal stresses in a band of landfast ice with width *L*_{y} along an infinite coastline running in the *x* direction with a surface forcing *τ*_{LFI}=*τ*_{y}, by integrating the 1D momentum equation in which the $\partial /\partial x$ terms are neglected, i.e.,

Placing the landfast ice edge (where *σ*_{yy}=0) at *y*=0, the largest compressive stresses will be located along the coast, at $y=-{L}_{y}$. Note that in this case, shear stress is zero in the entire landfast ice and the stress invariants are function of both *σ*_{xx} and *σ*_{yy}:

This allows us to write the Mohr–Coulomb criterion in terms of *σ*_{yy}:

The error *δ**F* associated with a function $F(X,Y,Z,\mathrm{\dots})$ with uncertainties $(\mathit{\delta}x,\mathit{\delta}y,\mathit{\delta}z,\mathrm{\dots})$ is given by

In the damage parameterization, the components of the corrected stress tensor used as the memory terms (*σ*_{ijM}) can be written in terms of the uncorrected stress tensor (${\mathit{\sigma}}_{ij}^{\prime}$) and the damage factor Ψ (Eq. 13):

Using Eq. (A2), this can be rewritten in terms of the uncorrected stress invariants (${\mathit{\sigma}}_{\mathrm{I}}^{\prime}$, ${\mathit{\sigma}}_{\mathrm{II}}^{\prime}$):

Assuming that the model has converged to a solution within an error on the stresses $\mathit{\delta}{\mathit{\sigma}}_{ij}^{\prime}=\mathit{\u03f5}{\mathit{\sigma}}_{ij}^{\prime}$, $\mathit{\delta}{\mathit{\sigma}}_{\mathrm{I}}^{\prime}=\mathit{\u03f5}{\mathit{\sigma}}_{\mathrm{I}}^{\prime}$, $\mathit{\delta}{\mathit{\sigma}}_{\mathrm{II}}^{\prime}=\mathit{\u03f5}{\mathit{\sigma}}_{\mathrm{II}}^{\prime}$, where *ϵ* is a small number, the model convergence error propagates on the stress memory with an error of

Substituting ($\mathit{\delta}{\mathit{\sigma}}_{ij}^{\prime}$, $\mathit{\delta}{\mathit{\sigma}}_{\mathrm{I}}^{\prime}$, $\mathit{\delta}{\mathit{\sigma}}_{\mathrm{II}}^{\prime}$) for *ϵ* and using Eq. (C3), we obtain

or

Assuming that the error on the stress memory components (*ϵ*_{M}) has the form *δ**σ*_{ijM}=*ϵ*_{M}*σ*_{ijM}, we can express the relative error of the stress memory components as a function of the stress invariants as follows:

where

Our sea ice model code and output are available upon request.

MP coded the model, ran all the simulations, analyzed results and led the writing of the manuscript. BT participated in regular discussions during the course of the work and edited the manuscript. ML participated in regular discussions during MP research stays in Germany and edited the manuscript. JFL participated in the implementation of the MEB rheology and provided edits to the manuscript.

The authors declare that they have no conflict of interest.

Our sea ice model code and output are available upon request. We acknowledge the use of imagery from the NASA Worldview application (https://worldview.earthdata.nasa.gov, last access: 12 August 2019), part of the NASA Earth Observing System Data and Information System (EOSDIS).

Mathieu Plante would like to thank the Fonds de recherche du Québec – Nature et technologies (FRQNT) for the financial support received during the course of this work. Bruno Tremblay is grateful for the support from the Natural Science and Engineering and Research Council (NSERC) Discovery Program and the Office of Naval Research (N000141110977). This work is a contribution to the research program of Québec‐Océan and of the ArcTrain International Training Program.

This paper was edited by Daniel Feltham and reviewed by Jérôme Weiss and two anonymous referees.

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

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

Barber, D. and Massom, R.: Chapter 1 The Role of Sea Ice in Arctic and Antarctic Polynyas, in: Polynyas: Windows to the World, edited by: Smith, W. and Barber, D., vol. 74 of Elsevier Oceanography Series, 1–54, Elsevier, https://doi.org/10.1016/S0422-9894(06)74001-6, 2007. a

Barry, R., Moritz, R., and Rogers, J.: The fast ice regimes of the Beaufort and Chukchi Sea coasts, Alaska, Cold Reg. Sci. Technol., 1, 129–152, https://doi.org/10.1016/0165-232X(79)90006-5, 1979. a

Beatty, C. K. and Holland, D. M.: Modeling Landfast Sea Ice by Adding Tensile Strength, Am. Meteorol. Soc., 40, 185–198, https://doi.org/10.1175/2009JPO4105.1, 2010. a, b

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

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

Cowie, P. A., Vanneste, C., and Sornette, D.: Statistical physics model for the spatiotemporal evolution of faults, J. Geophys. Res.-Solid Earth, 98, 21809–21821, https://doi.org/10.1029/93JB02223, 1993. a

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

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

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

Divine, D. V., Korsnes, R., and Makshtas, A. P.: Temporal and spatial variation of shore-fast ice in the Kara Sea, Cont. Shelf Res., 24, 1717–1736, https://doi.org/10.1016/j.csr.2004.05.010, 2004. a, b

Dumont, D., Gratton, Y., and Arbetter, T. E.: Modeling the Dynamics of the North Water Polynya Ice Bridge, J. Phys. Oceanogr., 39, 1448–1461, https://doi.org/10.1175/2008JPO3965.1, 2008. a, b, c, d, e, f, g

Dumont, D., Gratton, Y., and Arbetter, T. E.: Modeling Wind-Driven Circulation and Landfast Ice-Edge Processes during Polynya Events in Northern Baffin Bay, J. Phys. Oceanogr., 40, 1356–1372, https://doi.org/10.1175/2010JPO4292.1, 2010. a

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

Hannah, C. G., Dupont, F., and Dunphy, M.: Polynyas and Tidal Currents in the Canadian Arctic Archipelago, Arctic, 62, 83–95, https://doi.org/10.14430/arctic115, 2009. a

Hata, Y. and Tremblay, L. B.: A 1.5-D anisotropic sigma-coordinate thermal stress model of landlocked sea ice in the Canadian Arctic Archipelago, J. Geophys. Res.-Oceans, 120, 8251–8269, https://doi.org/10.1002/2015JC010820, 2015a. a

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

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

Hibler, W. D., Hutchings, J., and Ip, C.: Sea-ice arching and multiple flow States of Arctic pack ice, Ann. Glaciol., 44, 339–344, https://doi.org/10.3189/172756406781811448, 2006. a

Hunke, E. C.: Viscous–Plastic Sea Ice Dynamics with the EVP Model: Linearization Issues, J. Comput. Phys., 170, 18–38, https://doi.org/10.1006/jcph.2001.6710, 2001. a, b, c, d, e

Hunke, E. C. and Dukowicz, J.: An Elastic – Viscous – Plastic Model for Sea Ice Dynamics, J. Phys. Oceanogr., 27, 1849–1867, 1997. a, b, c

Ip, C. F.: Numerical investigation of different rheologies on sea-ice dynamics, PhD thesis, 1, 242 pp., 1993. a

Kozo, T. L.: The hybrid polynya at the northern end of Nares Strait, Geophys. Res. Lett., 18, 2059–2062, https://doi.org/10.1029/91GL02574, 1991. a

Kubat, I., Sayed, M., Savage, S., and Carrieres, T.: Flow of Ice Through Converging Channels, Int. J. Offshore Polar, 16, 268–273, 2006. a

Kwok, R.: Variability of Nares Strait ice flux, Geophys. Res. Lett., 32, L24502, https://doi.org/10.1029/2005GL024768, 2005. a, b

Kwok, R. and Cunningham, G. F.: Contribution of melt in the Beaufort Sea to the decline in Arctic multiyear sea ice coverage: 1993–2009, Geophys. Res. Lett., 37, L20501, https://doi.org/10.1029/2010GL044678, 2010. a

Langleben, M.: Young's modulus for sea ice, Can. J. Phys., 40, 1–8, 1962. a

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

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

Lemieux, J.-F., Dupont, F., Blain, P., Roy, F., Smith, G. C., and Flato, G. M.: Improving the simulation of landfast ice by combining tensile strength and a parameterization for grounded ridges, J. Geophys. Res.-Oceans, 121, 7354–7368, https://doi.org/10.1002/2016JC012006, 2016. a, b, c, d, e

Lemieux, J.-F., Lei, J., Dupont, F., Roy, F., Losch, M., Lique, C., and Laliberté, F.: The Impact of Tides on Simulated Landfast Ice in a Pan-Arctic Ice-Ocean Model, J. Geophys. Res.-Oceans, 123, 1–16, https://doi.org/10.1029/2018JC014080, 2018. a

Lewis, J. K.: A model for thermally-induced stresses in multi-year sea ice, Cold Reg. Sci. Technol., 21, 337–348, 1993. a

Losch, M., Menemenlis, D., Campin, J.-M., Heimbach, P., and Hill, C.: On the formulation of sea-ice models. Part 1: Effects of different solver implementations and parameterizations, Ocean Modell., 33, 129–144, https://doi.org/10.1016/j.ocemod.2009.12.008, 2010. a

Mahoney, A., Eicken, H., Gaylord, A. G., and Shapiro, L.: Alaska landfast sea ice: Links with bathymetry and atmospheric circulation, J. Geophys. Res., 112, C02001, https://doi.org/10.1029/2006JC003559, 2007. a

Mahoney, A. R., Eicken, H., Gaylord, A. G., and Gens, R.: Landfast sea ice extent in the Chukchi and Beaufort Seas: The annual cycle and decadal variability, Cold Reg. Sci. Technol., https://doi.org/10.1016/j.coldregions.2014.03.003, 2014. a

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

McPhee, M. G.: The effect of the oceanic boundary layer on the mean drift of pack ice: application of a simple model, J. Phys. Oceanogr., 9, 388–400, 1979. a

Melling, H.: Sea ice of the northern Canadian Arctic Archipelago, J. Geophys. Res., 107, 3181, https://doi.org/10.1029/2001JC001102, 2002. a, b, c

Moore, G. W. K. and McNeil, K.: The Early Collapse of the 2017 Lincoln Sea Ice Arch in Response to Anomalous Sea Ice and Wind Forcing, Geophys. Res. Lett., 45, 8343–8351, https://doi.org/10.1029/2018GL078428, 2018. a, b

Olason, E.: A dynamical model of Kara Sea land-fast ice, J. Geophys. Res.-Oceans, 121, 3141–3158, https://doi.org/10.1002/2016JC011638, 2016. a, b, c, d, e, f, g, h

Pritchard, R. S.: Mathematical characteristics of sea ice dynamics models, J. Geophys. Res.-Oceans, 93, 15609–15618, https://doi.org/10.1029/JC093iC12p15609, 1988. a

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

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

Rasmussen, T., Nicolai, K., and Kaas, E.: Modelling the sea ice in the Nares Strait, Ocean Modell., 35, 161–172, https://doi.org/10.1016/j.ocemod.2010.07.003, 2010. a

Reimnitz, E., Toimil, L., and Barnes, P.: Arctic continental shelf morphology related to sea-ice zonation, Beaufort Sea, Alaska, Mar. Geol., 28, 179–210, 1978. a, b

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

Richter-Menge, J. A. and Elder, B.: Characteristics of pack ice stress in the Alaskan Beaufort Sea, J. Geophys. Res., 103, 21817–21829, 1998. a, b

Richter-Menge, J. A., McNutt, S. L., Overland, J. E., and Kwok, R.: Relating arctic pack ice stress and deformation under winter conditions, J. Geophys. Res.-Oceans, 107, SHE15-1–SHE15-13, https://doi.org/10.1029/2000JC000477, 2002. a

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

Ryan, P. A. and Münchow, A.: Sea ice draft observations in Nares Strait from 2003 to 2012, J. Geophys. Res.-Oceans, 122, 3057–3080, https://doi.org/10.1002/2016JC011966, 2017. a, b

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

Schulkes, R. M. S. M.: A Note on the Evolution Equations for the Area Fraction and the Thickness of a Floating Ice Cover, J. Geophys. Res., 100, 5021–5024, 1995. a

Schulson, E. M., Fortt, A. L., Iliescu, D., and Renshaw, C. E.: Failure envelope of first-year Arctic sea ice: The role of friction in compressive fracture, J. Geophys. Res.-Oceans, 111, C11S25, https://doi.org/10.1029/2005JC003235, 2006. a, b

Selyuzhenok, V., Krumpen, T., Mahoney, A., Janout, M., and Gerdes, R.: Seasonal and interannual variability of fast ice extent in the southeastern Laptev Sea between 1999 and 2013, J. Geophys. Res.-Oceans, 120, 7791–7806, https://doi.org/10.1002/2015JC011135, 2015. a

Selyuzhenok, V., Mahoney, A., Krumpen, T., Castellani, G., and Gerdes, R.: Mechanisms of fast-ice development in the south-eastern Laptev Sea: a case study for winter of 2007/08 and 2009/10, Polar Res., 36, 1411140, https://doi.org/10.1080/17518369.2017.1411140, 2017. a

Shroyer, E. L., Samelson, R. M., Padman, L., and Münchow, A.: Modeled ocean circulation in Nares Strait and its dependence on landfast-ice cover, J. Geophys. Res.-Oceans, 120, 7934–7959, https://doi.org/10.1002/2015JC011091, 2015. a

Sodhi, D. S.: Ice arching and the drift of pack ice through restricted channels, Cold Reg. Res. Eng. Lab. (CRREL) Rep. 77-18, 11 pp., 1997. a, b, c

Sukhorukov, K.: Experimental investigations of relaxation properties of sea ice internal stresses, in: The proceedings of the sixth (1996) international offshore and polar engineering conference, 354–361, 1996. a, b

Sulsky, D. and Peterson, K.: Toward a new elastic – decohesive model of Arctic sea ice, Physica D, 240, 1674–1683, https://doi.org/10.1016/j.physd.2011.07.005, 2011. a, b, c, d

Tabata, T.: A measurement of Visco-Elastic Constants of Sea Ice, J. Oceanogr. Soc. JPN., 11, 185–189, 1955. a, b, c, d

Tang, C.: Numerical simulation of progressive rock failure and associated seismicity, Int. J. Rock Mechan. Min. Sci., 34, 249–261, https://doi.org/10.1016/S0148-9062(96)00039-3, 1997. a

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

Tran, H. D., Sulsky, D. L., and Schreyer, H. L.: An anisotropic elastic-decohesive constitutive relation for sea ice, Int. J. Num. Analy. Methods Geomechan., 39, 988–1013, https://doi.org/10.1002/nag.2354, 2015. a, b, c, d

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

Tremblay, L.-B. and Mysak, L. A.: Modeling Sea Ice as a Granular Material , Including the Dilatancy Effect, J. Phys. Oceanogr., 27, 2342–2360, 1997. a, b, c

Tucker, W. B. and Perovich, D. K.: Stress measurements in drifting pack ice, Cold Reg. Sci. Technol., 20, 119–139, 1992. a, b, c

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

Vincent, R. F.: A Study of the North Water Polynya Ice Arch using Four Decades of Satellite Data, Sc. Rep., 9, 20278, https://doi.org/10.1038/s41598-019-56780-6, 2019. a, b, c

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

Weeks, W. F. and Assur, A.: The mechanical properties of sea ice, USA Cold Regions Research and Engineering Laboratory, Monograph II-C3, 80 p., 1967. a, b, c, d

Weiss, J., Schulson, E. M., and Stern, H. L.: Sea ice rheology from in-situ, satellite and laboratory observations : Fracture and friction, Earth Planet. Sci. Lett., 255, 1–8, https://doi.org/10.1016/j.epsl.2006.11.033, 2007. a, b, c, d, e, f

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

Yu, Y., Stern, H., Fowler, C., Fetterer, F., and Maslanik, J.: Interannual Variability of Arctic Landfast Ice between 1976 and 2007, J. Climate, 27, 227–243, https://doi.org/10.1175/JCLI-D-13-00178.1, 2014. a

Zhang, J. and Hibler, W. D.: On an efficient numerical method for modeling sea ice dynamics, J. Geophys. Res., 102, 8691–8702, 1997. a

- Abstract
- Introduction
- Maxwell elasto-brittle model
- Results
- Discussion
- Conclusions
- Appendix A: Damage factor Ψ
- Appendix B: Analytical solutions of the 1D momentum equation
- Appendix C: Error propagation analysis
- Code availability
- Author contributions
- Competing interests
- Acknowledgements
- Financial support
- Review statement
- References

- Abstract
- Introduction
- Maxwell elasto-brittle model
- Results
- Discussion
- Conclusions
- Appendix A: Damage factor Ψ
- Appendix B: Analytical solutions of the 1D momentum equation
- Appendix C: Error propagation analysis
- Code availability
- Author contributions
- Competing interests
- Acknowledgements
- Financial support
- Review statement
- References