Simulating the processes controlling ice-shelf rift paths using damage mechanics

Abstract

Rifts are full-thickness fractures that propagate laterally across an ice shelf. They cause ice-shelf weakening and calving of tabular icebergs, and control the initial size of calved icebergs. Here, we present a joint inverse and forward computational modeling framework to capture rifting by combining the vertically integrated momentum balance and anisotropic continuum damage mechanics formulations. We incorporate rift–flank boundary processes to investigate how the rift path is influenced by the pressure on rift–flank walls from seawater, contact between flanks, and ice mélange that may also transmit stress between flanks. To illustrate the viability of the framework, we simulate the final 2 years of rift propagation associated with the calving of tabular iceberg A68 in 2017. We find that the rift path can change with varying ice mélange conditions and the extent of contact between rift flanks. Combinations of parameters associated with slower rift widening rates yield simulated rift paths that best match observations. Our modeling framework lays the foundation for robust simulation of rifting and tabular calving processes, which can enable future studies on ice-sheet–climate interactions, and the effects of ice-shelf buttressing on land ice flow.

1. Introduction

Ice-shelf rifting weakens ice shelves and precedes calving of tabular icebergs, which comprise the vast majority of calved Antarctic ice volume (Tournadre and others, 2016). Calving and submarine melting are the two major causes of the recent Antarctic losses of ice-shelf mass and buttressing of upstream grounded ice (Greene and others, 2022). Decreased buttressing can affect the discharge of grounded ice into the ocean and ice-sheet contribution to sea-level rise (Haseloff and Sergienko, 2022), and calved icebergs can transport freshwater to lower latitudes to influence ocean circulation and sea-ice growth (e.g. Jongma and others, 2009; Martin and Adcroft, 2010; Merino and others, 2016), as well as the marine biosphere (e.g. Arrigo and others, 2002; Laufkötter and others, 2018).

Rifts are typically oriented with depth within 5$^\circ$ of the vertical plane (Walker and Gardner, 2019), and propagate horizontally across ice shelves. The processes that control rifting are poorly understood, and it remains a challenge to capture rifting within computer simulations of ice-shelf evolution. Some rift activity has been correlated with environmental forcings, including surface temperature and tides (Olinger and others, 2019), and the arrival of tsunamis (Walker and others, 2013). However, past observational evidence and modeling suggest that rifting is primarily driven by viscous glaciological stresses associated with gravity-driven ice flow (e.g. Bassis and others, 2005, 2007, 2008; Joughin and MacAyeal, 2005; Borstad and others, 2016). In turn, these stresses are sensitive to the history of rift behavior (Wang and others, 2022), so that the rift state at one point in time directly feeds back to future rifting. Previous studies have also established that rift propagation is sensitive to crucial rift–flank boundary processes such as the ‘gluing’ together of rift flanks by mechanically coherent mélange (i.e. mélange capable of transmitting stress between flanks) (Larour and others, 2004), backpressure on rift–flank walls from ice mélange (Larour and others, 2014, 2021), and contact between flanks (Lipovsky, 2020). However, these studies only examined whether or not a rift would propagate based on the sharp fracture assumption, which is consistent with the discrete representation of a zero-thickness crack in the model geometry; but they did not model the propagation of rifts.

To capture rift propagation and its time-varying effects on ice-shelf stresses, we require advanced computational modeling approaches that account for the coupling between ice flow and fracture. Modeling the propagation of rifts and crevasses under the sharp fracture assumption, using the finite-element method (FEM) and linear elastic fracture mechanics (LEFM), introduces algorithmic complexities (Yu and others, 2017). While the displacement correlation method can simplify the implementation of LEFM with the FEM, it is still restricted to linear elastic ice rheology (Jimenez and Duddu, 2018). In contrast, continuum damage mechanics combined with the FEM exploits the diffuse fracture assumption, which is consistent with a smeared representation of the crack over the model mesh or grid. This diffuse approach obviates the need for complicated crack tracking algorithms and simplifies the incorporation of fracture processes within an ice flow model based on the non-linear Stokes equations. Recently, damage mechanics has been used to simulate glacier-scale crevasse propagation (Jiménez and others, 2017; Duddu and others, 2020; Sun and others, 2021; Clayton and others, 2022), and ice-shelf-scale mechanical weakening (Albrecht and Levermann, 2012, 2014; Borstad and others, 2016; Sun and others, 2017) and rift propagation (Huth and others, 2021b).

Here, we develop methods to incorporate rift–flank boundary processes within the computational framework based on anisotropic ‘creep damage’ and vertically integrated ice-shelf flow models (Huth and others, 2021a, 2021b). This quasi-three-dimensional (3-D) (ice-shelf viscous stresses are simulated in the horizontal plane and damage is evolved in a 3-D space) framework represents the initiation and time-dependent propagation of crevasses and rifts. We then investigate how mélange fill, mélange strength, and rift–flank contact influence the rift path through several parametric sensitivity studies. To demonstrate that the viability of the proposed framework, we simulate the observed final years of rifting on the Larsen C Ice Shelf that resulted in the calving of iceberg A68 in 2017 (Fig. 1). The paper is organized as follows: in Section 2, we describe the model equations and rift–flank boundary scheme; in Section 3, we present the parametric studies; in Section 4, we discuss the results; and in Section 5, we offer concluding remarks.

Figure 1. NASA MODIS images of Larsen C ice shelf on (a) 3 December 2014 and (b) 11 November 2017, 4 months after calving of iceberg A68. The blue star in (a) marks the initial tip position of the rift that propagated to calve iceberg A68. The yellow arrow in (a) indicates a damaged region, shown in detail in (c), that was not captured in the inversion (Fig. 5c). BIR, Bawden Ice Rise; GIR, Gipps Ice Rise; KP, Kenyon Peninsula.

2. Model equations

In this section, we summarize the ice-flow model, the anisotropic damage model, and the rift–flank boundary scheme. All model equations are presented in indicial notation: vectors are notated as ${\boldsymbol a} = a_i{{\hat e}}_i$, where i are the spatial indices of the Cartesian coordinate system (x ₁,  x ₂,  x ₃) = (x,  y,  z) and ${\hat e}_i$ are orthonormal basis vectors; second-order tensors are denoted as ${\boldsymbol A} = A_{ij}\ {\hat e}_i\otimes {\hat e}_j$, where ⊗ is the dyadic product of the Cartesian base vectors; and principal values of the tensor are denoted as 〈A _i〉. We adopt Einstein's convention where repeated spatial indices imply summation.

2.1 Ice flow model

We simulate ice flow with the two-dimensional (2-D) shallow shelf approximation (SSA), which is most appropriate for ice shelves and ice streams that have minimal or no basal drag, so that vertical shear is negligible (Macayeal, 1989; Huth and others, 2021a). The SSA is derived by assuming incompressibility and that the vertical normal stress is equal to the overburden pressure, and neglecting vertical shear stress from the Stokes equations and vertically integrating, yielding:

(1)$${\partial M_{ij}\over \partial x_j} - ( \tau_{\rm b}) _i = \rho gH{\partial s\over \partial x_i},\; $$

