The material properties of ice bridges in the Maxwell Elasto-Brittle rheology

. The shape and break-up of landfast ice arches in narrow channels depend on the material properties of the sea-ice. The effect of the material parameters on ice arches in a sea ice model with the Maxwell Elasto-Brittle (MEB) rheology is investigated. The MEB rheology, which includes a damage parameterization, is implemented using the numerical framework of a Viscous-Plastic model. This conﬁguration allows to study their different physics independently of their numerical implementation. 5 Idealized ice bridge simulations show that the elastic part of the model together with the damage parameterization allows the propagation of fractures in space at very short time-scales. The fractures orientation is sensitive to the chosen angle of internal friction, but deviates from theory. It is speculated that these deviations stem from the absence of a ﬂow rule in the rheology. Downwind of a channel, the MEB model easily forms ice arches and sustains an ice bridge. Using a material cohesion in the range of 15-21 kPa is most consistent with the ice bridges commonly observed in the Arctic. Upstream of the channel, the 10 formation of ice arches is complicated by the absence of a relationship between the ice strength and the ice conditions, and by the presence of numerical errors associated with the damage parameterization. Results suggest that the formation of ice arches upwind of a channel is highly dependent on the rheology and calls for more analysis to determine the necessary conditions for their formation. In this paper, we investigate the role of the damage parameterization and the material strength parameters in the formation of ice arches using a MEB rheology that we implemented in an Eulerian ﬁnite difference VP model. To our knowledge, it 75 is the ﬁrst time the MEB rheology is implemented in a ﬁnite difference framework, such that the different physics can be assessed independently from the numerical implementation. Using idealised channel simulations, we show that the stress concentration, and thus memory, is dominant in determining the localisation of the ice fractures. We show that a rheology with a dependency between the ice strength and the ice thickness is needed for realistic ridge building. We also show that the simple stress correction used in the damage parameterization corresponds to a ﬂow rule that is independent of the orientation of the 80 failure surface. The stress correction scheme also ampliﬁes numerical convergence errors by orders of magnitude for large compressive states, introducing numerical artifacts that accumulate in the simulated ﬁelds.


Momentum and continuity equations
The simplified momentum equation for the 2D motion of sea ice used in this model is : 90 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, τ a (= τ axî +τ ayĵ ) is the air-ice surface stress and τ w (= τ wxî +τ wyĵ ) is the ice-ocean stress. The advection of momentum is neglected, being orders of magnitude smaller than the other terms in the momentum balance (Zhang and Hibler, 1997;Hunke and Dukowicz, 1997). The ocean current, and thus the sea surface tilt term, is set to zero and the Coriolis term is ignored, as it is identically zero for landfast ice. These simplifications are reasonable for landfast ice and keep the axial 95 symmetry of the problem. The air-ice stress τ a and ice-water surface stress τ w terms are defined using standard bulk formula (McPhee, 1979), with the air and ocean turning angles set to zero and assuming that the wind velocity is orders of magnitude larger that of sea ice: where C da and C dw are the air and water drag coefficients, u a is the surface wind and ρ a and ρ w are the air and water density (see values in Table 2).
The continuity equations used for the temporal evolution of mean ice thickness h (volume per grid cell area) and concentration A are written as: where S h and S A are thermodynamic sink and source terms for ice thickness and compactness respectively. For simplicity, S h and S A are set to zero in this study.
where E is the Elastic Stiffness defined as the vertically integrated Young Modulus of sea ice, λ is the viscous time relaxation (λ = η E , η being the viscosity) and˙ is the strain rate tensor. The advection of stress is neglected in this study as we focus on the landfast ice. In the linear elastic regime (when the main balance is between the first term on the lhs and the rhs of Eq. 6), deformations are relatively small, occur over a small time-scale and are reversible. As such, the elastic strains return to zero if 120 all loads are removed and the model includes a material memory of the strain history. In the viscous-elastic regime (when all three terms are important), viscous deformations, or creep, relax the elastic stresses over a longer time-scale if the forcing is sustained, which results in permanent deformations. The brittle regime occurs when the internal stresses reach a material yield criterion. In this case, brittle fractures (or cracks) progressively reduce the elastic stiffness and visco-elastic ratio, allowing for larger visco-elastic deformations with the same amount of stress. This brittle component is thus comparable to the plastic 125 regime of the standard VP model of Hibler (1979), with deformations that are relatively large, and stresses that are strain rate independent.
Note that while Eq. 6 is similar in form to the stress-strain relationship of the Elastic Viscous Plastic (EVP) model (Hunke and Dukowicz, 2002), the elastic component in the EVP model does not represent the true elastic behavior of sea ice. Rather, the elastic term was introduced to improve the model convergence and computational efficiency of the VP model, allowing for 130 an explicit numerical scheme and simple parallelization (Hunke and Dukowicz, 1997). In fact, in the new implementation of the EVP model (Hunke, 2001), elastic waves are faster than the observed elastic waves unless a time-step of several hours or more is used (Williams et al., 2017). In the MEB model, the elastic component represents the true 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 135 (Tucker and Perovich, 1992;Sukhorukov, 1996;Hata and Tremblay, 2015b) in field experiments. The viscous relaxation term is also analogous to the viscous term in the thermal stress models of Lewis (1993) and Hata and Tremblay (2015a).

