MEB constitutive equation

The MEB rheology consists of a linear elastic part of the constitutive equation for a continuous solid, a viscous part of the constitutive equation for irreversible deformations, a local Mohr–Coulomb (MC) criterion for brittle failure, and an isotropic progressive damage mechanism that rescales the viscous and elastic dynamics to initiate avalanches of damage (Dansereau and others, Reference Dansereau, Weiss, Saramito and Lattes2016). We repeat the main aspects here for clarity.

The constitutive equation for vertically integrated internal stress ${\boldsymbol \sigma }$ (here in $\hbox{Pa}\, \hbox{m} = \hbox{N}\, \hbox{m}^{-1}$) and strain rates $\dot {\boldsymbol \varepsilon}$ for a 2-D compressible, visco-elastic, continuous solid is

(1)$$\dot{\boldsymbol \sigma} + \lambda^{-1}{\boldsymbol \sigma} = E( d) \, {\bf C}\cdot\dot{{\boldsymbol \varepsilon}}$$

with the elastic modulus tensor ${\bf C}$ (a function of the Poisson ratio ν) and the viscous relaxation timescale λ. Note that the stress and strain rate tensors are reduced in their order by using the Voigt notation for symmetric tensors. The relaxation timescale λ is written as the ratio of the viscosity ξ, the elastic modulus E and a damage parameter d representing the amount material degradation from accumulating micro (sub-grid) cracks in the sea ice (Dansereau and others, Reference Dansereau, Weiss, Saramito and Lattes2016):

(2)$$\lambda = {\xi( d) \over E( d) } = \lambda_0 \, ( 1-d) ^{\alpha-1}$$

where ξ and E depend on the fractional ice cover, mean ice thickness (i.e. using the formulation for the VP ice strength of Hibler, Reference Hibler III1979) and α > 1 is a parameter ruling the transition from elastic to viscous behaviour. E ₀ and ξ ₀ are the undamaged mechanical parameters and λ ₀ = ξ ₀/E ₀. In contrast to some previous work (e.g. Dansereau and others, Reference Dansereau, Weiss, Saramito and Lattes2016), we define damage so that d = 0 for undamaged ice and d = 1 for maximally damaged ice.

The damage increases when the stress states exceed the yield curve (Fig. 1) and it contains the history of the previous damaging events (Dansereau and others, Reference Dansereau, Weiss, Saramito and Lattes2016; Plante and others, Reference Plante, Tremblay, Losch and Lemieux2020). The increase depends on the scaling factor d _crit (critical damage, which is determined by the requirement to bring the overshooting stress state back to the yield curve).

Figure 1. Illustration of elliptic yield curve (VP, black dotted and solid ellipses) and MC yield curve (MEB, black piecewise linear lines). Invariant stress axes (σ _I, σ _II) in black and principal stress axes (σ ₁, σ ₂) in grey. σ _c is the critical uniaxial compressive stress (Eqn (3)) and σ _t is the critical tensile stress (Eqn (4)). The maximum tensile stress T _m (Eqn (6)) is indicated by the green dashed line. a and b denote the semi-major axes of the elliptic yield curve. Grey shading marks the cohesive stress states.

One possible yield curve for the MEB rheology is the MC criterion with a tensile cut-off (Fig. 1) (Dansereau and others, Reference Dansereau, Weiss, Saramito and Lattes2016). The critical uniaxial compressive stress σ _c at the intersection of the MC yield curve with the principal stress σ ₁ axis (Fig. 1) is

(3)$$\sigma_{\rm c} = 2\, ch\sqrt{q}$$

where c is the cohesion and q = ((μ ² + 1)^1/2 + μ)² is the slope defined by the internal friction coefficient μ. In contrast to the standard elliptic yield curve of a VP rheology (Hibler, Reference Hibler III1979), this yield curve permits isotropic tensile stresses. The critical tensile stress σ _t is defined as the intersection of the principal stress σ ₂ axis with the MC criterion (Fig. 1) so that

(4)$$\sigma_{\rm t} = -{\sigma_{\rm c}\over q}.$$

Implementation details

The finite-volume implementation on the C-grid of the MITgcm sea-ice model follows for the most part the implementation of the MEB rheology in the finite-differences C-grid implementation in the McGill sea-ice model (Plante and others, Reference Plante, Tremblay, Losch and Lemieux2020).