where i,  j ∈ {1,  2} are the spatial indices in the horizontal x ₁–x ₂ plane, ρ is the ice density, g is the gravitational acceleration, H is the ice thickness, and s is the surface height above sea level. Parameter (τ _b)_i is the basal traction, which is non-zero for grounded ice only, and is described here using a linear friction law:

(2)$$( \tau_{\rm b}) _i = {\hat \beta}^2 v_i ,\; $$

where ${\hat \beta }^2$ is a positive basal friction coefficient, and v _i is the velocity. In (1), the vertically integrated stress tensor M _ij is defined as

(3)$$M_{ij} = 2\bar{\eta}H\left(\dot{\varepsilon}_{ij} + ( \dot{\varepsilon}_{11} + \dot{\varepsilon}_{22}) \delta_{ij}\right),\; $$

where $\dot {\varepsilon }_{ij} = ( {1}/{2}) ( {\partial v_i}/{\partial x_j} + {\partial v_j}/{\partial x_i})$ is the strain-rate tensor defined as the symmetric gradient of the velocity field and δ _ij is the Kronecker delta. The depth-averaged viscosity $\bar {\eta }$ is defined as

(4)$$\bar{\eta} = {1\over 2}\bar{B}\dot{\varepsilon}_{\rm II}^{( 1-n) /n},\; $$

where n is the Glen's flow law exponent (Glen, 1955) and $\dot {\varepsilon }_{\rm II}$ is the second invariant of the strain rate tensor. Herein, we define the depth averaged ice rigidity $\bar {B}$ as

(5)$$\bar{B} = E^{-1/n} \bar{B}_{\rm T},\; $$

where E is an enhancement factor commonly associated with fabric variations that can vary spatially in the horizontal plane. Parameter $\bar {B}_{\rm T}$ is the vertical average of the temperature-dependent ice rigidity B _T(z):

(6)$$\bar{B}_{\rm T} = {1\over H}\int_{b}^{s}B_{\rm T}( z) \, {\rm d}z,\; $$

where b is the ice-shelf draft and B _T(z) is calculated from the depth-varying temperature field T(z) using the standard Arrhenius relation for ice (Cuffey and Paterson, 2010). The boundary condition at the ice front is

(7)$$M_{ij}{\hat n}_j = \left({1\over 2}\rho g H^2 - {1\over 2}\rho_{\rm w} g b^2\right){\hat n}_i,\; $$

where ${\hat {\boldsymbol n}}$ is the unit (outward) normal to the ice front and ρ _w is the seawater density. The first and second terms in the parentheses are the depth integrals of the pressures associated with ice and seawater, respectively. Because the ice pressure exceeds seawater pressure, this boundary condition acts such that it ‘pulls’ the ice shelf seaward. Appropriate Dirichlet conditions for velocity are enforced at all other boundaries. We solve the SSA using the finite-element routine available in Elmer/Ice (Gagliardini and others, 2013), which we modify to incorporate damage as described below.

2.2 Anisotropic damage model

We use an SSA parameterization (Huth and others, 2021b) of the anisotropic creep damage model that was calibrated for glacier ice according to laboratory tests of ice creep to failure under uniaxial tension (Pralong and Funk, 2005). Damage gradually accumulates with time according to a stress-based evolution function, and is incorporated into the vertically integrated momentum balance (1), where it acts to decrease ice viscosity and increase deformation rates for a given stress. The gradual evolution of damage can represent micro/meso-scale crack formation to macro-scale brittle fracture driving the propagation of full-thickness crevasses or rifts, which is consistent with seismic observations (Bassis and others, 2007). We evaluate the damage rate and track the damage variable on the integration points (defined by Gaussian quadrature) within each finite element, because the stresses are most accurate at integration points used during the SSA finite-element solution. We ignore advection for simplicity, which is justified given the short timescale of our simulations.

2.2.1 Creep damage evolution

Damage is represented as a second-order tensor, D, which has three real principal values, 〈D _i〉. Each principal value represents the ratio of the area of cracks to the originally undamaged area along a principal plane normal to the respective principal direction, where the value of 〈D _i〉 ranges from 〈D _i〉 = 0 for undamaged ice to a maximum value of 〈D _i〉 = 1 for fully damaged ice. In the SSA formulation, 〈D ₃〉 = D ₃₃ is always aligned with the vertical x ₃ axis and the other two principal components lie in the horizontal plane.

As described in Pralong and others (2006), a linear transformation between the effective stress $\tilde{\boldsymbol \sigma}$ (i.e. force per unit ice area, ignoring cracks and voids) and the applied stress σ (force per unit area of ice, including cracks and voids) can be defined based on the tensorial damage variable as

(8)$$\tilde{\sigma}_{ij} = {1\over 2}( \sigma_{ik}w_{kj} + w_{ik}\sigma_{kj}) ,\; \quad w_{ij} = ( \delta_{ij}-D_{ij}) ^{-1} .$$

Similarly, an effective strain rate can be defined as

(9)$$\tilde{\dot\varepsilon}_{ij} = {1\over 2}( \dot\varepsilon_{ik}w_{kj}^{-1} + w_{ik}^{-1}\dot\varepsilon_{kj}) ^\prime ,\; $$

where the prime denotes the deviatoric part obtained by subtracting the mean of the diagonal components from each diagonal component of the second-order tensor.

The rate of damage accumulation $\dot {D}_{ij}$ can be obtained based on the objective (Jaumann) rate of damage D as given by (Pralong and Funk, 2005)

(10)$$\dot{D}_{ij} = {\partial{D}_{ij}\over \partial t} = f_{ij} + W_{ik}D_{kj}-D_{ik}W_{kj},\; $$

where t is the time, W is the spin (vorticity) tensor with its Cartesian components W _ij = (1/2)(∂v _i/∂x _j − ∂v _j/∂x _i), and f is the objective damage rate function

(11)$$f_{ij} = B^{\ast}\langle\langle\chi- \sigma_{{\rm th}}\rangle\rangle^r \Bigg(w_{kl} {\hat \xi}^{( 1) }_k {\hat \xi}^{( 1) }_l \Bigg)^k \Bigg({\hat \xi}^{( 1) }_i {\hat \xi}^{( 1) }_j\Bigg).$$

In the above equation, χ is the Hayhurst's equivalent stress

(12)$$\chi = \alpha\langle{\tilde{\sigma}}_1\rangle + \ \beta\sqrt{{3\over 2}\tilde{\sigma}^\prime_{mn}\tilde{\sigma}^\prime_{mn}} + \lambda\tilde{\sigma}_{kk},\; $$

which weights the damage response based on the maximum (most tensile, with the convention that tension is positive) effective principal stress (weighted by α), the Von Mises stress (weighted by β) and the effective hydrostatic stress (weighted by λ = 1 − α − β), where 0 ≤ α,  β,  λ ≤ 1. Damage accumulation is restricted to where χ exceeds the stress threshold, σ _th, according to the Macaulay brackets 〈〈 ⋅ 〉〉 in (11), defined as