Linear Elasticity
The elastic component of the MEB rheology implies that the internal stresses are related to the strains (units of m/m) rather than to the strain rates (units of s −1 ). We first write the stress-strain relationship of a 2-D linear elastic solid, that is, ignoring 140 the viscous and brittle components of the rheology (Rice, 2010;Bouillon and Rampal, 2015)) : where E is the Elastic Stiffness of sea ice, σ ij are the vertically integrated stresses acting in the j th direction on a plane perpendicular to the i th direction, ν ( = 0.33) is the Poisson ratio and δ ij is the Kronecker delta. The linear elastic stress component of Eq. 6 is obtained by taking the time derivative of Eq (7), assuming negligible variations in E: where ":" denotes an inner double tensor product and where The components of the 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 2010). This form is also used in Bouillon and Rampal (2015); Rampal et al. (2015), but differs from Dansereau et al. (2015) who used the components derived from the plane strain approximation, in which case the vertical components of strain are zero but the vertical normal stress component (σ 33 ) is nonzero. This should be avoided, as it implies that the vertical stress is of significant order (Rice, 2010) and determined by the 155 horizontal normal stress: σ 33 = ν(σ 11 + σ 22 ).

Brittle fracture parameterization
In this model, the critical stress at which the ice fails is determined by the Mohr-Coulomb failure criterion. This yield criterion is based on laboratory experiments (Schulson et al., 2006) and was found to agree with field observations (Weiss et al., 2007).
A criterion is further used on the second principal stress (σ 2 = σ I − σII) to limit the compression (see Fig. 2). The yield 160 function can be written in terms of the stress invariant σ I and σ II , as: where σ I is the isotropic normal stress (compression defined as negative), σ II is the maximum shear stress, c is the cohesion, µ = sin φ is the coefficient of internal friction of ice, φ is the angle of internal friction and σ c is the compressive strength.
The cohesion, angle of internal friction and compressive strength are material properties derived from in-situ observations (see 165   Table 1 for values and references) and laboratory experiments (Timco and Weeks, 2010). The model values used in our control run are listed in Table 2 along with a summary of the value used in other studies in Table 3.
The fracture of ice is parameterized as a local decrease in elastic stiffness and consequently as a local increase in the magnitude of the elastic deformation associated with a given stress state. The local increase in deformations in the material results in the concentration of internal stresses in adjacent grid cells, leading to the propagation of the brittle fracture in space. The elastic stiffness is written in function of a damage parameter d representing the weakening of the ice upon fracturing (0 < d ≤ 1, Bouillon and Rampal 2015), with a dependency on the ice thickness and sea ice concentration inspired by the ice strength parameterization of Hibler (1979) : where Y (= 1 GPa) is the Young Modulus of undeformed sea ice. This value is smaller than that used in the neXtSIM model 175 (10 GPA, Bouillon and Rampal 2015) and similar to that of the MEB model (0.8 GPa, Dansereau et al. 2016). The damage parameter d has a value of 0 for undamaged sea ice and 1 for fully damaged ice, and represents the amount of fractures and material degradation present in the ice. Note that the yield parameters (c, µ, σ c ) in the MEB model are not a function of the damage parameter (nor of the ice thickness or concentration). The decrease in elastic stiffness when the ice fails yields larger deformations (and potentially changes in h and A) for the same state of stress, but does not change the critical stress (i.e., the ice 180 strength). This is in contrast with a strain-weakening or strain-hardening material, a property that has been observed in sea ice (Richter-Menge et al., 2002). The lack of strain hardening in the MEB model leads to non-physical results in convergence with the absence of ridge propagation in the direction parallel to the second principal strain (maximum axial compressive strain).
This will be discussed in section 3.1.4.
Following Rampal et al. (2015), the introduction of damage upon failure is proportional to the local stress in excess of the 185 yield criterion. A damage factor Ψ (0 < Ψ < 1) is used to keep the stress 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. In the following, prime quantities are used to indicate stress terms that are not corrected by the damage factor Ψ. Note that the stress 190 components are all scaled by the same damage factor, following a line from the uncorrected stress state to the origin (see Fig.   2). This correction does not corresponds to a flow rule: only the magnitude of the excess stress is used to increase the damage parameter. This determines the magnitude of the strain associated to a stress state (and indirectly, the strain rate tensor), but not its orientation. This is different from the normal flow rule of Hibler 1979 (see for instance Bouchat and Tremblay 2017), in which the yield curve surface determines the relative importance of shear and divergence in the deformations.