We note the structural similarity of the MEB and VP constitutive equations: the product of the elastic modulus tensor ${\it C}$ and the strain rate tensor $\dot {{\bf \varepsilon }}$ in Eqn (1) is (e.g. Dansereau and others, Reference Dansereau, Weiss, Saramito and Lattes2016)

(5)$$\left[{\boldsymbol C}\cdot\dot{{\boldsymbol \varepsilon}}^n\right]_{ij} \, = \, {\nu\over ( 1 + \nu) ( 1-\nu) }\dot{\varepsilon}_{kk}\delta_{ij} + {1\over 1 + \nu} \dot{\varepsilon}_{ij}.$$

This is the same form as the VP constitutive equation $\sigma _{ij} = 2\eta \dot {\varepsilon }_{ij} + [ ( \zeta -\eta ) \dot {\varepsilon }_{kk} -{P}/{2}] \delta _{ij}$ with P = 0 and shear and bulk viscosities η = (1/2)(1/(1 + ν)) and ζ = (1/2)(1/(1 − ν)). After re-interpreting these variables, we can re-use most of the VP code without additional changes. For details of the discretization we refer to Plante and others (Reference Plante, Tremblay, Losch and Lemieux2020).

On the staggered C-grid, some variables are naturally defined at centre (C) points (e.g. σ ₁₁), while others are naturally defined at corner (Z) points (e.g. σ ₁₂ and $\dot {\varepsilon }_{12}$) (Losch and others, Reference Losch, Menemenlis, Campin, Heimbach and Hill2010). Numerical stability requires that σ ₁₂, d _crit, d, λ ⁻¹ and E are defined on both C- and Z-points of the C-grid cell. The associated averaging is reduced to a minimum, so that only d _crit, d, h, a are linearly averaged to Z-points and only $( \dot {\varepsilon }_{12}) ^2$ is averaged to C-points. E, λ ⁻¹ and σ ₁₂ are computed for centre and corner points with the averaged variables.

Validation

We confirmed the plausibility of our MEB implementation with analytic solutions and symmetry tests (not shown, see chs 6 and 7 in Bourgett, Reference Bourgett2022) and with a reproduction of an idealized ice channel (Plante and others, Reference Plante, Tremblay, Losch and Lemieux2020, not shown).

The general behaviour of the dynamics is identical to previous results (Plante and others, Reference Plante, Tremblay, Losch and Lemieux2020). At the beginning of the simulation the tensile stresses downstream of the channel increase and damage develops downstream of each channel boundary. After 3300 s (τ = 0.06 N m⁻²) a concave shape at the downstream end of the channel indicates the ice arching effect. The stress values agree with the previous results (Plante and others, Reference Plante, Tremblay, Losch and Lemieux2020, their fig. 9). The divergent stress in the middle of the channel is small. The tensile stresses and the shear stresses in the downstream corners of the channel increase so that damage extends over the channel . The ice detaches from the upstream coastline but does not move yet (it remains landfast, land-locked by the islands). Both shear and divergent stress fields downstream of the ice channel drop to zero when the ice downstream of the channel detaches (not shown, see ch. 7 in Bourgett, Reference Bourgett2022).

Channel with idealized ice bridge

Inspired by previous ice arch simulations (Dumont and others, Reference Dumont, Gratton and Arbetter2009; Dansereau and others, Reference Dansereau, Weiss, Saramito, Lattes and Coche2017), we use an idealized channel set-up modified from Plante and others (Reference Plante, Tremblay, Losch and Lemieux2020). A 800 km × 200 km domain with a grid resolution of 2 km and closed boundaries with a no-slip boundary condition in the x -direction features a channel in the y -direction. The channel itself is 200 km long and 60 km wide. The domain has open boundaries at y = 0 km and y = 800 km with Neumann conditions for all variables. The Neumann conditions ensure that sea ice can drift freely into and out of the domain and does not need to detach from a solid boundary at y = 800 km, so that slowing down of the ice upstream the channel is solely determined by the ice arching. The sea-ice cover is forced by surface stress in the negative y -direction (‘southwards’) that increases linearly from 0 to 0.625 N m⁻² within 10 h. The simulation is run for 240 h with no further increase of the forcing. The slowing down of the sea ice upstream of the channel due to the formation of ice arches is used for comparison between VP and MEB.

