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, Reference Macayeal1989; Huth and others, Reference Huth, Duddu and Smith2021a). 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, Reference Glen1955) 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, Reference Cuffey and Paterson2010). 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, Reference Gagliardini2013), which we modify to incorporate damage as described below.

2.2 Anisotropic damage model

We use an SSA parameterization (Huth and others, Reference Huth, Duddu and Smith2021b) 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, Reference Pralong and Funk2005). 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, Reference Bassis2007). 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 (Reference Pralong, Hutter and Funk2006), 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, Reference Pralong and Funk2005)

(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 (Reference Pralong and Funk2005) and Duddu and Waisman (Reference Duddu and Waisman2012).

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, Reference Glen1955)

(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, Reference Keller and Hutter2014)

(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, Reference Greve and Blatter2009) 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 (Reference Huth, Duddu and Smith2021b). 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, Reference Duddu and Waisman2013) 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, Reference Duddu and Waisman2013). 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 (Reference Huth, Duddu and Smith2021b), 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, Reference De Rydt, Gudmundsson, Nagler, Wuite and King2018; Lipovsky, Reference Lipovsky2020). 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 (Reference Larour2014), except that we introduce the H _c and b _c terms that account for pressure from rift–flank contact. Larour and others (Reference Larour2014) 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, Reference Weis2001; Greve and Blatter, Reference Greve and Blatter2009; Huth and others, Reference Huth, Duddu and Smith2021a). 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.

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, Reference Pope2016) with minimal infilling of gaps in coverage using the 2015–16 MEaSUREs data mosaic (Rignot and others, Reference Rignot, Mouginot and Scheuchl2017). 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, Reference Borstad, Rignot, Mouginot and Schodlok2013, Reference Borstad2016; Sun and others, Reference Sun, Cornford, Moore, Gladstone and Zhao2017). 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, Reference Fürst2015), 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, Reference Walker and Gardner2019). 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, Reference Albrecht and Levermann2012, Reference Albrecht and Levermann2014; Sun and others, Reference Sun, Cornford, Moore, Gladstone and Zhao2017), 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, Reference Pope2016) 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, Reference Huth, Duddu and Smith2021b). 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.