195
After each fracture event, the damage parameter is updated such that (1 − d f ) = Ψ(1 − d ). That is, the damage factors are accumulated in the damage parameter during the simulation, such that the damage gradually increases with each fracture. As the amount of stress overshoot is time-step dependent, the temporal evolution of the damage parameter is parameterized as a simple relaxation with a damage time scale T d (Dansereau et al., 2016):

200
where T d is set to the advective time scale associated with the propagation of elastic waves in undamaged ice where ∆x is the spatial resolution of the model and c e the elastic wave speed). Note that the relaxation time scale in Eq.
(14) (T d /(Ψ − 1)) is time-step dependent via its dependency on the damage factor Ψ: 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 lead Dansereau et al. (2016) to argue that the model time step should be 205 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 are 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 should be favored if possible in order to resolve the elastic waves.

210
Note that 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.

Maxwell viscosity
The viscosity η used to define the viscous time relaxation (λ = η E ) is parameterized to change faster with the damage parameter than the elastic stiffness (Dansereau et al., 2016) : with η 0 being the viscosity of undeformed sea ice and α an integer set to 4 that determines the smoothness of the transition from linear elastic behavior to viscous behavior. This definition ensures that the viscous term is negligible in undamaged ice, but important in heavily damaged ice (see Eq. 6, where η appears in the denominator). Note that as the damage increases asymptotically to 1, the singularity in Eq. 15 (if d = 1) never occurs. The relaxation time λ in Eq. 6 can then be written using 220 Eq. 12 and 15 : where λ 0 (= η 0 /Y = 10 5 ) is a parameter that corresponds to the viscous relaxation time scale in undamaged sea ice with 1m mean thickness. The value of about one day falls into the range of observations (see Table 1). Note that if λ 0 is sufficiently high, the MEB rheology reduces to the Elasto-Brittle rheology Rampal et al., 2015). The set 225 dependency of λ on the mean ice thickness and concentration ensures that both the total stress tends toward zero for low mean thickness or concentration (i.e. in free drift), while a continuous (A ∼ 1, h > 0) but heavily damaged ice cover behaves as a viscous material.

Numerical approaches
In this section, we discuss the numerical implementation of the MEB rheology exploiting the code framework of a standard VP 230 model. In this framework, the stress components in the momentum equation do not appear explicitly; instead they are written in terms of viscous coefficients and strain-rates. This is done in our implementation of the MEB model by treating the stress memory as an additional forcing. The damage parameterization is therefore the only new module that needs to be implemented.
n can then be written using Eq. 6, 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 17 is then substituted in the stress 240 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. 9) and where the stress memory terms have been included in 245 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 = 250 ξ C1+C2 2 and η V P = ξC 3 respectively, setting the pressure term to P = 0 and adding the stress memory terms.
The variable h n , A n and d n (and thus E n and λ n ) in Eq. 17 to 20 are discretized explicitely, as: where the damage factor is computed from Eq. 13 and Eq. 17.

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 265 (see Fig. 3). The x-component of the momentum equation are written as : where :