(13)$$\langle\langle x \rangle\rangle = \left\{\matrix{ x,\; \hfill & {\rm if}\ x\ge0,\; \hfill \cr 0,\; \hfill & {\rm if}\ x < 0. \hfill }\right.$$

In (11), ${ {\hat {\boldsymbol \xi} }}^{( 1) }$ is the eigenvector corresponding to the maximum effective principal stress in the horizontal x ₁–x ₂ plane, so that damage only accumulates on the plane normal to the ${\hat {\boldsymbol \xi }}^{( 1) }$ direction. The other model parameters, B*,  r and k, are empirical constants. All parameter values are listed in Table 1, and their physical interpretation is described in full in Pralong and Funk (2005) and Duddu and Waisman (2012).

Table 1. Ice and damage parameters common to simulations 1–5

2.2.2 Implementation of damage within the SSA

The SSA only yields vertically integrated deviatoric stresses, whereas damage evolution is defined in terms of the 3-D Cauchy stress field. Therefore, we approximate the Cauchy stress and calculate damage over 21 evenly spaced vertical layers associated with each 2-D integration point (Fig. 2), where we also store the 3-D temperature field. The vertical average of this 3-D damage field is incorporated into the SSA to account for damage-induced weakening of ice. Thus, our quasi-3-D modeling framework accounts for the coupling between the 3-D stress field determined from the 2-D ice flow model and the 3-D damage field describing crevasse and rift propagation.

Figure 2. Flowline depiction of integration points (red), which are each associated with a series of vertical layers (blue) that are distributed evenly along their thickness. Here, we use 21 vertical layers, where 3-D variables such as damage and temperature are represented.

To calculate the Cauchy stress, we first calculate deviatoric stress at the vertical coordinate z of each layer using the flow law (Glen, 1955)

(14)$$\sigma^\prime_{ij}( z) = 2\eta( z) \tilde{\dot{\varepsilon}}_{ij}( z) ,\; $$

where ${ \tilde {\dot {\boldsymbol \varepsilon }}}( z)$ is determined from (9), and the depth-dependent isotropic viscosity is

(15)$$\eta( z) = {1\over 2} E^{-1/n} B_{\rm T}( z) {\dot{\varepsilon}}_{\rm II}^{( 1-n) /n} .$$

Next, we calculate the Cauchy stress as

(16)$$\sigma_{ij}( z) = \sigma^\prime_{ij}( z) -p_{\rm eff}( z) \delta_{ij},\; $$

where p _eff(z) is an ‘effective’ pressure parametrization that accounts for the opposing seawater pressure that penetrates into basal crevasses (Keller and Hutter, 2014)

(17)$$p_{\rm eff}( z) = p_{\rm i}( z) -p_{\rm w}( z) .$$

In the above equation, $p_{\rm i}( z) = \rho g( s-z) -\sigma _{11}^\prime ( z) -\sigma _{22}^\prime ( z)$ is the ice pressure used to derive the SSA under the hydrostatic assumption (Greve and Blatter, 2009) and p _w is the basal water pressure. If a layer is above sea level or is only associated with a surface crevasse (i.e. it is not the basal layer and at least one deeper layer is undamaged) then p _w(z) = 0; else, p _w(z) = ρ _wg(z _sl − z), where z _sl is the sea level elevation. Here, we set z _sl = 0.

Using these Cauchy stresses, the damage tensor components may be updated on integration point layers as detailed in Section 2.2.1. The complete numerical implementation of the 3-D damage update procedure is described in Huth and others (2021b). It includes a Runge–Kutta–Merson scheme for updating D, adaptive time stepping to restrict large changes in damage between time steps, and a non-local integral damage scheme (Duddu and Waisman, 2013) that alleviates mesh dependence by spatially smoothing the change in D each time step over a non-local length scale, l _c. Additionally, we account for rapid damage growth associated with brittle rupture by setting the maximum principal component of 3-D damage to its maximum value of one wherever it meets or exceeds a critical threshold D _crit (Duddu and Waisman, 2013). Subsequent damage evolution on ruptured layers is only allowed through rotation of the damage tensor via the spin terms in (10).

The effect of the 3-D damage field can be incorporated into the vertically integrated SSA stress tensor using the effective strain rate definition as

(18)$$M_{ij} = \int_{b}^{s} 2\eta( z) \left[\tilde{\dot{\varepsilon}}_{ij}( z) + ( \tilde{\dot{\varepsilon}}_{11}( z) + \tilde{\dot{\varepsilon}}_{22}( z) ) \delta_{ij}\right]\, {\rm d}z.$$

However, the above equation rewritten in a simplified form to resemble (3) as

(19)$$M_{ij} = 2\bar{\eta}H\left(\bar{\tilde{\dot{\varepsilon}}}_{ij} + ( \bar{\tilde{\dot{\varepsilon}}}_{11} + \bar{\tilde{\dot{\varepsilon}}}_{22}) \delta_{ij}\right),\; $$

where the effective strain rate ${\bar {\tilde {\dot {\boldsymbol \varepsilon }}}}$ depends on the depth averaged damage ${\bar {\boldsymbol D}}$ as

(20)$$\bar{\tilde{\dot{\varepsilon}}}_{ij} = {1\over 2}( \dot\varepsilon_{ik}\bar{w}_{kj}^{-1} + \bar{w}_{ik}^{-1}\dot\varepsilon_{kj}) ^\prime ,\; \quad \bar{w}_{ij} = ( \delta_{ij}-\bar{D}_{ij}) ^{-1} .$$

By equating (18) and (19), the expression for the depth-averaged damage can be obtained as

(21)$$\bar{D}_{ij} = {\int_{b}^{s}D_{ij}( z) B_{\rm T}( z) \, {\rm d}z\over \int_{b}^{s}B_{\rm T}( z) \, {\rm d}z} .$$

We enforce an additional brittle rupture criterion on the depth-averaged damage ${\bar {\boldsymbol D}}$ to capture rapid damage growth leading to full opening of a rift. This depth averaged rupture criterion uses a unique critical damage threshold $\bar {D}_{{\rm crit}}$ and maximum damage value $\bar {D}_{\rm max}$, typically set close to 1. In practice, $\bar {D}_{\rm max}$ must be set less than the theoretical maximum damage value of one to prevent the SSA from becoming ill-posed. Following Huth and others (2021b), we set $\bar {D}_{{\rm crit}} = 0.8$. At any integration point, if $\langle \bar {D}_1 \rangle \ge \bar {D}_{{\rm crit}}$ or all vertical layers have ruptured, we update all principal components of ${\bar {\boldsymbol D}}$ to reflect that the point has fully failed and now represents a rift. Here, we perform this update by setting $\langle \bar {D}_1 \rangle$ to $\bar {D}_{\rm max}$ and the other principal components of ${\bar {\boldsymbol D}}$ to $\bar {D}_{{\rm max}}-0.05$. This adjustment retains a unique maximum principal component $\langle \bar {D}_1 \rangle$ that allows us to determine rift orientation, which is required to track rift–flank contact (see Section 3). We also set the gravitational driving stress to zero (ρgH (∂s/∂x _i) = 0) for rifted integration points.

2.3 Rift–flank boundary scheme

In this section, we discuss the derivation of the rift–flank boundary condition, its implementation within the FEM–SSA framework, and its modification for mechanically coherent mélange that transmits stress between flanks. The rift–flank boundary condition accounts for the surface forces arising from the contact between rift–flank walls and the presence of mélange and seawater.

2.3.1 Derivation of rift–flank boundary

We derive the boundary condition for rift–flank walls that takes a similar form to the ice-front boundary condition (7), but also accounts for the pressure on flank walls from ice mélange within the rift (Fig. 3a) and full or partial contact between opposite rift–flank walls (Fig. 3b). Partial contact can occur, for example, near the top of rifts due to flexure and rotation of rift flanks (De Rydt and others, 2018; Lipovsky, 2020). We denote the mélange thickness as H _m and the corresponding ice mélange draft as

(22)$$b_{\rm m} = H_{\rm m}{\rho\over \rho_{\rm w}},\; $$

assuming that ice mélange is in floatation and has the same density as glacial ice. Similarly, we define a thickness of contact between flank walls as H _c, which may also have a portion below sea level, b _c. The depth integrated boundary condition for pressure on the rift–flank walls then takes the form

(23)$$\eqalign{M_{ij}{\hat n}_j = \left[{1\over 2}\rho g \left(H^2-H_{\rm c}^2-H_{\rm m}^2\right) - {1\over 2}\rho_{\rm w} g\left( b^2-b_{\rm c}^2-b_{\rm m}^2\right) \right]{\hat n}_i ,\;} $$

where mélange cannot co-exist at the same depth as contact between rift flanks. Like the ice-front boundary condition (7), this boundary force is oriented along the outward normal to the rift–flank wall. Note that this expression is similar to one derived by Larour and others (2014), except that we introduce the H _c and b _c terms that account for pressure from rift–flank contact. Larour and others (2014) also considered friction between rift flanks, as detected for longitudinal rifts along the shear margins where ice shelves meet the bay walls that constrain them. Because this scenario is not applicable to the lateral rifting of interest on the Larsen C Ice Shelf, we do not parameterize friction between flanks here.

Figure 3. Schematic of the mechanics within an open rift that are parameterized by the rift–flank boundary condition. (a) The pressures from seawater (blue) and ice mélange (red) with thickness, H _m, partially oppose the pressure from ice shelf rift flanks (gray). (b) Contact between rift flanks over a thickness, H _c, imparts a similar opposing pressure (not shown) to mélange. We assume H _c is always aligned with the rift–flank surface.

2.3.2 Implementation within the FEM–SSA damage framework

Typically, in the FEM framework, the rift–flank boundary may be embedded into the mesh as a one-dimensional (1-D) interface (i.e. comprised of the edges of 2-D finite elements). The corresponding boundary condition over a 1-D rift–flank boundary element can be applied, similarly to the SSA ice-front boundary condition as discussed in the literature (e.g. Weis, 2001; Greve and Blatter, 2009; Huth and others, 2021a). This involves integrating (23) over each 1-D boundary element Γ_rf so that its contribution to the residual force vector f _iI for the node I is

(24)$$\int_{\Gamma_{\rm rf}}{\phi_I\left[{1\over 2}\rho g\left(H^2-H_{\rm c}^2-H_{\rm m}^2\right) -{1\over 2}\rho_{\rm w} g\left(b^2-b_{\rm c}^2-b_{\rm m}^2\right) \right]{\hat n}_i \, {\rm d}\Gamma},\; $$

where ϕ _I are the standard nodal basis functions. However, evaluating this integral requires explicitly defining the 1-D rift–flank boundary and remeshing as the rift propagates. Instead, we evaluate the contributions to the residual force vector over the 2-D rift zone defined by fully damaged integration points, so that the internal boundary condition can be enforced at run time as the rift propagates, without requiring remeshing. For each 2-D element, we map the contribution of the internal boundary to f _iI as

(25)$${\sum_{r = 1}^{n_{r}}{-{\partial{{\phi_{I}( {\boldsymbol x}_r) }}\over \partial{x_i}} \Bigg[{1\over 2}\rho g( H^2-H_{\rm m}^2-H_{\rm c}^2) -{1\over 2}\rho_{\rm w} g( b^2-b_{\rm m}^2-b_{\rm c}^2) \Bigg]_I}\ w_{r}\vert J_{r}\vert },\; $$

where n _r is the number of fully damaged integration points within the element, x_r is the spatial coordinates, w _r is the weight corresponding to the integration point and |J _r| is the determinant of the Jacobian matrix for the transformation between local (isoparametric) coordinates and global coordinates. To get an intuitive sense of the mapping in Eqn (25), note that it closely resembles (24) if it was converted into a volume integral with the divergence theorem, and evaluated using Gaussian quadrature. However, the difference is that here, the bracketed term containing the depth integrals of the pressures from seawater, mélange, and rift–flank contact is written as a nodal term that describe conditions at rift–flank walls. We discuss how we determine the values of the nodal mélange and contact thicknesses and drafts, used within the bracketed term of (25), in Section 3.2.

A simple example of how (25) enforces the internal rift–flank boundary condition on an ice shelf is given in Figure 4. The red and blue dots within the finite elements (square gridcells) represent fully damaged (rifted) and undamaged integration points, respectively. There are six elements that are fully rifted or fully failed (i.e. all integration points in these elements are fully damaged), while all other elements are undamaged (i.e. they only contain undamaged integration points). The arrows indicate the direction and magnitude of the total contribution from the internal boundary condition to f for each node, by evaluating (25) over all elements. Note that this magnitude decreases as x ₁ increases because the ice thickness decreases as x ₁ increases. mélange and rift–flank contact are both absent in this example, so that there is an open-water boundary condition, and $\bar {D}_{{\rm max}} \approx 1$ so that effectively no stress is transmitted between rift flanks. Recalling that we also remove the gravitational driving stress from rifted integration points, then the only non-negligible contribution of a rift integration point to the model in Figure 4 is through (25). In this case, our scheme behaves similarly to an element-deletion scheme for any fully failed element, wherein the failed element is removed from the mesh and (24) is applied at the new boundaries that appear in its place (i.e. the edges that the failed element had shared with non-failed elements). Note that both our scheme and element-deletion schemes require that the rift width, as represented by integration points, must span at least one element in order to approximate the forces associated with inserting a sharp crack into the mesh. This requirement is satisfied in the non-local damage formulation by using a mesh resolution that is sufficiently smaller than the characteristic non-local damage length, l _c. Furthermore, note that this scheme is most accurate when using consistent element sizes throughout the mesh, e.g. a regular grid. For irregular meshes, the metric term or weighting in (25) may need to be modified to accurately calculate force at the nodes.

Figure 4. Example of the direction (arrows) and magnitude (arrow color and size) of the total contribution from the internal rift–flank boundary condition to the nodal residual force vector, by evaluating the mapping in Eqn (24) over all elements. Each arrow is associated with a node (black points) in the mesh, which is a regular grid of square finite elements (gridcells). Each element edge (black lines) has a length of 1 km. The four dots within each element are integration points. Red dots represent fully damaged (rifted) integration points and blue dots represent undamaged integration points. Here, the domain is a floating ice shelf. There is no mélange or rift–flank contact in this example, so that the rift flanks have an open-water boundary condition like at the ice front. Ice and seawater density match those given in Table 1. Thickness decreases in the x ₁ direction from 410 m at the far left side of the domain to 290 m on the far right side.

2.3.3 Modification for coherent mélange

The above rift–flank boundary scheme can be easily modified to account for mechanically coherent mélange that transmits stress between flanks. For example, in our Larsen C simulations (see Section 3.3), we consider decreasing $\bar {D}_{{\rm max}}$ in some regions as an ad hoc approach to assess the influence of a coherent mélange, in lieu of implementing a more complicated granular rheological model (e.g. Amundson and Burton, 2018). Stress transmission between flanks is also possible without mélange, where horizontal compressive stress could be transmitted between rift flanks that are in contact. While not applicable to our Larsen C simulations, such a situation could occur if a rift is actively closing. This effect could be accounted for using tension/compression asymmetry schemes (e.g. Murakami, 1988).

3. Simulations of rifting on Larsen C Ice Shelf

We perform parametric studies on the rift propagation on the Larsen C Ice Shelf that led to the calving of iceberg A68 in 2017. We start from an initial rift configuration that roughly corresponds to its state in late 2014, which was held through early 2015. At this point, the rift had already propagated from Gipps Ice Rise (GIR), as indicated by the star in Figure 1a marking the rift tip. We run several simulations of the subsequent rift propagation and ice flow evolution from this initial configuration with and without the new rift boundary scheme. The simulations with the rift boundary scheme differ in the application of mélange and flank contact conditions, in order to investigate their role in controlling the rift path. By performing this study on the Larsen C Ice Shelf, we also aim to demonstrate that our damage model can simulate observed rifting. In the following sections, we describe the initial model configuration, the approach used to track rift–flank contact and assign rift–flank boundary conditions during the simulations, and the setup and results for each rifting simulation.

3.1 Initial configuration

To develop the initial model state, we establish the ice geometry, solve for 3-D temperature, and determine fields for the basal friction coefficient, the enhancement factor and an initial damage. While our study focuses on ice-shelf processes, the model domain also comprises all grounded ice within the Larsen C ice-sheet–ice-shelf system. Inclusion of grounded ice is necessary to capture advection into the ice shelf during the temperature solution; it is also necessary because rift propagation during the prognostic (i.e. forward-in-time) simulations can affect ice velocity both throughout the ice shelf and upstream of the grounding line. We determine the initial ice geometry from satellite observations, as described in Appendix A. We use the same 0.5 km node spacing for both this initialization procedure and the rifting simulations.

We determine a 3-D temperature field as it is required to calculate B _T(z) using the standard Arrhenius relation for ice and its vertical average $\bar {B}_{\rm T}$. Recall that B _T(z) influences the 3-D viscosity field in Eqn (15) and $\bar {B}_{\rm T}$ influences the vertically averaged viscosity through Eqns (4) and (5). We summarize our procedure to determine the 3-D temperature field in Appendix B, and the resulting $\bar {B}_{\rm T}$ field is shown in Figure 5a.

Figure 5. Results from the inversion scheme used to separate the three field variables contributing to $\bar {B}$: (a) contribution to $\bar {B}$ from temperature, $\bar {B}_{\rm T}$; (b) $\bar {B}$ from the first inversion; (c) extracted isotropic damage field, $\bar {D}$; (d) $\bar {B}$ from the second inversion; (e) extracted enhancement factor, E, and (f) velocity field from the second inversion.

We determine the basal friction coefficient, enhancement factor and initial damage fields using an inversion procedure that minimizes mismatch between observed and modeled velocities. The observed velocities are derived from a smoothed compilation of 2015 Landsat-8 data (Pope, 2016) with minimal infilling of gaps in coverage using the 2015–16 MEaSUREs data mosaic (Rignot and others, 2017). In lieu of an anisotropic inversion scheme, we define the initial damage field as an isotropic, vertically averaged field $\bar {D}$, to simply incorporate it into the SSA solution as part of the vertically averaged ice rigidity field, $\bar {B}$ in (4) (e.g. Borstad and others, 2013, 2016; Sun and others, 2017). Therefore, we express $\bar {B}$ based on the isotropic vertically averaged damage as

(26)$$\bar{B} = ( 1-\bar{D}) E^{-1/n}\bar{B}_{\rm T} .$$

We designed our inversion procedure to optimize both the basal friction coefficient ${\hat \beta }^2$ (in grounded regions), and the vertically averaged ice rigidity $\bar {B}$ (e.g. Fürst and others, 2015), with additional treatment to separate the contributions of E and $\bar {D}$ to $\bar {B}$. While there is no unique solution for how to separate these variables, we aim to determine a $\bar {D}$ field with sharp gradients aligned with observed fractures, and a smoother E field that describes gradual changes to ice fabric over the domain. A consequence of the smooth E field is that during the prognostic simulations, the rift propagates into a region with smooth spatial variations of ice stiffness, which was inferred without overfitting. This effect helps ensure that inferred ice stiffness influences the simulated rift paths less than changing rift–flank boundary treatments between simulations does, so that we can more clearly investigate how these boundary treatments alone affect rifting. We fully describe the inversion scheme in Appendix C, and the relevant $\bar {B}$, $\bar {D}$ and E fields are plotted in Figure 5.

While the inferred initial damage is 2-D and isotropic, we emphasize that all new damage accumulation during the prognostic simulations is 3-D and fully anisotropic, and is incorporated into the SSA as described in Section 2.2.2. It is possible to convert the inferred 2-D damage into a 3-D field by assuming a vertical distribution of damage, so that the resulting 3-D field can then accumulate additional damage during forward modeling. However, during the prognostic simulations, we do not allow subsequent damage accumulation over areas where non-zero isotropic 2-D damage was inferred, for two reasons: (1) besides the rifting in question, imagery does not show major changes in fracture on the shelf during our time frame of interest; and (2) we are only focused on the propagation of the A68 rift into the undamaged ice near the ice front, where we aim to isolate the effect that varying the rift–flank boundary treatments between simulations has on the rift paths. Isolating this effect requires that damage elsewhere on the shelf remains consistent between simulations, as new damage anywhere on the shelf changes stresses throughout the shelf, impacting rifting. Therefore, we disallow subsequent damage evolution in regions with inferred damage, and also in the immediate vicinity of Gipps and Bawden ice rises, which are pinning points where changes in damage could substantially influence stress throughout the shelf.

Before performing the prognostic simulations, we make two modifications to the initial damage field, which are reflected in Figure 6: first, the inferred damage is unrealistically diffuse so that it does not clearly represent the initial A68 rift, so we redraw it as a sharper rift of fully damaged points, along which we assess the effects of varying rift–flank boundary conditions during the simulations. Second, we initialize an additional region of damage that is observed near the center of the ice front, but not captured during the inversion possibly because it has too minimal of an impact on the smoothed velocity observations. In this region, we assign anisotropic damage corresponding to crevasses oriented normal to the ice flow direction, i.e. opening in the direction of flow. In agreement with observations (yellow arrow in Figs 1a, c), this damage acts to arrest spurious rifting that can otherwise originate from this section of the ice front due to radial spreading during some prognostic simulations. For each point in this region, we assign a vertical damage profile described by fully damaged surface and basal crevasses, separated by an undamaged region consistent with some thickness of ice that is floating in hydrostatic equilibrium. We interpolate this profile to the vertical layers of each integration point, where the undamaged thickness is calculated so that the depth-averaged maximum-principal damage at each point equals 0.5.

Figure 6. Initial damage field used in the prognostic model of the Larsen C Ice Shelf rift propagation. The redrawn initial rift is plotted here with $\bar {D}_{\rm max} = 0.995$. The arrow identifies the additional damage initialized along the front. BIR, Bawden Ice Rise; GIR, Gipps Ice Rise.

3.2 Implementation of the rift–flank boundary scheme

Implementing our rift–flank boundary scheme (Section 2.3) within a prognostic, time-varying simulation requires a method to track the evolution of the rift–flank contact with changes in the rift width. For the simulations here that use the rift–flank boundary scheme, we initially assign full contact between flanks for any new rifting. We assume that the flanks gradually separate as the rift widens because bending effects should cause them to remain in contact near the surface for longer than the base. Other processes may also contribute to enhanced contact between flanks, such as fully or partially calved ice blocks, refreezing of seawater or perhaps a combination of a few degrees of misalignment of the vertical rift plane with the z-axis and buoyancy forces (Walker and Gardner, 2019). However, we assume for simplicity that bending is the primary mechanism of enhanced contact here so that the contact region is always aligned with the top of the rift flanks and b _c = max( − (s − z _sl − H _c),  0). While the bending of rift flanks is not captured within the SSA, we approximate its effect here by linearly decreasing the contact thickness as the rift widens. To do this, we first save the orientation of the rift at full-thickness rupture of an integration point by setting the maximum principal 2-D damage component, $\langle \bar {D}_1 \rangle$, to $\bar {D}_{{\rm max}}$, while the other principal 2-D damage components are set to slightly lower values of $\bar {D}_{{\rm max}}-0.005$. As a proxy for rift widening, we track the accumulated strain, $\varepsilon _{\rm r}$, in the $\langle \bar {D}_1 \rangle$ direction (i.e. the rift-opening direction) on each rifted integration point. A new rift point is initially assigned $\varepsilon _{\rm r} = 0$, and $\varepsilon _{\rm r}$ evolves on subsequent time steps as

(27)$$\varepsilon_{\rm r}^{m + 1} = {\rm max}\left(\varepsilon_{\rm r}^m + \dot{\bar{\varepsilon}}_{\rm r}^m \Delta t^m,\; 0\right),\; $$

where m is the time-step counter and Δt is the size of the time step. Parameter $\dot {\bar {\varepsilon }}_{\rm r}$ is the non-local strain rate in the rift-opening direction at the integration point, calculated as the average of all neighboring integration points within a radius l _c. Without this non-local averaging, $\varepsilon _{\rm r}^{m + 1}$ tends to increase at the center of the rift width compared to the edges, potentially causing error in how the pressure on rift–flank walls is applied. Then, we convert $\varepsilon _{\rm r}$ to a fraction of contact, $F_{\rm c} = {\rm max}( 1-\varepsilon _{\rm r} / \varepsilon _{\rm r}^{\rm max},\; \, 0)$, so that contact linearly varies between 100% at initial full-thickness rupture ($\varepsilon _{\rm r}$ = 0) to 0% when $\varepsilon _{\rm r} \geq \varepsilon _{{\rm r}}^{\rm max}$. Here, we set $\varepsilon _{\rm r}^{\rm max} = 0.04$ for all simulations. We discuss the selection of the value of $\varepsilon _{\rm r}^{\rm max}$ and its influence on rift behavior in Section 4. An example $\varepsilon _{\rm r}$ field is given in Figure 7, which is taken from simulation 3 in Section 3.3. After calculating F _c, we linearly interpolate it to nodes and calculate the nodal contact thickness in (25) as (H _c)_I = (F _c)_I H _I, as well as the corresponding (b _c)_I.

Figure 7. Accumulated strain, $\varepsilon _{\rm r}$, used as a proxy for tracking rift widening, in the rift-opening direction (blue arrows of the inset) as the rift propagates in experiment 3, with $\varepsilon _{\rm r}^{\rm max} = 0.04$.

The nodal mélange thickness is determined similarly to the thickness of rift–flank contact as (H _m)_I = (F _m)_I H _I, where F _m is a constant mélange fraction that we assign at specified integration points and interpolate to nodes. We determine (b _m)_I assuming that the mélange is freely floating. In the simulations, we never allow mélange and rift–flank contact to coexist at the same point, thereby guaranteeing that mélange and contact thicknesses do not overlap within a vertical rift profile.

3.3 Rifting simulations and results

We present five prognostic rifting simulations that demonstrate the model performance under different ‘what-if’ scenarios. The results are reported in Figure 8, where each row (S1–S5) corresponds to one of the simulations (e.g. row S1 corresponds to simulation 1). The columns provide a description of the simulation setup, the rift widening rate ($\dot {\varepsilon }_{\rm r}$) averaged over 0.01 years (~4 d) after the rift begins propagating, and the final damage fields (i.e. rift paths), upon calving.

Figure 8. Results of the five rifting simulations (S1–S5), including (left column) a summary of the initial setup for each experiment; (middle column) the rate of rift widening, $\dot {\varepsilon }_{{\rm r}}$, averaged over the first 0.01 years of rift propagation and (right column) the final maximum principal damage fields, $\langle \bar {D}_1 \rangle$, upon calving. Note $\dot {\varepsilon }_{{\rm r}}$ for S1 is plotted on a different scale than the other simulations. The damage plots use the same colormap as Figure 6. GIR, Gipps Ice Rise.

We describe the motivation, setup and results for each simulation in the subsections below. Most of these simulations test either an extreme end-member of range of possible rift-boundary treatments (e.g. 100% vs 0% mélange fill), or realistic conditions for mélange fill within the rift (e.g. partial mélange fill that is inviscid vs mechanically coherent). Rift treatments are varied between simulations along the initialized portion of the rift, and in some cases for any newly propagated portion of the rift. However, all simulations are assigned an open-water (i.e. no mélange) boundary condition for the portion of the initialized rift that borders Gipps Ice Rise (GIR), which is indicated by the small blue region next to GIR in the Description column of Figure 8 in row S4. In addition to varying the rift treatment between simulations, we also assign each simulation a unique damage stress threshold (σ _th), set low enough to allow rift propagation but high enough to avoid excess damage accumulation elsewhere. Adjusting σ _th in this way yields the sharpest and most realistic rifting possible, thereby optimizing each simulation to potentially match the observed rifting. This adjustment is not strictly necessary, as we can set the same σ _th between all simulations and obtain similar, but more diffuse, rift paths (Supplementary Fig. S1). However, note that if we set σ _th much lower, the resulting excess damage accumulation throughout the domain can change ice-shelf stresses enough to influence the rift path.

Simulation 1: No rift boundary scheme, $\bar {D}_{\rm max} \approx 1$

In simulation 1 (Fig. 8, row S1), we implement the damage model without the rift boundary scheme (except for the open-water boundary next to GIR) and set $\bar {D}_{\rm max} = 0.995 \approx 1$ so that effectively no stress is transmitted between rift flanks; this approach is equivalent to implementing the rift boundary scheme under the end-member assumption that rift flanks are always in contact, or alternatively, always fully filled with inviscid mélange that does not transmit stress. This approach is also consistent with many previous SSA damage models (e.g. Albrecht and Levermann, 2012, 2014; Sun and others, 2017), though the underlying damage models differ. We set σ _th = 0.7 MPa, and the rift propagates to the ice front much more acutely than observed (Fig. 1b). Notably, as rift propagation begins, the simulated maximum ice velocity downstream of the rift is ~80 km a⁻¹ while the respective observed velocities (Pope, 2016) were under ~1 km a⁻¹ (Fig. 5f). The rift widening rate from this simulation is also much greater than the other simulations below.

Simulation 2: No rift boundary scheme, smaller $\bar {D}_{\rm max}$

The setup of simulation 2 (Fig. 8, row S2) is identical to simulation 1 except that we lower $\bar {D}_{\rm max}$ to 0.86 and set σ _th to 0.21 MPa. This simulation tests an ad hoc approach to controlling rifting by adjusting $\bar {D}_{\rm max}$, as performed in a previous study on an idealized geometry (Huth and others, 2021b). Due to the smaller $\bar {D}_{\rm max}$, some stress is transmitted between flanks, which restrains the nascent berg from separating from the ice shelf as quickly as in simulation 1; at the start of rift propagation, the simulated maximum ice velocity downstream of the rift is ~1.2 km a⁻¹, greatly reducing the rate of rift widening as compared to simulation 1. Simulation 2 yields a final rift path that matches observations more closely than simulation 1, illustrating this ad hoc rift scheme as a simple alternative to the internal rift boundary scheme for achieving more realistic rift paths. However, the simulated rift path is not as arcuate (curved) as the observed rift path. Furthermore, this ad hoc scheme represents rifts by means of ‘damage softening’, as opposed to modeling rifts as a discontinuity when using the internal rift–flank boundary scheme. This ad hoc scheme lacks a physical interpretation, so tuning to account for specific rift boundary conditions is challenging.

Simulation 3: Rift boundary scheme, no mélange

In simulation 3 (Fig. 8, row S3), we implement the rift boundary scheme with ‘no mélange’ conditions both within the initialized rift and newly propagated portions of the rift. Thus, this simulation tests the opposite end-member scenario to simulation 1. In this case, $\bar {D}_{\rm max} = 0.995 \approx 1$, and we start the simulation with 100% rift–flank contact at the initialized rift tip that linearly decreases to 0% contact over 30 km from the tip, as indicated by the black-to-white gradient in the Description column of Figure 8 in row S3. We also set σ _th = 0.154 MPa. The simulated maximum ice velocity downstream of the rift at the start of rift propagation matches observations well, at ~0.9 km a⁻¹, resulting in a slowly widening rift. However, unlike simulation 2, essentially no stresses are transmitted between flanks. Instead, these velocities and widening are smaller compared to simulation 1 because the open water boundary condition along much of the rift reduces the net force pulling the flanks apart.

Simulation 4: Rift boundary scheme, weak mélange

Simulation 4 (Fig. 8, row S4) tests the effect of a realistic and inviscid mélange. The setup is identical to simulation 3 except for two modifications: (1) we permanently assign 40% inviscid mélange fill where the initial rift is colored red in the Description column of Figure 8, row S4; and (2) we set σ _th = 0.22 MPa. The mélange effectively does not transmit stress because $\bar {D}_{\rm max} = 0.995 \approx 1$. The inviscid mélange fill reduces the ability of the rift to resist opening, yielding maximum velocities downstream of the rift at the start of propagation of ~1.8 km a⁻¹, which is roughly twice the respective observed velocities. This simulation has an increased rate of rift widening around GIR as compared to simulations 2, 3 and 5 (Fig. 8, column 2), which better simulate the observed rift path. In other words, the nascent berg is rotating away from GIR faster. The resulting rift path lies between the end cases of simulation 1 (i.e. no rift boundary condition, or 100% mélange fill that does not transmit stress) and simulation 3 (i.e. rift boundary condition with no mélange).

Simulation 5: Rift boundary scheme, mechanically coherent mélange

Simulation 5 (Fig. 8, row S5) tests the effect of a realistic and mechanically coherent mélange that transmits stress between rift flanks. The setup of this simulation is identical to simulation 4, except that we lower $\bar {D}_{\rm max} = 0.98$ wherever the 40% mélange fill is applied and set σ _th = 0.165 MPa. The usual $\bar {D}_{\rm max} = 0.995$ is set everywhere else. In the mélange regions, decreasing $\bar {D}_{\rm max}$ from 0.995 to 0.98 locally quadruples the minimum ice stiffness, which scales with $( 1-\bar {D}_{\rm max})$. Thus, the mélange can transmit some stress between flanks and acts to ‘hold’ them together. Simultaneously, the rift pressure boundary condition is active in this simulation, and reduces the net force pulling the flank walls apart from each other, so long as the rift–flank contact is <100%. Thus, simulation 5 is a hybrid of simulations 2–4. It achieves velocities downstream of the rift and a rift path that are consistent with observations.

4. Discussion

Our results demonstrate how rift–flank boundary conditions and mélange strength can influence the rift path. A greater amount of weak (inviscid) mélange fill or contact between rift flanks decreases the rift–flank boundary force (23), which is oriented normal and outward to each rift flank. As demonstrated in simulations 1 and 4, this effect increases rift-widening rates (i.e. increases velocities downstream of the rift), especially near Gipps Ice Rise, which diverts the rift path toward the ice front at a more acute angle than for simulations characterized by smaller rift-widening rates (i.e. smaller downstream velocities). Smaller rift-widening rates result from the opposite conditions – a lower amount of weak mélange fill or contact between flanks – or from stronger mélange fill that can transmit sufficient stresses between flanks to slow them from separating. The rift paths in simulations 3 (‘no mélange’) and 5 (‘strong mélange’), which are associated with smaller rift-widening rates, closely matched the observed rift paths; whereas, the rift paths for simulations associated with weak mélange and increased rift-widening rates (simulations 1 and 4) did not. Though varying amounts of mélange fill were measured within the rift near GIR (Larour and others, 2021), we cannot fully confirm that our ‘strong mélange’ simulation is the most accurate representation of the observed rifting, for two reasons: (1) for simplicity, we approximated the rift system near GIR as a single rift, but satellite imagery suggests it was actually a system of two rifts separated by a thin strip of intact ice until calving, which could have contributed to the overall stress regime of the rift; and (2) the observed mélange fill possibly separated from the rift–flank walls or stretched thin as the rift widened, thereby transitioning to the ‘no mélange case’ over time, but we hold mélange conditions constant over time. To improve the accuracy of our approach for determining the processes that drove the Larsen C rifting, we would need to implement the observed complex rift geometry and spatiotemporally varying mélange conditions.

Our simulations only vary mélange fill and $\bar {D}_{{\rm max}}$, while holding all other damage and rift–flank boundary parameters constant. However, there are likely other combinations of these constant parameters that may yield similar modeled rift paths. For example, there is little observational guidance for choosing an appropriate value for $\varepsilon _{\rm r}^{\rm max}$. Nevertheless, only the most extreme values of $\varepsilon _{\rm r}^{\rm max}$ seem to have a strong impact on rifting. Setting $\varepsilon _{\rm r}^{\rm max}$ too close to its lower limit of zero will effectively eliminate rift–flank contact. Such lack of rift–flank contact will prevent the rift in the ‘no mélange’ simulation (simulation 3) from propagating for any damage stress threshold, σ _th. Conversely, setting $\varepsilon _{\rm r}^{\rm max}$ to a large value effectively prevents rift flanks from separating, resulting in greatly increased velocities downstream of the rift and rapid rift propagation, which can influence the rift path like in simulation 1. In simulations 3–5, we aimed to set $\varepsilon _{\rm r}^{\rm max}$ so that the only effect of rift–flank contact was to consistently enable rift propagation by locally increasing stress at the rift tip, without excessive flank contact that could noticeably influence the rift path. This approach allowed us to solely attribute any differences in rift paths between the simulations to their individual mélange conditions, rather than also having to consider the effects of rift–flank contact on the rift path.

Even though the exact extent of rift–flank contact does vary during and between simulations 3–5, the resulting influence on rift-widening or velocities downstream of the rift is smaller than that from varying the mélange conditions between the simulations. For example, as compared to simulations 3 (‘no mélange’) and 5 (‘strong mélange’), simulation 4 (‘weak mélange’) averages about half the extent of rift–flank contact, but has consistently greater rift-widening rates and velocities downstream of the rift, with roughly twice as high rift-widening rates at the onset of propagation as shown in Figure 8. Therefore, these rift-widening rates and downstream velocities must have increased more by the presence of weak mélange in simulation 4, than decreased by the relatively small extent of rift–flank contact. In fact, it is likely in this case the extent of rift–flank contact is small because weak mélange increases rift-widening rates and causes rift flanks to separate and lose contact more quickly.

While a range of $\varepsilon _{\rm r}^{\rm max}$ may be appropriate, decreasing $\varepsilon _{\rm r}^{\rm max}$ may prevent rifting unless σ _th is decreased as well. Conversely, if increasing $\varepsilon _{\rm r}^{\rm max}$, it may be advantageous to increase σ _th to prevent excess diffuse damage from growing around the rift tip. Unfortunately, both the extent of rift–flank contact, as controlled by $\varepsilon _{\rm r}^{\rm max}$ in our parameterization, and σ _th are poorly constrained. There are few ground-penetrating radar profiles of rift–flank contact available to guide our rift–flank contact parameterization (De Rydt and others, 2018), which are unlikely to be representative of all ice shelves. Moreover, it is not so clear how to even use radar profiles to calibrate $\varepsilon _{\rm r}^{\rm max}$. Another set of potentially poorly constrained parameters are the weights in the Hayhurst stress criteria, a wide range of whose values seem capable of producing similar results. For example, Figure 9 shows how similar results for simulation 5 can be obtained when weighting the Hayhurst criteria entirely by the tensile effective stress (α = 1), rather than using the Von Mises-dominant weighting that we apply otherwise (Table 1). However, it is possible that the Hayhurst weights may have a more substantial effect on diffuse damage accumulation far from the rift tip, which is typically associated with crevassing. We do not assess this effect here because our focus is on understanding the role of rift–flank boundary processes on rift propagation. We only selected one set of model parameters (e.g. $\varepsilon _{\rm r}^{\rm max}$ and the damage constants) sufficient to study these rift–flank processes and simulate observed rifting, but other combinations of model parameters may be equally effective.

Figure 9. Final vertically averaged maximum principal damage field (at 1.57 years) when running simulation 5 (S5) with α = 1,  β = 0 and σ _th = 0.1 MPa.

5. Conclusions

We successfully simulated observed rifting in Larsen C Ice Shelf that led to the calving of iceberg A68, using a combined inverse and forward computational framework based on vertically integrated viscous ice-shelf flow and anisotropic damage formulations. The inversion scheme separates the contributions from damage and the enhancement factor to the inferred ice rigidity. This scheme gives a depth-averaged isotropic damage field that largely resembles observed major rifting and fracture features, and a smoother enhancement factor that may better represent gradual changes in fabric. The forward ice flow and damage model incorporates rift–flank boundary processes as an internal boundary condition, which may be enforced at run time as the rift propagates, within a finite-element framework. These boundary processes include contact between rift flanks and backpressure on rift–flank walls from ice mélange that may also transmit stress between flanks. Through our Larsen C rifting simulations, we found that rift–flank contact, mélange thickness, and mélange strength inside rifts can influence the rift path. Increased contact or weak mélange resulted in a smaller iceberg, and decreased contact or strong mélange resulted in a larger iceberg. Furthermore, the results of our rifting simulations lend support for the argument that gravity-driven viscous stress is sufficient to drive rifting consistent with observations, even without including other mechanical processes, such as ice-shelf flexure in response to the impact of ocean swells.

The current modeling framework is based on the shallow shelf approximation (SSA) of ice flow. While computationally efficient, the SSA ignores bending effects and vertical shear that may influence damage evolution, and a parameterized pressure is required to calculate the 3-D Cauchy stress field needed to evolve damage. It may be advantageous to implement a different flow approximation or an alternative pressure parameterization. Future studies may also consider modifying our inversion procedure to extract anisotropic damage, and convert to 3-D damage according to observed crevasse depths (e.g. from ICESat-2); this could allow further evolution of the inferred damage field within prognostic simulations, which we did not consider herein. Finally, future research should focus on developing more physically based and climate-coupled representations of the rift–flank boundary processes, such as a granular rheology model for mélange and a parameterization of how it grows and decays over time based on environmental forcings. Developing these representations and implementing them within our SSA-damage approach would be a major step toward a comprehensive modeling framework that can simultaneously represent ice flow, melt, rifting and tabular calving. Such a modeling framework would better simulate how ice-shelf weakening may progress, thereby improving projections of ice-sheet evolution and sea-level rise associated with changes in ice-shelf buttressing.