Different parameters of the yield curves were tested to allow cohesive stress states. We choose the parameters so that the maximum tensile stress T _m (Eqn (6)) of the VP rheology is equal to the critical tensile stress σ _t (Eqn (4)) of the MEB rheology, since both represent the maximum positive value of the principal stress σ ₁ (Fig. 1). Specifically, we choose c = 10 and 30 kN m⁻² leading to σ _t = 10.4 and 31.24 kN m⁻¹ for MEB. The corresponding T _m are computed with P* = 49.92 kN m⁻² and P* = 149.92 kN m⁻², k = 0.05 and a small value for e = 1.2 (Kubat and others, Reference Kubat, Sayed, Savage and Carrieres2006; Lemieux and others, Reference Lemieux, Dupont, Blain, Roy, Smith and Flato2016).

Except for the VP simulations with T _m = 31.24 kN m⁻¹, the effect of ice arching to the upstream ice drift velocities can be observed and the ice slows down for both VP and MEB simulation (Fig. 2). The parameter set with T _m = 31.24 kN m⁻¹ makes the ice so stiff that it does not start to move at all. The ice drift in the MEB simulation with c = 30 kN m⁻² decreases within 40 h. For 10 kN m⁻², the ice drift upstream increases quickly and then slows down gradually with rates that are very similar between the MEB (m _meb = 1.45 × 10⁻⁷ m s⁻²) and VP simulations (m _vp = 1.43 × 10⁻⁷ m s⁻²) (Fig. 2, solid lines). Also, the velocity fields upstream (Figs 2, 3) are very similar with c = 10 kN m⁻² for MEB and its correspondent mechanical parameters for VP. The maximum ‘southwards’ velocity upstream is reached after approximately 15 h.

Figure 2. Averaged ice velocities parallel to channel upstream of channel. The sea ice does not move at all (VP) or rapidly stops (MEB) for the high cohesion case (c =30 kN m⁻², P* =149.92 kN m⁻², dash-dotted lines). There is a slow and very similar stopping effect by the formation of an ice arch in both the MEB simulation and the VP simulation for the low cohesion case (c =10 kN m⁻², P* =49.92 kN m⁻², dashed lines). The solid lines are the linear regression of the ice velocities.

Figure 3. Snapshots of the effective ice thickness h and the ice drift velocity (arrows) for the VP rheology (left two panels) and the MEB rheology (right two panels) at t = 12 and 24 h. Note that the colour scale is chosen to emphasize deviations from the initial state (h =1 m).

The effective ice thickness is generally similar for both rheologies (Fig. 3). In both cases, leads form downstream of the channel and ridging occurs upstream of the channel. Some differences in the exact location and shape of the leads and ridges are attributed to the different failure processes, namely the damage propagation and the associated normal flow rule for the MEB and VP rheologies, respectively. For instance, some ice remains attached to the islands downstream of the channel in the VP simulation as the deformation transitions from lead opening downstream of the islands to pure shear on the sides, while in the MEB simulation the damage propagation is directly along the coastlines. Upstream of the channel, the ridging area contains additional diagonal patterns in the MEB simulations due to the formation of secondary fracture lines, while the ice thickness is smoother and more uniform in the VP simulations.

Our results agree with other ice arch simulations (Dumont and others, Reference Dumont, Gratton and Arbetter2009; Dansereau and others, Reference Dansereau, Weiss, Saramito, Lattes and Coche2017; Plante, Reference Plante2021, ch. 5) and demonstrate that the cohesive strength of the ice plays an important role in ice arching so that corresponding mechanical parameters lead to similar results between the different rheologies.

Quadratic domain with cyclonic winds

A quadratic box with closed boundaries, constant anticyclonic (clockwise) ocean circulation and a moving cyclonic wind system was suggested to compare different sea-ice models (Mehlmann and others, Reference Mehlmann, Danilov, Losch, Lemieux, Hutter, Richter, Blain, Hunke and Korn2021). This ‘benchmark’ problem was used to analyse how different VP models simulate sea-ice deformation, in particular LKFs. Here, the ‘benchmark’ problem is repeated with the MITgcm using different grid spacings (Δx = 2,  4, 8 km) to analyse spatial heterogeneity in both the MEB and VP models. Note that for all grid resolutions the simulation is produced with the same time step (Δt = 120 s for VP, Δt = 0.5 s for MEB). In the MEB case, this value is chosen to ensure that the constitutive equation is well resolved at the highest resolution (i.e. according to the CFL criterion for resolving the elastic waves).