275
The shear terms in Eq. 28 and 32 (˙ xy , ξ z and γ z ) are thus defined at the lower-left grid node rather than at the grid center.
This 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 ) as a simple average of the neighbouring cell centres, i.e. : 280 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. 12, 16 and 18.
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 (as in Eq. 33 for the scalars) cause a checker board instability to appear in the solution, because of the staggered shear stress corrections and memories. To avoid this, we approximate the shear 285 stress at the cell center using a different shear stress memory term, defined at the grid center and only used in the damage parameterization, and only average the shear stress increments. That is: where σ n xy i,j | C is the uncorrected shear stress at the grid center and σ n−1 xy i,j | C is the corrected shear stress at the grid center from the previous time step. This approximation is only used to calculate the damage factor ψ.

290
Note that the averages in Eq. 33 and 34 cause a smoothing of the variables used to define the shear stress state. This may yield significant differences with the previous implementation of the damage parameterization, which was developed for models using a Lagrangian  or Finite Element Method (Dansereau et al., 2016) schemes, but would allow insightful comparisons with VP and EVP model simulations using similar discretization schemes.

Numerical solution 295
The discretized equations described above are solved simultaneously using an IMplicit-EXplicit (IMEX) approach (Lemieux et al., 2014). This process 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 (momentum) 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 until the L2norm of the solution residual falls below a set tolerance of res = 10 −10 N/m 2 . The uncorrected stresses 300 σ n is then scaled by the damage factor Ψ n and stored as the stress memory σ n for the following time step.
To summarize, the IMEX stepping scheme can be written as : 1. Start time level n with initial iterate u 0 do k = 1, k max 2. Linearize the momentum equation using u n,k−1 , h n,k−1 , A n,k−1 and d n,k−1 305 3. Calculate u k,n by solving Eq. 19 and 20 with GMRES 7. If the Picard solver converged to a residual < res , stop.

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. 5) and space (see Fig. 6 and 12), 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. 5). The first corresponds to the failure of ice in tension with the development of an ice arch on the downwind side of the channel (Fig. 6). The second event corresponds to the collapse of 330 the landfast ice bridge with the break up of ice within and upwind of the channel (Fig. 12). The three remaining periods during which few new brittle fractures occur correspond to an elastic landfast ice regime (stage A), a stable downwind ice arch regime (stage C), and a drift ice regime when ice flows within, downstream and upstream of the channel (stage E).

Elastic regime: stage A
In the first stage of the simulation, elastic stress builds up but remains inside the yield curve in the entire domain such that 335 there is no brittle failure activity (Fig. 5, 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 wind with respect to the coastlines: there are large tensile stresses on the downwind coastlines, compressive stresses on the upwind coastlines and shear stresses on the four corners of the channel (Fig. 8). At the vertical line of symmetry (away from channel openings,  Figure   9). 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. 8c) is along the x-direction downwind of the coastlines (where the ice is in uniaxial tension), and along the y-direction upwind of the coastlines (where the ice is in uniaxial compression). In the downwind end of channel, the second 345 principal stress alignment takes an arching shape, transitioning to a vertical alignment towards the upwind channel entrance.

Downstream ice arch: stage B
The formation of the downstream ice arch (Fig. 6) is initiated at a wind forcing of ∼ 3.5 m/s. The initial fractures are located at the downwind corners of the channel where the stress state reaches the critical shear strength for positive (tensile) normal stresses (Fig. 10). The fractures then propagate from these locations and form an arch (Fig. 6, point 3). The progression of the 350 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. This results in the propagation of the fractures in space through regions where the internal stress state was originally sub-critical. This process occurs on very short time scales (within minutes), and preconditions the formation of an arching flaw polynya over longer time scales (Fig. 6b).

