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

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.


Introduction
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 20 landfast ice areas can form and remain stable due to the presence of islands or by the grounding of ice keels on the ocean floor 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., 2008Lemieux et al., , 2014.
To our knowledge, it is the first time the MEB rheology is 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 60 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 organised as follows. In section 2, we present the implementation of the Maxwell Elasto-Brittle rheology in our FD numerical framework. A detailed analysis of the break up of the ice bridge simulated by the MEB rheology is presented in 65 section 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 section 4, with summarized conclusions in section 5.

Momentum and continuity equations
The 2D momentum equation describing the motion of sea ice is 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 τ (= τ a + τ 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 75 stresses (i.e., ∇ · (hσ)) 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 three orders of magnitude smaller than a characteristic air or ocean surface stresses (Zhang and Hibler, 1997;Hunke and Dukowicz, 1997). At the edge of an ice arch 80 where a discontinuity in sea ice drift is present at the grid scale (2 km in our case), it remains two orders of magnitude smaller than other terms in the momentum equation.
The total surface stress is defined in terms of an effective stress (τ LF I ) 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 τ = ρ a C da |u a |u a + ρ w C dw |u w − u|(u w − u), ≈ ρ a C da |u a |u a − ρ w C dw |u w |u w − ρ w C dw |u|u, 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 90 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 τ LF I 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 (0 < A < 1) are written as: 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.

100
As the mean ice thickness h is allowed to grow, the capping increases the actual ice thickness (Schulkes, 1995).

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 dash-pot in series if the stresses are relatively small. The corresponding stress-strain relation is that of a Maxwell visco-elastic 105 material: where λ is the viscous time relaxation (λ = η 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˙ are written as: 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 for instance Rice 115   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 time scales (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 120 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 Elastic Viscous Plastic (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 (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).

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 non-linear (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) 135 at the model element (or local) scale. If heterogeneity is present in the material, the damage parameterization simulates the self-organisation 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 140 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.

145
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 cut-off to limit the uniaxial compression (σ 2 = σ I − σII, see Fig. 2). In terms of the stress invariants σ I and σ II , this can be written as: σ I is the isotropic normal stress (defined as negative in compression), σ II is the maximum shear stress, c is the verticallyintegrated cohesion, µ (= sin φ) is the coefficient of internal friction of ice, φ is the angle of internal friction and σ c is the 155 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 σ c0 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.

160
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 Ψ (0 < Ψ < 1) is used to scale the stress back on the yield curve. It is defined as (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 165 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 time scale T d (Dansereau et al.,170 2016): where T d is set to the advective time scale associated with the propagation of elastic waves in undamaged ice (i.e., T d = ∆x/c 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 break up of ice bridges at small time scales. For the same reason, the advection of damage is neglected. The relaxation time scale (T d /(Ψ − 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 time scale 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 180 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 CFL criterion associated to the elastic waves.
The Elastic stiffness E and Maxwell viscosity η are written as a non-linear function of d, with a dependency on the ice 185 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 190 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:

195
where λ 0 (= η 0 /Y = 10 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 time scale 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,

200
the non-linear 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 205 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 non-linear relationship of the viscous relaxation time scale 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 210 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.

Numerical approaches 215
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 independently from 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 scheme (Rampal et al., 2016). In the standard VP numerical framework, 220 the stress components do not appear explicitly in the momentum equation. Instead they are written in terms of the non-linear 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.

Time discretization 225
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: where n − 1 is the previous time level and where: Note that σ 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 : 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, that is : τ n y = ∂ γ n σ n−1 yy ∂y + ∂ γ n σ n−1 xy ∂x + τ n ay + τ n wy .
The MEB rheology equations can then be implemented in a VP model by setting the VP bulk and shear viscosity to ζ V P = 240 ξ C1+C2 2 and η V P = ξC 3 respectively, setting the pressure term P = 0 and adding the stress memory terms.
The variable E n and λ n in Eq. 18 to 21 are discretized explicitely, as:

Space discretization
The model equations are discretized in space using a centered finite different scheme on an Arakawa C-grid. In this grid, the diagonal terms of the σ and˙ tensors are naturally computed at the cell centers and the off-diagonal terms at the grid nodes.
The x-component of the momentum equation are written as : where : The shear terms in Eq. 32 and 36 (˙ 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 neighbouring cell centres, 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 Eq. 15, 17 and 19.
The shear stress at the cell centre 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 270 checker board 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 (ξ n z˙ n xy ), which are integrated in another shear stress memory term, defined at the grid center. That is: where σ n xy i,j | C is the uncorrected shear stress at the grid center, ξ n z˙ n xy i,j is the shear stress increment averaged as in Eq. 37 275 and σ n−1 xy i,j | 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).

Numerical solution
With nx tracer points in the x-direction and ny in the y-direction, the spatial discretization on our C-grid leads to a system of 280 N = (ny(nx+1)+nx(ny+1)) non-linear 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 non-linear 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 non-linear system of equations is linearized and solved using a preconditioned 285 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 res = 10 −10 N m −2 . The uncorrected stresses σ 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. 290 1. Start with initial iterate u 0 do k = 1, k max 2. Linearize the non-linear system of equations by using u n,k−1 , h n,k−1 , A n,k−1 and d n,k−1 3. Calculate u n,k by solving the linear system of equations with FGMRES 7. If the Picard solver converged to a residual < res , stop. enddo 8. Update the stress memory σ n = Ψ n σ n 300 where 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.

Results
In the following, we present a series of idealized simulations to document the formation and break-up of ice arches with 305 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 x 200 km with a spatial resolution of 2 km (Fig. 3). The boundary conditions are periodic on the state at all time, which allows us to determine the critical forcing associated with a fracture event.

Control run
The break up 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 Fig. 5 and 6), separated by periods of elastic stress build up (low brittle failure activity). Two major fracture events are seen in the simulation (stage B and D in Fig. 4). The first corresponds to the failure of ice in tension with 320 the development of an ice arch on the downstream side of the channel (Fig. 5). The damage occurs on very short time scales (within minutes), and preconditions the formation of an arching flaw lead downstream of the ice bridge over longer time scales (Fig. 5b), in accord with results from Dansereau et al. (2017). The second event corresponds to the collapse of the landfast ice bridge with the break up of ice within and upstream of the channel (Fig. 6). As for the downstream ice arch, the lines of fractures are formed on short time scales and precondition the location of ridging on the advection time scale (Fig. 6b). The 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 340 channel entrance.

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 345 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 sub-critical. To first order, the arching progression of the fracture from the channel corners follows the second 7c). This differs from the expected angle of fracture in a coulombic material of θ = ±(π/4 − φ/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 non-zero 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, 355 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 changes from stage A, except for their increase in magnitude due to higher forcing ( Fig. 9). As the yield parameters (c, σ c ) are not 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.  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 localised (see Fig. 6b), but slowly expands in the direction perpendicular to the lines of fracture. This follows from the increase in material strength with 375 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 localised in time (see peaks in stage E, in Fig. 4).

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 380 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 (σ t = c/µ). 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.

Cohesion
Changing only the cohesion c 0 (with a fixed internal angle of friction φ) moves the entire yield curve along the first stress 385 invariant (σ I ) axis. For example, a higher cohesion increases the isotropic tensile strength σ t0 = c 0 / sin φ 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 section 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 break up can be related to the cohesion using the 1D steady state 390 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 τ LF Ic can be related to the cohesion as: In the simulations, the critical forcing for the complete decohesion of ice bridges (point 5 in Fig. 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.

400
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 exceeds 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 405 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 (40kN 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.

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

Tensile strength
The yield curve modifications discussed above (varying c 0 and φ) also change the tensile strength (both uniaxial and isotropic)

430
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 ν = 1/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 (5 < c < 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 a tens of kilometers unless anchor points are provided by stamukhi (Mahoney et al., 2014).

Compressive strength criterion
Not used in other MEB implementations (Dansereau et al., 2016(Dansereau et al., , 2017, the compressive cut-off offers a limit on the simulated uniaxial compression, which can reach unrealistically large values and cause numerical instabilities (see section 4). Including a compressive strength criterion (σ I −σ II > σ c ) can modify the upstream fracture event (stage D) by the development of uniaxial compression fractures along the upstream coast of the islands, if the uniaxial compressive stress upstream of the islands exceeds 445 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: 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 450 ideal case, the compression strength criterion is: The compression criterion can thus be written as a function of the surface forcing, as: Whether the ice will fail in shear (Mohr-Coulomb criterion) or in compression can be evaluated by substituting τ LF I from Eq.

455
(39) into Eq. 46, yielding the 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 , compression fracture occurs before the ice bridge break up and a ridge forms along the upstream coastlines, propagating in the channel 460 entrance while the ice in the channel remains landfast (Fig. 13b). If the terms are of 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).

Discussion
In the Arctic, ice arches are commonly observed upstream of narrow channels, where granular floes jam when forced into the 465 narrowing passage. This requires the ice 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 470 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 θ(= π/2 − φ/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 475 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 . This raises the question 480 whether the lines of fracture may be influenced by the stress correction path used in the damage parameterization, which determines the stress state associated to 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 485 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 loose their longitudinal symmetry 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 non-linear solver converges rapidly, within 6 iterations, given the small time step required by the CFL criterion to resolve the elastic waves. The errors are rather related 490 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-step, 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 495 stress state has a relative error , the error on the corrected stress is (see derivation in Appendix C): where If σ I > 0 (tensile stress state), 0 < R < 1 (triangle inequality) and the error on the memory terms ( M ) is of the same order 500 as that of the uncorrected stress state ( ≤ M ≤ √ 2 ). If σ I < 0 (compressive stress state), we have R ≥ 1, and 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 time scale, and this issue could be mitigated by decreasing the viscous coefficients η 0 . Using a compressive strength cut-off capping also offers a limit to the 505 uniaxial compression and reduces this instability. Another solution could be using a non-linear 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.

Conclusions
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 515 rheology is implemented in the same framework of 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 time scales. 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 520 yield curve influences the localisation 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  showing that the fracture orientation is determined by the planes of damage instability. Preliminary results suggests 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 525 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 to the material heterogeneity used to initialise the model (Dansereau et al., 2016), such that 530 this instability difficult to detect. A possible solution to this problem would be to use a non-linear yield curve which 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.
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 in 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 are in the range of 5-10 kN m −2 , lower 540 than stress buoys which measure both dynamical and thermal stresses at smaller scales but in the range of values previously associated to ice arch observations.
Code availability. Our sea-ice model code and outputs are available upon request.
Appendix A: Damage factor Ψ Let σ I and σ II be the stress invariant at time level n before the correction is applied, and σ If and σ IIf the corrected stress 545 invariant lying on the yield curve. Following Bouillon and Rampal (2015) we use a damage factor Ψ (0 < Ψ < 1) to reduce the elastic stiffness and bring the stress state onto the yield curve. I.e. : Substituting these relations into the Morh Coulomb criterion (σ IIf + µσ If = c) we solve for Ψ: Note that this relation implies that the stress correction is done following a line from the stress state (σ I ,σ II ) to the origin (see Fig. 2). This scheme stems from applying the damage factor to each individual stress components. 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).

555
Appendix B: Analytical solutions of the 1D momentum equation Considering an infinite channel of landfast ice (u = 0) along the y-direction with forcing τ LF Y = τ y , we write the 1D steady state momentum equation as: where we have neglected the ∂/∂y terms. In this case, the normal stress is zero in the entire channel and the stress invariants 560 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 = ± W 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 τ LF I = τ y , by integrating the 1D momemtum equation in which the ∂/∂x terms are neglected. That is: 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 land-fast ice and the stress invariants are function of both σ xx and σ yy : This allows to write the Mohr-Coulomb criterion in terms of σ yy : σ II + sin φσ I = 1 + 2 sin φ 3 σ yy < c,

Appendix C: Error propagation analysis
The error δF associated with a function F (X, Y, Z, ...) with uncertainties (δx, δy, δz, ...) 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 (σ ij ) and the damage factor Ψ (Eq. 13): Using Eq. A2, this can be rewritten in terms of the uncorrected stress invariants (σ I , σ II ): Assuming that the model has converged to a solution within an error on the stresses δσ ij = σ ij , δσ where is a small number, the model convergence error propagates on the stress memory with an error of : Substituting (δσ ij , δσ I , δσ II ) 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 of the stress invariants as : is grateful for 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. We acknowledge the use of imagery from the NASA Worldview application (https://worldview.earthdata.nasa.gov), part 610 of the NASA Earth Observing System Data and Information System (EOSDIS).     tensile strength (σt) and uniaxial tensile strength (σ * I , where the second principal stress invariant σ2 is zero, or σI = σII = σ * I ). The stress before and after the correction (see Eq. 13) is denoted by σ , and σ f respectively. The correction from σ to σ f is done following a line going through the origin.