To choose similar yield curve parameters for MEB and VP as in the channel experiment, we have to consider the following: the VP parameters of this benchmark P* = 27.5 kN/m² with e = 2 and no tensile stress (k = 0) lead to a very low cohesion of 1.56 kN m² (Eqn (6)). Using the large P* values implied by the cohesion of the channel experiment in the VP rheology would change the benchmark dramatically from previously published results (Mehlmann and others, Reference Mehlmann, Danilov, Losch, Lemieux, Hutter, Richter, Blain, Hunke and Korn2021), so that to compare the different rheologies we instead adjust the MEB parameters to match the low VP cohesion. The low cohesion of 1.56 kN m² leads to very low stress states. For comparison, we also use a high value of c = 25 kN m². These cohesion values cover the range of previously reported values (Dansereau and others, Reference Dansereau, Weiss, Saramito and Lattes2016; Plante and others, Reference Plante, Tremblay, Losch and Lemieux2020). The model parameters for the experiments are summarized in Table 1.

Results from Mehlmann and others (Reference Mehlmann, Danilov, Losch, Lemieux, Hutter, Richter, Blain, Hunke and Korn2021) are reproduced by our VP simulations, with more radial features in the compressive stress field than circular ones and without tensile stress states (Fig. 4). In both models, there are fewer identifiable deformation patterns and the deformation fields also become smoother with decreasing resolution (not shown, see ch. 8 in Bourgett, Reference Bourgett2022).

Figure 4. Snapshot of the stress invariant σ _I at t = 2 d and with Δx = 2 km of the VP simulation on the left and the MEB rheology with low (centre) and high (right) cohesion. Positive values mean convergence. Divergent (negative) stress state are only allowed in the MEB model. The size of the stress invariant depends on the choice of the cohesion.

The presence of radial or circular features and the range of the stress values depends to some degree on the choice of cohesion for the MEB rheology. Using the MEB rheology with a cohesion similar to that of the VP simulations yields much smaller stresses but otherwise similar features as in the VP simulations, with mostly radial features and only a few circular stress features. Increasing the cohesion in the MEB model to get similar stress states as in the VP simulation (Fig. 4), however, changes these patterns and the features are mostly circular. Using the VP rheology with analogous values to match the cohesive stress states results in the same dependency (not shown). This dependency suggests that the shape of the features are sensitive to the shear strength; fewer cohesive stress states (grey area in Fig. 1) strongly result in smaller shear stresses, which favours radial features, while more cohesive stress states result in larger shear stresses, which favours the production of circular features.

There are other yield-curve related causes for the different stress fields: for example, the different yield curve shape of the MEB rheology allows isotropic tensile stresses (see negative σ _I in Fig. 4), while the VP rheology with the standard elliptical yield curve does not; and the MEB model does not have a flow rule. Neither of these causes are explored here.

Further, the rheologies are compared by means of the number of LKFs as detected by a tracking algorithm (Hutter and Losch, Reference Hutter and Losch2020, with parameter modifications by Mehlmann and others, Reference Mehlmann, Danilov, Losch, Lemieux, Hutter, Richter, Blain, Hunke and Korn2021) (Table 2). The number of LKFs increases for both rheologies with increasing resolution. As expected because of the damage mechanism and long-range elastic interactions that produce sub-grid fracturing (Dansereau and others, Reference Dansereau, Weiss, Saramito and Lattes2016), the MEB simulation (independent of the choice of cohesion) has more LKFs than the VP simulation on all grids, especially on Δx = 2 km grid. Increasing the cohesion tends to lead to fewer LKFs (Table 2). Note that the decreased heterogeneity in the MEB simulations with high cohesion is associated with a much less extensive damage field. As the damage mechanism is known to be a numerical error integrator (Plante and Tremblay, Reference Plante and Tremblay2021), this raises a question about the impact of numerical noise in seeding the heterogeneity.

Table 2. Number of LKFs for both VP and MEB rheology for simulations with 2 km, 4 km, and 8 km grid spacing Δx

The index ‘st’ indicates that the simulation uses a stochastic parameterization of c (cohesion) for the MEB rheology or P* (ice strength) for the VP rheology. The MEB simulations are run with a high value of c = 25 kN m⁻² and with a low value of c = 1.56 kN m⁻².