355
To first order, the arching progression of the fracture from the channel corners follows the second principal stress direction (i.e. a failure in uni-axial tension on the plane perpendicular to the maximum tensile stress, see Fig. 8c). This differs from the expected angle of fracture in a coulombic material, at θ = ±(π/4−φ/2) from the second principal stress orientation (Ringeisen et al., 2019). This deviation results from the absence of a flow rule in the MEB model. That is, only the Elastic stiffness is changed by the damage parameter to scale back the uncorrected stress to the critical state. The strain rate tensor associated 360 with the fracture is hence determined by the change in stress state at the end of the non-linear solution, rather then by the yield curve surface (see Fig. 10). The flow rule discrepancy is discussed in more details in section 3.2.2.

Stable ice arch regime: Stage C
A second period of low brittle fracture activity follows the formation of the ice arch (period C in Fig. 5). In this stage, the ice downwind of the ice arch is detached from the land boundaries and starts to drift. Upstream of the ice arch, the elastic 365 stresses, except for their increase in magnitude due to higher wind forcing, show little changes from stage A. 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, the increasing wind forcing constantly increases the stress states beyond the yield criterion, leading to further damage. Note that unless the yield parameters depend explicitly on the damage, tensile fracturing does not reduce the critical stress, in contrast to real ice fractures. As such, large tensile and shear stresses persist along and north of the ice arch 370 after the ice arch is formed (Fig. 8b). The formation of a stress-free surface could be created by defining the cohesion as a function of the ice thickness and/or damage. The second break-up event (Stage D) corresponds to the fracture of ice upwind of the channel and the collapse of the ice bridge between and upstream of the islands (Fig. 12). This fracture starts at a wind speed of 5.67 m/s on the upwind corners of the 375 islands where the internal stress reaches the critical shear strength for negative (compressive) normal stresses (green point in Fig. 10). The propagation of damage from these points is composed of two separate fractures. First, a shear fracture progresses downwind along the channel walls (Fig. 12, point 5). This results in the decohesion of the landfast ice in the channel from the channel walls, increasing the load on the ice arch and in the landfast ice north of the channel. Second, a shear fracture propagates upwind from the channel corners at an angle 61.8 • from the coastline (Fig. 12, point 6). This orientation corresponds to an 380 angle θ = 28.2 • from the second principal stress orientation (Fig. 8c). Again, this angle deviates from the theoretical fracture orientation in a granular material with φ = 45 • , at 22.5 • from the second principal stress orientation (Ringeisen et al., 2019).
As for the downwind ice arch, the lines of fractures are formed at short time scales and precondition the location of ridging on the advection time scale (Fig. 12b).

Drift and ridge building: stage E 385
The last stage of the simulation (stage E) corresponds to a regime where most of the ice in the domain is drifting. As in Stage C, the non-zero brittle fracture activity corresponds to further damage being produced in the already damaged ice. In this stage, landfast ice only remains in two wedges of undeformed ice upwind from the islands (see Fig. 12b). 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. Note that the ridge building is highly localised, with no further stress build up elsewhere in the domain. This prevents 390 the formation of an ice arch upwind of the channel seen in observations (e.g. in the Lincoln Sea). This unrealistic behaviour of the model is another consequence of the use of constant strength parameters and yield criterion: instead of increasing the pressure with increasing ice thickness during ridging, the stress field is in a steady state set by the constant critical stress along the deformation lines (see Fig 13). The increasing wind forcing is then only balanced by the inertial term and water drag term, resulting virtually in a free drift mode where the ice velocity (and thus the rate of deformation at the ridging lines)  (Fig. 13c). This is left as future work.
Note that the damage field at the end of the simulations is highly sensitive to the solution residual tolerance res . With time, unless a very low res is used, the damage fields are no longer horizontally symmetrical about the center of the channel (Fig.   14). This indicates the presence of artifacts in the model, although the solutions always converged to the set precision. An error propagation analysis shows that these asymmetries are produced by 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 of , the error on the corrected stress is (see derivation in Appendix C): where

410
If σ I > 0 (tensile stress state), 0 < R < 1 (triangle inequality) and the error of the memory components ( M ) is of the same order 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, this issue can be improved 415 by decreasing the viscous coefficients η 0 . Another solution would be to use 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). A different stress correction path would furthermore allow the application of a flow rule based on granular physics. This would require a damage tensor that would scale the components of 420 the stress tensor independently, as commonly used in continuum damage mechanics. Implementing a damage tensor is left for future work.

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 425 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 a compression strength criterion on the simulated ice bridge. We place a particular focus on the propagation of the ice fractures in space both upwind and downwind of the channel.

430
Changing the cohesion c 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 σ t = c/ 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 wind forcing associated with the different stages of the simulations but does not change the series of events described in section 3 or the orientation of the ice fractures. This is in agreement with results from Dansereau et al. (2017).
The critical wind forcing associated with the ice bridge break up 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. 4). Using the yield criterion (Eq. 11) with σ I = 0 (i.e. σ II = c), the maximum sustain-440 able wind forcing τ ac can be related to the cohesion as: In the ice bridge simulations, the critical wind forcing for the break up of the ice bridge (stage D) follows the simple 1D model, although with lower wind forcing values (Fig. 15). This is expected considering the stress concentration occurring at the channel corners and the contribution of the ice upwind and downwind of the non-infinite channel, which pushes and pull 445 on the ice in the channel such that a smaller critical wind forcing is required to break the ice.
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 wind speed regularly exceeds 10 m/s (τ a > 0.15 N/m 2 ), this suggests a lower bound on the cohesion of sea ice of at least 15 kN/m (see green curve in Fig. 15). 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) suggests that the cohesion of sea ice should be 450 smaller than 21 kN/m (see red curve in Fig. 15). These values are lower than estimates from ice stress buoys measurements (40kPa, Weiss et al., 2007) which includes many driving forces neglected in our simulations, such as thermal stresses (Hata and Tremblay, 2015b) and, to a lesser extent, tides. Our values are similar to previous large scale estimates based on wind forcing alone (Tremblay and Hakakian, 2006) and to values used in the neXtSIM model (see Table 3).

455
The angle of internal friction φ, analogous to the static friction between two solids, determines the constant of proportionality (µ = sin φ) between the shear strength and the normal stress (see Eq. 11 and Fig. 2). In the following, we vary the angle of 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 orientation of the failure surface is determined by the point at which the Mohr circle reaches the yield criterion in the Mohr circle space. This point of failure is aligned at angle 470 θ = ±(π/4 − φ/2) from the first principal stress direction (Ringeisen et al., 2019). In the MEB model, the angle of fracture does not follow the theory (see Fig .16b). We speculate that the deviations are related to the absence of a flow rule linking the deformations to the yield curve and the angle of internal friction. Note however that the tendency of the fracture orientation remains consistent with the theory: decreasing the angle of internal friction increases the downwind ice arch curvature and the angle of the fractures (from the positive x-direction) upwind of the islands. This suggests that despite the absence of a flow 475 rule, the angle of internal friction determines the direction of crack propagation indirectly by changing the material strength.
Note that for angles of internal friction > 60 • , the upwind lines of fracture propagate away from the boundaries, and a second ice arch forms upwind of the channel (not shown). This angle of internal friction is un-realistically large: previous estimates of the coefficient of friction from ice in situ observations rather suggest values of 30 − 45 • (see Table 1). This difference stems from the fact that in the Arctic, the ice arches that are commonly observed upwind of a channel are formed when granular 480 floes jam when forced into a constricting channel in which the ice is not landfast. In our experiments, we rather simulate the propagation of ice fractures through the landfast ice upwind of a channel. It should be possible, however, to form these arches during the drift ice regime after the collapse of the ice bridges. In our simulations, two issues impede this process: the lack of compressive stress build up during ridge building and the flow rule (i.e. the path of the stress correction) which favors ridge building over sliding along the landfast edges, limiting the increase in compression directly upwind of the channel.

Tensile strength
The yield curve modifications discussed above (varying c and φ) also change the tensile strength (both uniaxial and isotropic) of ice. The tensile strength determines the magnitude of the critical wind forcing necessary for the formation of the downwind ice arch (stage B). Downwind from the islands, the tensile stresses can be approximated using the 1D version of the momentum equation as a function of the fetch distance F down (see Fig. 4) between the downwind coast of the islands and the bottom 490 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. 11) for the 1D case (Appendix B):

495
where ν = 1/3 was used. Substituting σ yy from Eq. 39 into Eq. 40, the yield criterion can be written in terms of the wind forcing and the material parameters: Using our cohesion estimates (15kP a < c < 21kP a) and angles of internal friction in the range of observations (30 and 45 • , this would suggest stable bands of landfast ice of extent F down ∼ 100 km should be sustainable. This is much larger than those 500 observed in the Arctic, where bands of landfast ice rarely exceed a few tens of kilometers unless anchor points are provided by stamukhi (Mahoney et al., 2014). This discrepancy highlights that neglected processes such as waves and tides are important in shaping the landfast ice cover, and may prevent such landfast extent from occuring.

Compressive strength criterion
Including a compressive strength criterion (σ I > σ c ) can modify the upwind fracture event (stage D) by the development 505 of compression fractures along the upwind coast of the islands. The compression strength only affects the simulation if the compressive stress upwind of the islands exceeds the compressive strength.
As for the downwind case, the critical wind stress for the development of a compressive fracture can be approximated using the 1D version of the momentum equation. The maximum normal stress at the upwind coast of the islands is: where F up is the distance between the top boundary of the domain and the upwind coasts of the islands (see Fig. 4). In the ideal case, the compression strength criterion is: The compression criterion can thus be written as a function of the wind forcing, as:

515
Whether the ice will fail in shear (Mohr-Coulomb criterion) or in compression can be evaluated by substituting τ a in Eq. (37) by Eq. 44, yielding the criterion: If this condition is met, the compression strength criterion does not influence the simulation, and the upwind shear fracture lines develop as in the control simulation (Fig. 17a). If the left hand side term of Eq. 45 is much smaller that σ c , compression 520 fracture occurs before the ice bridge break up and a ridge forms along the upwind coastlines, propagating in the channel entrance while the ice in the channel remains landfast (Fig. 17b). If the terms are of similar order, the decohesion of the ice bridge and the compression fractures are initiated simultaneously: compression fracture occurs along the upwind coastlines but not in the channel entrance, as the ice both in and upwind of the channel starts to drift (Fig. 17c).

525
The MEB rheology was implemented on the Eulerian, Finite Difference 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, this is the first time the MEB rheology is implemented on the same framework of a VP or EVP model. This will allow direct comparison of 530 these models using the same numerical platform in future work.
In idealised 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 cause a concentration of stresses close the preexisting damage. We also show that while the choice of yield curve influences the localisation and orientation of the ice fractures, the angles of fracture propagation differ from those expected in 535 a granular material such as sea ice. This is due to the simple stress correction scheme using a scalar damage parameter applied on all stress tensor components, such that the flow rule is determined by the stress state independently from the orientation at which the ice fracture is occuring. The angles of the fractures could be more physical by the use of a damage tensor, with a stress correction path derived from a physical flow rule. This is left as a future model development.
The stress correction scheme in the damage parameterization  is also found to cause a problematic 540 propagation of numerical errors in the stress memory terms. The magnitude of the error propagation 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. Note that these errors are hardly detectable when using material heterogeneity (Dansereau et al., 2016)  Based on these results, these are our recommendations for using the MEB model: -A very low solution residual tolerance res < 10 −10 should be use to limit the accumulation of errors associated to the correction scheme.

560
-The cohesion should be limited to c < 21 kPa.
the formation of unrealistically large ridges.
-An alternative stress correction scheme should be developed to limit the accumulation of errors in the stress memory terms.
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 the use of a single damage factor applied on the elastic stiffness, which is linear to each of the 575 stress components, i.e., the correction is applied equally to each 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 each component of the stress tensor (i.e. using a damage tensor, rather than a damage parameter). This could be used to cure the error propagation problem when large compressive stresses are present (see Appendix C).
Appendix B: Analytical solutions of the 1D momentum equation

580
Considering an infinite channel of landfast ice (u = 0) along the y-direction with wind forcing τ a = τ ay and water stress τ w = 0, 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 are σ I = 0, σ II = σ xy . The shear stress at any arbitrary point x across the channel can be determined by integrating Eq. B1 585 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 wind forcing τ = τ ay and water drag τ w = 0, by integrating the 1D momemtum equation in which the ∂/∂x terms are neglected. I.e. : ∂σ yy ∂y + τ ay = 0, 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 :

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