Theoretical stability of ice shelf basal crevasses with a vertical temperature profile

Abstract

Basal crevasses threaten the stability of ice shelves through the potential to form rifts and calve icebergs. Furthermore, it is important to determine the dependence of crevasse stability on temperature due to large vertical temperature variations on ice shelves. In this work, considering the vertical temperature profile through ice viscosity, we compare (1) the theoretical crack depths and (2) the threshold stress causing the transition from basal crevasses to full thickness fractures in several fracture theories. In the Zero Stress approximation, the depth-integrated force at the crevassed and non-crevassed location are unbalanced, violating the volume-integrated Stokes equation. By incorporating a Horizontal Force Balance (HFB) argument, recent work showed analytically that the threshold stress for rift initiation is only half of that predicted by the Zero Stress approximation. We generalize the HFB theory to show that while the temperature profile influences crack depths, the threshold rifting stress is insensitive to temperature. We compare with observations and find that HFB best matches observed rifts. Using HFB instead of Zero Stress for cracks in an ice-sheet model would substantially enlarge the predicted fracture depth, reduce the threshold rifting stress and potentially increase the projected rate of ice shelf mass loss.

1. Introduction

Basal crevasses, which are vertical fractures on the underside of ice, can play an important role in the calving process and thus the stability of ice shelves and marine-terminating glaciers (McGrath and others, 2012a; Colgan and others, 2016; Jeong and others, 2016). A reduction in buttressing through calving or other processes can increase the mass loss from ice sheets and raise global mean sea level by increasing the grounding line flux (Thomas and MacAyeal, 1982; Rott and others, 2002; Rignot and others, 2004; Dupont and Alley, 2005; Pritchard and others, 2012; Haseloff and Sergienko, 2018). Individual basal crevasses can induce surface crevasses, create surface depressions when sufficiently deep that may enable surface meltwater ponding, reduce the ability of ice shelves to provide back stress to upstream grounded ice and penetrate through the full ice thickness to form rifts (Pralong and Funk, 2005; McGrath and others, 2012b; Luckman and others, 2012; Child and others, 2021). When spaced periodically, basal crevasses can potentially stabilize ice shelves from breakup through dampening stresses transmitted through high-frequency elastic-flexural waves (Freed-Brown and others, 2012). The evolution of basal crevasses has been modeled by the balance of ice shelf and hydrostatic ocean stresses (Zarrinderakht and others, 2022), as well as idealized ocean dynamics and the mass balance of melting and freezing (Jordan and others, 2014). Basal crevasses may play a central role in the flexure-driven calving of marine-terminating glaciers seen in Greenland and Canada (Wagner and others, 2014; Murray and others, 2015; Wagner and others, 2016; Benn and others, 2017). In Antarctica, basal crevasses may initiate rifts (Jeong and others, 2016; Joughin and others, 2021) that can propagate across the ice shelf and calve icebergs (Lipovsky, 2020). These icebergs can transport fresh meltwater equatorwards and threaten biodiversity of islands in the Southern Ocean (Huth and others, 2022).

While the calving of icebergs is likely caused by multiple mechanisms, we focus on the transition from basal crevasses to full-thickness fractures, also referred to as rifts. In the absence of ample surface meltwater such as surface melt ponds, basal crevasses are theoretically more vulnerable than surface crevasses to full-thickness penetration and cause rift initiation (Lai and others, 2020). The depth-averaged deviatoric stress formulations of Nye (Nye, 1955), Weertman (Weertman, 1973) and the zero toughness or half-space formulations of Mode I Linear Elastic Fracture Mechanics (LEFM) (van der Veen, 1998a) all predict basal crevasses to be about nine times deeper than dry (water-free) surface crevasses. The magnitude of lithostatic stress only decreases as basal crevasses propagate upwards, yet increases as surface crevasses propagate downwards, making the initiation of rifts more likely due to basal crevasses than dry surface crevasses. Thus, we study basal crevasses as the precursors of rifts in the absence of strong atmospheric forcing (Morris and Vaughan, 2003; van Wessem and others, 2023).

Driven by the importance to predict rift initiation on ice shelves, here we extend analytical and numerical models to both predict the ice shelf threshold stress $R_{xx}^\ast$ for rifts to initiate from Mode I (tensile) basal crevasses, and include a vertical temperature profile, which we show can substantially affect fracture predictions, yet is currently neglected in fracture parameterization in ice-sheet models (Goelzer and others, 2020; Seroussi and others, 2020). We consider three theories for Mode I basal crevasses to propagate into rifts: (1) the Zero Stress approximation (Jezek, 1984; Benn and others, 2007; Nick and others, 2010), (2) LEFM for basal crevasses (van der Veen, 1998a; Tada and others, 2000; Lai and others, 2020) and (3) the Horizontal Force Balance (HFB) (Buck, 2023). On ice shelves, the Zero Stress approximation forms rifts from the intersection of basal crevasses with a surface crevasse, and requires the largest amount of tensile stress among the three theories for a rift to form. LEFM considers an isolated basal crevasse such that stress can concentrate at the crack tip and predicts the smallest amount of tensile stress required for rift initiation. The major limitation of current LEFM-based basal crevasse theory used in this paper (van der Veen, 1998a) is its neglection of the restoring force at the deformed ice-ocean interface, leading to non-physical predictions of the rifting stress. The HFB theory considers the union of surface and basal crevasses for rift formation in a material with zero strength, similar to the Zero Stress approximation. However, unlike the Zero Stress approximation, HFB theory does not violate static force balance. We find that the rifting stress predicted by HFB for a range of temperature profiles is the same as the tensile stress at the ice shelf front with no sea ice buttressing; a rifting stress slightly greater than that of LEFM but much less than that of the Zero Stress approximation.

Our fracture models’ deviations from standard implementations and validation of the Zero Stress approximation (Jezek, 1984; Benn and others, 2007; Nick and others, 2010) and Mode I LEFM for basal crevasses (van der Veen, 1998a) are summarized below. First, we include the depth variation of the resistive stress R _xx due to vertical temperature variation in all theories and compare the results. We chose a linear temperature profile for simplicity, and then show the effects of a different temperature profile. Second, we enforce the Zero Stress approximation to uphold HFB and include vertical temperature variations, following a procedure similar to Buck (2023) and obtaining a simple analytical result. Third, we shift the focus from crevasse depth prediction to rift formation prediction, analyzing results in terms of stress required for basal crevasses to unstably propagate and initiate rifts or calving events. Fourth, we validate rift formation predictions with an existing rift catalog (Walker and others, 2013) on the Ross Ice Shelf (RIS) and Larsen C Ice Shelf (LCIS). We verified that the deviation in resistive stress between the 1D fracture theory and the 2D Shallow-Shelf Approximation (SSA) (MacAyeal, 1989) is less than $10\%$, thus ensuring validity of the 1D flow assumption in the regions of interest.

2. Fracture models of basal crevasses with vertically varying ice temperature

The basal crevasse evolution problem is illustrated schematically in Figure 1a, where crack growth may either terminate at a stable length, or propagate unstably to form a rift, with outcome subject to the chosen fracture theory. Figure 1b shows the differences between existing isothermal basal crevasse depth prediction theories: the Zero Stress approximation, LEFM and a new HFB model in (Buck, 2023). While these theories differ for a range of crack depths, the largest discrepancy is near the sea level, where basal crevasses can unstably propagate and form rifts. Although the basal crevasse to rift transition is challenging to precisely measure given the current spatiotemporally sparse observations (see Appendix B), Figure 1c through f motivates our study of rift initiation associated with basal crevasse vertical propagation through observed increases in resistive stress upon rift formation. Importantly, ice shelves are not isothermal, and basal crevasse depths are sensitive to the ice shelf vertical temperature profile (Rist and others, 2002; Borstad and others, 2012; Lai and others, 2020). We analyze several fracture models (van der Veen, 1998a; Nick and others, 2010; Buck, 2023) and the stress required to form rifts, considering a vertically linear temperature profile for simplicity. We assume that the base of the ice shelf is held at T _b = −2 °C, and take a linear temperature profile up to the surface temperature T _s as predicted by RACMO (van Wessem and others, 2018). To conclude this section, we show that having a Robin temperature profile (Robin, 1955) does not change the results of this study.

Figure 1. (a) Schematic of crevasses propagating stably or unstably and forming a full-thickness fracture called a rift. Crevasse depths d _s,  d _b, thickness H, coordinate system and piezometric head at sea level z _h are illustrated. (b) Several previously existing isothermal crevasse depth predictions versus depth-averaged resistive stress $\bar {\! R}_{xx}$ normalized by the analytical depth-averaged ice tongue resistive stress $\bar {\! R}_{xx}^{IT}$ for H = 300 m. Circular red dots are LEFM basal crevasse depth numerical predictions, solid lines are basal crevasse depth analytical theory and dash-dotted lines are surface crevasse depth analytical theory. Rifts initiate either where d _b = H for LEFM or at the stars that denote the intersection of surface and basal crevasse tips. Horizontal Force Balance (HFB) and Zero Stress require the existence of both surface and basal crevasses on ice shelves for theoretical consistency, whereas LEFM treats an isolated basal crevasse. Subfigures (c) to (f) show a potential instance of the basal crevasse-to-rift transition (Jeong and others, 2016; Joughin and others, 2021) over Pine Island Ice Shelf during January to May 2019. (c) An estimate of the ratio $\bar {\! R}_{xx}/\bar {\! R}_{xx}^{IT}$ over Pine Island Ice Shelf found using ice velocity data averaged over the month of January 2019 (Wuite and others, 2021). (d) SAR backscatter image at 50 m resolution from 5th January 2019, and a close-up showing the terminus region of Pine Island Ice Shelf where a fracture (can be surface crack, surface expression of basal crack or rift) is dimly visible. (e) The equivalent of (c) for May 2019. (f) A backscatter image from May 2019 where a rift is clearly visible. Grounded ice (Haran and others, 2013), grounding lines (Mouginot and others, 2017b) and calving fronts are denoted by masking, a dashed black line and a solid black line, respectively.

The historical focus on crevasses was centered on surface crevasses, with Nye's Zero Stress (Nye, 1955) and LEFM (Smith, 1976). However, the use of radar confirmed the presence of crevasses extending upwards from the bottom of the ice (Jezek and others, 1979; Jezek and Bentley, 1983). The surface crevasses theories were later extended to study basal crevasses (Jezek, 1984; Benn and others, 2007; Nick and others, 2010; van der Veen, 1998a). The Zero Stress approximation we refer to in this paper considers these extensions (Jezek, 1984; Benn and others, 2007; Nick and others, 2010) of Nye's Zero Stress approximation that consider both a surface crevasse and a hydrostatic seawater pressure inside a basal crevasse.

Since the setup of these basal crevasse depth theories share many basic assumptions, we will first present the commonalities, and then explore each individually. As shown in Figure 1a, given a 2D coordinate system with x as the horizontal dimension and z as the vertical dimension, assuming incompressibility and a stress-free upper surface, we can write the net longitudinal stress σ _n of van der Veen (1998a) at the approximately vertical basal crevasse interface as

(1)$$\sigma_{\rm n}\left(z \right) = R_{xx} \left(z \right)- p_{{\rm l}} \left(z \right) + p_{{\rm w}} \left(z \right),\; $$

where R _xx(z), p _l(z) and p _w(z) are the along-flow component of ice shelf resistive stress defined in Cuffey and Paterson (2010), ice lithostatic pressure and hydrostatic water pressure, respectively. As in the Shallow Shelf Approximation of Stokes flow (MacAyeal, 1989), we assume that there is negligible vertical shear stress in the ice due to negligible shear stress on the surface and basal boundaries. By setting z = 0 at the ice shelf base and positive upwards as shown in Figure 1a, we define the pressure terms as p _l = ρ _ig (H − z) and water pressure p _w = ρ _wg max(z _h − z,  0), with gravitational acceleration g = 9.8 m s⁻², vertically integrated ice density ρ _i = 917 kg m⁻³ and ice thickness H. The piezometric head $ z_{\rm h} = H \rho_{\rm i} / \rho_{\rm w} $ as defined in Nick and others (2010) depends on the water density; for this study, we assume constant saltwater density ρ _w = 1028 kg m⁻³.

The way that we account for vertical temperature profile T(z) is through the ice hardness B(T) in the effective viscosity. Modeling ice as a Non-Newtonian power-law fluid (Glen, 1958), the effective viscosity is,

(2)$$\mu = {B\left(T \right)\over 2} \dot{\epsilon}_e^{({1}/{n})\, - 1} \approx {B\left(T \right)\over 2} \dot{\epsilon}_{xx}^{({1}/{n})\, - 1},\; $$

where the effective strain rate $\dot {\epsilon }_e = \sqrt {\dot {\epsilon }_{ij} \dot {\epsilon }_{ji} /2}$, i.e. the second invariant of the strain rate tensor $\dot {\epsilon }_{ij}$, is dominated by the along-flow component $\dot {\epsilon }_{xx}$ in comparison to the across-flow and shear terms. Denoting τ _xx as the along-flow component of the deviatoric stress, the along-flow resistive stress $R_{xx} = 2\tau _{xx} = 4\mu \dot {\epsilon }_{xx}$ can be written as

(3)$$R_{xx} = 2 B \left(T \right)\dot{\epsilon}_{xx}^{{1}/{n}}.$$

Based on lab experiments, the Glen's flow law gives n = 3 (Glen, 1958), and B(T) (LeB Hooke, 1981) can be expressed as

(4)$$B \left(T \right) = B_0 \exp \left[{T_0\over T} - {C\over \left(T_{\rm r} - T \right)^k} \right],\; $$

where B ₀ = 2.207 Pa ⋅ yr^1/n, T ₀ = 3155 K, T _r = 273.39 K, k = 1.17 and C = 0.16612 K^k are constants determined from empirical fit (LeB Hooke, 1981). We discuss sensitivity to rheology choices in the Robin Temperature Profile subsection, the Discussion and Conclusions section, and Appendix F.

2.1. Zero stress approximation

Nye's theory (Nye, 1955) assumes that a crevasse will stop propagating when further infinitesimal crack growth would put the net longitudinal stress σ _n at the crack tip into compression. Originally conceived for surface crevasses only, it has been extended to model basal crevasses by approximating the stress exerted from the ocean on crevasse walls as hydrostatic (Jezek, 1984; Benn and others, 2007; Nick and others, 2010; Duddu and others, 2020). This so-called Zero Stress approximation generates crevasses where there is net tensile stress; on ice shelves, surface and basal crevasses are both allowed to propagate. While previous Zero Stress applications assume isothermal ice (Nye, 1955; Jezek, 1984; Benn and others, 2007; Nick and others, 2010; Duddu and others, 2020), we have extended the theory to consider depth-varying temperature. As derived in Appendix C, a basal crevasse on a freely floating ice shelf with vertical temperature variation can unstably propagate to the surface when the resistive stress is above the threshold

(5)$${\bar{R}_{xx} \over \bar{R}_{xx}^{IT} } \geq 2 {\rho_{\rm w}\over \rho_{\rm i}} {d_{\rm b}^\ast \over H} {\overline{\tilde{\! B} \left(T \right)} \over \tilde{\hskip-2pt B} \left (T \left(d_{\rm b}^\ast /H \right)\right)},\; $$

where the overline $\bar {q}$ represents a depth-averaged value for the variable q. The dimensionless numbers involved are the ratio of densities ρ _w/ρ _i; the unstable basal crevasse depth $d_{\rm b}^\ast$ relative to ice thickness H; the non-dimensional ice hardness $\tilde {\! B} ( T ) \equiv {B( T) }/{B( T( z = 0) ) }$; and the non-dimensional ice hardness at the unstable basal crevasse depth $\tilde {B} ( T ( d_{\rm b}^\ast /H ) )$; the depth-averaged resistive stress $\bar {\! R}_{xx}$ relative to the resistive stress of a 1D unbuttressed ice tongue $\bar {\! R}_{xx}^{IT} \equiv \left (1-\rho_i / \rho_w \right ) \rho_{i} g H/2 $ as derived by Weertman (1957) from the force balance for an ice tongue . In the Zero Stress approximation, all basal crack depths less than the unstable basal crevasse depth $d_{\rm b}^\ast$ are stable, while those greater than or equal to $d_{\rm b}^\ast$ will result in rift formation.

As shown in Appendix C, we write the result of depth-dependent resistive stress in terms of depth-averaged resistive stresses such that we can use the same ratio $\bar {R}_{xx}^\ast / \bar {R}_{xx}^{IT}$ between all theories in this paper. Note that when we treat the ice as isothermal or utilize a depth-averaged resistive stress instead of the depth-varying resistive stress, the ice hardness ratio ${{\overline{\tilde{\! B} \left(T \right)}}}/\tilde {B} ( T ( d_{\rm b}^\ast /H ) )$ is 1, and the maximum basal crevasse depth is from the base to sea level $\left ( \rho_w d_b^* \right ) / \left ( \rho_i H \right ) = 1 $, making the right-hand side of Eqn (5) equal to 2. Thus, under the Zero Stress approximation, a basal crevasse will propagate to the sea level, intersecting a surface crevasse to form a rift in isothermal ice when the depth-averaged resistive stress is twice that of an isothermal ice tongue, as shown in Figure 1b.

The simplicity of the Zero Stress approximation comes with limitations. The rifting stress threshold we derived in Appendix C to get Eqn (5) does not include stress concentration near crack tips, the material strength of ice, accumulation and melt, nor crevasse-induced stresses in ice. The lack of stress concentration and zero material strength is applicable in the limit of closely spaced crevasses (De Robin, 1974; Weertman, 1974), where the spacing between crevasses is much less than the individual crevasse depths (Weertman, 1977). Observations indicate basal crevasse spacing to be roughly one to several ice thicknesses (Luckman and others, 2012; McGrath and others, 2012a; Lawrence and others, 2023), breaking the densely spaced crevasses assumption. Recent extensions of the Zero Stress approximation include non-zero material strength (Benn and others, 2007) and accumulation and melting effects (Bassis and Ma, 2015; Huth and others, 2021). In the limit of densely spaced crevasses with negligible flexural stress, the effects of crevasse-induced stress on crack depth has been included to satisfy a HFB argument (Buck, 2023). We present an approach to generalize this approximate crevasse-induced stress for vertically varying ice temperature in the Horizontal force balance section.

2.2. Linear elastic fracture mechanics

In contrast to the Zero Stress approximation, the LEFM framework, first applied to surface crevasses by Smith (1976) and later applied to basal crevasses by van der Veen (1998a), considers an isolated basal crevasse with stress concentration near the crack tip and assumes small-scale yielding (Zehnder, 2012). It has been shown that LEFM agrees with the analytical approach of including stress concentration near crack tips by Weertman (1973) for small crack depths (Buck and Lai, 2021). Unlike Weertman's infinite thickness assumption (Weertman, 1973), LEFM comes with the advantage of accounting for a prescribed finite thickness (van der Veen, 1998a). The LEFM rifting threshold $\bar {\! R}^\ast _{xx}$ for an isothermal ice shelf with traction-free upper and lower surfaces has been reported by Zarrinderakht and others (2022) and Lai and others (2020). Building on previous work (van der Veen, 1998a; Tada and others, 2000; Lai and others, 2020), we extended the LEFM analysis across a range of surface temperatures applicable to Antarctica and present the rifting threshold $\bar {R}^\ast _{xx}$ for a vertically varying ice temperature.

The criterion for Mode I (tensile) crevasse propagation in LEFM is that a crack can propagate so long as the stress intensity factor K _I is at least as large as the fracture toughness K _Ic, K _I ≥ K _Ic. Following Tada and others (2000); van der Veen (1998a); Lai and others (2020), the stress intensity factor for a Mode I basal crack can be written as

(6)$$K_{{\rm I}} = \int_{0}^{d_{\rm b}} {2 \sigma_n( z) \over \sqrt{\pi d_{\rm b}}} {G( \zeta,\; \, \tilde{\!\! d}_{\rm b}) \over \left(1 - \, \tilde{\!\! d}_{\rm b} \right)^{3/2} \left(1 - \zeta^2 \right)^{1/2}} {\rm d}z,\; $$

where d _b is basal crack depth, H is ice thickness, $\zeta \equiv z/{d_{\rm b}}$ is a dimensionless height that is one at the crack tip, $\ {\tilde{\!\! d}_{\rm b}} \equiv {d_{\rm b}}/{H}$ is the dimensionless basal crack depth relative to the ice thickness, σ _n is the net longitudinal stress defined in Eqn (1), and $G( \zeta ,\; \, \, \tilde {\!\! d}_{\rm b})$ is a dimensionless weight function given by Tada and others (2000) (page 71) which is 99% accurate for any crack depth d _b. For the fracture toughness of glacial ice, we use the experimental value of K _Ic = 150 kPa$\sqrt {\rm{m}}$, which is found to be independent of temperature (Litwin and others, 2012). For rifting, we follow the non-dimensionalization and stability criteria in the Supplemental material of Lai and others (2020): rifts are formed from unstable basal cracks when the tensile stress is large enough such that $\tilde {K}_{\rm I} \equiv {K_{{\rm I}}}/{\rho _{\rm i} g H^{3/2}}$ increases monotonically with $\tilde {d}_{\rm b}$ (see Fig. 13). One caveat is the assumption that the ice includes pre-existing initial flaws with size $d_{\rm b}^{\rm i}$, which depends on the stress state, thickness and fracture toughness. The theoretical estimates of $d_{\rm b}^{\rm i}$ are 3 orders of magnitude smaller than the ice thickness for H = 300 m, as discussed in Appendix D, and thus one would expect the existence of these pre-existing flaws.

We can approximately explain the dependence of rifting stress on ice temperature by extending a torque equilibrium argument by Zarrinderakht and others (2022) with isothermal ice to account for the temperature structure effects. As in Appendix C of Zarrinderakht and others (2022), the rifting stress threshold can be determined by zero torque (or moment) into the page associated with a deep basal crevasse. Mathematically,

(7)$$\eqalign{\tau_y \,=\, \hat{y} \cdot \left(\underline{r} \wedge \underline{F} \right) & = 2 \int_b^s \left(R_{xx} \left(z \right)- \rho_{\rm i} g \left(s - z \right)\right)\left(s - z \right)\rm{d}{\it z} \cr & - 2 \int_b^0 \rho_{\rm w} g z \left(s - z \right){\rm d} z = 0,\;} $$

where s and b are the freely floating ice surface and base, the origin z = 0 is at sea level, and the ocean restoring force due to the vertically deformed ice-ocean interface is neglected for equivalent comparison. The only difference from Zarrinderakht and others (2022) is the depth-variation of the resistive stress due to temperature R _xx(z). For simplicity, we approximate the ice hardness function in Eqn (4) with B _a in Eqn (E8) (see Fig. 2a),

(8)$$B_{\rm a}\left(T\left(\tilde{z}\right)\right)\approx B_0 \exp \left[{T_0\over T_b} \right]\exp \left[{\tilde{z} + {{\rho_{\rm i}}\over{\rho_{\rm w}}} \over \tilde{z}_0} \right],\; $$

with $\tilde {z} = {z}/{H} = 0$ at sea level, constants B ₀ = 2.207 Pa ⋅ yr^1/n, n = 3, T ₀ = 3155 K, surface temperature T _s, basal temperature T _b and the dimensionless e-folding length scale assuming a linear temperature profile $\tilde {z}_0 \equiv ( ({T_0}/{T_{\rm b}}) ( 1- {T_{\rm s}}/{T_{\rm b}}) ) ^{-1}$. Computation of the moments with vertical temperature structure as in (7) were also presented in Buck (2024). Substituting (8) into the resistive stress (3) and calculating the torque equilibrium (7), the temperature-structure-dependent LEFM rifting stress threshold can be analytically derived,

(9)$${\bar{\! R}_{xx}^\ast \over \bar{\! R}_{xx}^{IT}} = {{{2}\over{3}} \left(2 - {{\rho_{\rm i}}\over{\rho_{\rm w}}} \right)\over 2 \tilde{z}_0 \left(1 - \left(\tilde{z}_0 \left(\exp{\left({1}\over{\tilde{z}_0} \right)} - 1 \right)\right)^{-1} \right)}.$$

The numerator of the equation is the same as that in Zarrinderakht and others (2022), while the denominator is the contribution due to a linear temperature profile. While this equation is based on an approximate ice hardness function, it approximately explains the numerical LEFM solution (see Fig. 2b), capturing the role of linear vertical temperature dependence through one new dimensionless variable, the dimensionless e-folding length scale $\tilde {z}_0$. Figure 2b shows that the dimensionless LEFM rifting stress threshold decreases with warmer surface temperature, consistent with Figure 4c.

Figure 2. (a) Comparison of depth-averaged ice hardness value $\bar {B}$ normalized by the constant B ₀ from LeB Hooke (1981) versus surface temperature, assuming a linear temperature profile. The red curve numerically integrates (4), while the blue curve is an analytical integral of the approximated ice hardness profile B _a (8). (b) Comparison of the rifting stress threshold versus surface temperature, assuming a linear temperature profile. The LEFM numerical solution calculated using (6) is shown with red dots, while the analytical solution (9) based on torque equilibrium with an approximated ice hardness B _a (8) is in solid blue. We plot the limit of $T_{\rm s} \xrightarrow {} -2$ °C for the rifting stress from torque equilibrium as it converges to the isothermal case. The difference between the analytical and numerical results is non-negligible, thus for the rest of the paper we use the numerical LEFM results. That said, the analytical result (9) offers interpretable insight by capturing the role of the linear temperature profile with one additional dimensionless variable, the e-folding length scale $\tilde {z}_0$.

One of the largest limitations of this version of LEFM is neglecting the oceanic restoring force associated with the vertical displacement at the ice-ocean boundary (Jimènez and Duddu, 2018; Huth and others, 2021; Zarrinderakht and others, 2022), which is expected to be particularly important for deep basal crevasses. Accounting for this oceanic restoring force in a LEFM framework is the subject of future work. However, the effect of buoyancy is included in the numerical simulations of an isolated basal crevasse in Buck and Lai (2021), which predicts the same rifting stress threshold as the HFB theory presented in the next section.

2.3. Horizontal force balance

In the standard Zero Stress approximation, the depth-integrated horizontal force per unit width (F _x) at the crevassed (x = x _c) and non-crevassed location (x = x _c + Δx) are unbalanced (see Figs 3b,d). If this force is unbalanced ($\sum F_x\neq 0$), according to the Newton's Second Law this net force acting on the ice in the control volume enclosed by the black dashed lines in Figure 3a should lead to acceleration of the ice masses, which is inconsistent with our Stokes flow assumption. The main assumption in the Zero Stress approximation that led to the inconsistency is the neglection of crevasse-induced stresses in the unbroken ligament between the surface and basal crevasse tips (Figs 3a,b). Recent work by Buck (2023) shows a new model that includes the crevasse-induced stresses in the unbroken ice (Fig. 3c) by incorporating a HFB argument. In this section, we apply the same argument for an ice shelf with vertically varying temperature. We use an Eulerian control volume approach (the box enclosed by the black dashed lines in Fig. 3a). Note that by satisfying HFB we mean that the total force acting on the control volume argument is zero ($\sum F_x = 0$), satisfying Newton's Second Law. As assumed in the Zero Stress approximation, we consider the limiting case of a densely crevassed ice shelf where flexural stresses are negligible.

Figure 3. (a) Schematic showing force balance between a downstream, unperturbed location x = x _c + Δx and the crevassed location at x = x _c with crevasse depths, surface d _s and basal d _b. The unbroken ligament is depicted of length L = H − d _s − d _b. The unperturbed background longitudinal stress at x = x _c + Δx in this example is $\sigma _{xx}( x_{\rm c} + \Delta x ) = R_{xx} - p_{\rm l} = 0.9 \bar {R}_{xx}^{IT} - \rho _{\rm i} g ( H - z )$. For the remaining subfigures, the top and bottom row correspond to the Zero Stress approximation and HFB (Buck, 2023), respectively. Figures (b) and (c) show the longitudinal stress profiles normalized by ρ _ig H at the crevassed location in the Zero Stress approximation and HFB (Buck, 2023). Figures (d) and (e) correspond to the stress difference profiles Δσ _xx ≡ σ _xx(x _c + Δx) − σ _xx(x _c) normalized by $\bar {\! R}_{xx}^{IT}$. The Zero Stress approximation does not uphold force balance because the stress difference is positive or zero for all depths, yet HFB is defined to uphold the horizontal force balance constraint as written in (a). We include crack depths for the Zero Stress approximation, where the blue curve intersects zero and labeled with the superscript 0, as well as the deeper HFB crack depths labeled with the superscript HFB.

Our fixed control volume in Figure 3 has vertical boundaries at the ice shelf surface and base, and horizontal boundaries at the symmetry plane of a basal crevasse, x = x _c, and at a nearby downstream location x = x _c + Δx with the same ice thickness. The temperature profile that dictates the vertical profile of R _xx(T(z)) is depicted in Figure 1. We start with Stokes flow, yet inertial terms may be non-negligible (Bassis and Kachuck, 2023). Taking glaciostatic balance as in Lindstrom and MacAyeal (1987) as the vertical force balance,

(10)$$\partial_z \sigma_{zz} = \rho_{\rm i} g,\; $$

we integrate from a level z to the surface H, i.e. σ _zz = ρ _ig (z − H), and solve for pressure

(11)$$p\left(z\right)\equiv -\sigma_{zz} + \tau_{zz} = \rho_i g \left(H -z\right) + \tau_{zz}$$

using a stress-free upper surface boundary condition. From the downstream side of the control volume at x = x _c + Δx, the longitudinal ice shelf Cauchy stress can be written as

(12)$$\eqalign{\sigma_{xx} \left(x_{\rm c} + \Delta x,\; z \right)& \equiv -p + \tau_{xx} = -\rho_{\rm i} g \left(H - z \right)- \tau_{zz} + \tau_{xx} \cr & = - \rho_{\rm i} g \left(H - z \right) + R_{xx}( z) ,\;}$$

where the deviatoric stress tensor τ _ij = σ _ij + pδ _ij by definition has zero trace tr (τ _ij) = 0 (e.g. see Section 12.1 of Rudnicki (2014)), and thus τ _xx = −τ _zz in a 2D (x,  z) system without assuming incompressibility. Note that this stress distribution (12) is valid for an uncrevassed location, where there is no basal crevasse-induced flexural stress.

At x = x _c, in order to satisfy HFB ($\sum F_x = 0$), the extra compressive stresses in the unbroken ligament of length L between the surface and basal crevasse tips, induced by the crevasses themselves, need to be considered. We parameterize this crevasse-induced stress to satisfy three conditions:

1. 1. the zero material strength assumption of the Zero Stress approximation (Nye, 1955; Jezek, 1984; Benn and others, 2007; Nick and others, 2010),

2. 2. continuity of stress at crack tips (Buck and Lai, 2021) and

3. 3. HFB (Buck, 2023).

Previous zero stress models (Nye, 1955; Benn and others, 2007; Nick and others, 2010; Duddu and others, 2020) use condition 1 to solve for crack depth; we maintain that the stress field obey condition 3 for static equilibrium, and take conditions 1 and 2 in the limit of a densely crevasses ice shelf with negligible flexural stresses.

The longitudinal stress that is consistent with conditions 1 and 2 can be written in a piecewise fashion corresponding to a dry surface crevasse of depth d _s (H − d _s ≤ z ≤ H), a water-filled basal crevasse of depth d _b (0 ≤ z ≤ d _b) and an ice ligament between the two (d _b ≤ z ≤ H − d _s),

(13)

Here, cR _xx(z) represents the sum of the background resistive stress and the crevasse-induced compressive stress in the unbroken ice ligament, d _b ≤ z ≤ H − d _s.

We note that while a precise numerical profile of the stresses in the ice ligament can differ from the parameterized form cR _xx, the deviation would be negligible in the limit of a very small unbroken ice ligament and so the rifting stress threshold will be accurate despite the parameterized stress profile (13). Investigations of material strength effects, crevasse-induced flexural stresses and more complicated physical descriptions of the stress within the unbroken ice ligament are subjects for future work.

HFB on an Eulerian or fixed control volume, defined with negligible inertial term as the Stokes equation, can be written as

(14)$$0 = \int_{0}^{H} \int_{x_{\rm c}}^{x_{\rm c} + \Delta x} \left[\partial_x \sigma_{xx} + \partial_z \tau_{zx} \right]{\rm d}x {\rm d}z.$$

The negligible shear stress on the upper and lower surfaces simplifies Eqn (14) such that the sum of the horizontal forces per unit width into the page on our control volume is zero,

(15)$$\int_0^H [ \sigma_{xx}\left(x_{\rm c} + \Delta x \right)- \sigma_{xx} \left(x_{\rm c} \right)] {\rm d}z = 0.$$

Solving for the constant, c, in Eqn (13) and finding a relation between surface and basal crack depths will close the system of equations, allowing us to determine crevasse depths given a normalized resistive stress $\bar {\! R}_{xx} / \, \bar {\! R}_{xx}^{IT}$.

2.3.1. Crevasse depths

Following Buck and Lai (2021), we use continuity of stress at crack tips, z = H − d _s and z = d _b, as a constraint to determine the constant c and then a relation between the crevasse depths. At the tip of the surface crevasse, we have that

(16)$$0 = - \rho_{\rm i} g d_{\rm s} + c ( R_{xx}) \vert _{z = H-d_{\rm s}},\; $$

which easily resolves the constant as

(17)$$c = {\rho_{\rm i} g d_{\rm s}\over ( R_{xx}) \vert _{z = H-d_{\rm s}}}.$$

At the basal crevasse tip, we have that

(18)$${- \rho_{\rm i} g H + \rho_{\rm w} g d_{\rm b} = - \rho_{\rm i} g \left(H - d_{\rm b} \right) + {\rho_{\rm i} g d_{\rm s}\over ( R_{xx}) \vert _{z = H-d_{\rm s}}} ( R_{xx}) \vert _{z = d_{\rm b}},\;} $$

which can readily be simplified to a dimensionless relation between basal and surface crack depth,

(19)$${d_{\rm b}\over d_{\rm s}} = {\rho_{\rm i}\over \rho_{\rm w} - \rho_{\rm i}} {B \left(T\vert _{z = d_{\rm b}} \right)\over B \left(T\vert _{z = H-d_{\rm s}} \right)} \quad \hbox{for the general case}.$$

Thus, the temperature profile T(z) affects the relative surface to basal crevasse depth through B(T (z)), with colder surface temperatures creating larger crevasse depth ratio d _s/d _b.

Finally, having solved for the constant, we may now impose the force balance constraint of Eqn (15) with the stress expressions in Eqns (12) and (13). Defining the dimensionless variables $\tilde {\! d}_{\rm b} = d_{\rm b}/H$, $\tilde {\!\! d}_{\rm s} = d_{\rm s}/H$ and $\tilde {z} = z/H$, the dimensionless force balance can be written as

(20)$$\eqalign{ \frac{\bar{\! R}_{xx} }{\bar{\! R}_{xx}^{IT} } &= \frac{\rho_w}{\rho_w - \rho_i} \,\ {\tilde{\hskip-3pt d}_s^{\,2}} + \frac{\rho_w}{\rho_i}\ \tilde{\hskip-3pt d}_b^{\, 2} + \frac{ \tilde{\hskip-4pt d}_s }{ \frac{1}{2} \left ( 1 - \frac{\rho_i}{\rho_w} \right ) \tilde{\! B} \left ( T |_{\tilde{z}=1-\hskip 2.5pt \tilde{\hskip-2pt d}_s} \right ) } \cr &\quad\times \int^{1-\hskip 2.5pt \tilde{\hskip -2.5pt d}_s}_{\hskip 2pt \tilde{\hskip-2.5pt d}_b} \tilde{\! B}(T(\tilde{z})) d \, \tilde{\! z} \ \text{for the general case}.}$$

Equations (19) and (20) form a system of two equations with the surface and basal crack depths as the two unknown variables, given that $\bar {R}_{xx}/\bar {R}_{xx}^{IT}$ is known.

For isothermal ice shelves, the ice hardness functions of Eqn (19) become constants as in Eqn (E6), and we have

(21)$${d_{\rm b}\over d_{\rm s}} = {\rho_{\rm i}\over \rho_{\rm w} - \rho_{\rm i}} \quad \hbox{for the isothermal case}.$$

The equation has an analytical solution (Buck, 2023),

(22)$$\eqalign{\tilde{\hskip-3pt d}_b & = \tilde{\hskip-3pt d}_s \frac{\rho_i}{\rho_w - \rho_i}, \quad \text{and}, \quad \tilde{\hskip-3pt d}_b = \frac{\rho_i}{\rho_w} \left ( 1 - \sqrt{1 - \frac{\bar{\! R}_{xx}}{\bar{\! R}_{xx}^{IT} }} \right ) \cr & \text{for the isothermal case},}$$

which predicts that rifts would form when the background resistive stress reaches that of a freely floating ice shelf without buttressing, i.e. $\bar {\! R}^\ast _{xx}/\bar {\! R}_{xx}^{IT} = 1$. Importantly, we did not set $\ \bar {\! R}_{xx}^{\ast } = \bar {\! R}_{xx}^{IT}$; this arises naturally as the force balance solution to rift formation where $\ \tilde{\hskip-3pt d}_b + \ \tilde{\hskip-3pt d}_s = 1$,

(23)$$\tilde{\hskip-3pt d}_b + \ \, \tilde{\hskip-3pt d}_s = 1 - \sqrt{1 - \frac{\bar{\! R}_{xx}}{\bar{\! R}_{xx}^{IT} }} = 1 \quad \text{for the isothermal rift case}.$$

When we include the vertical temperature profile, we compute the ice hardness function given vertical temperature variation numerically. We iterate through temperature profiles and basal crevasse depths to solve for surface crevasse depth through the equation residual of Eqn (19). Having numerically obtained a relation between $\tilde{d_{\rm b}}$ and $\tilde{d_{\rm s}}$, we use these values to solve for $\bar {\! R}_{xx} / \bar {\! R}_{xx}^{IT}$ in Eqn (20). We plot dimensionless basal crack depth as a function of $\bar {\! R}_{xx}/\bar {\! R}_{xx}^{IT}$ for the linear vertical temperature case in Figure 4. When the crack is stable, the crevasse depth is instantaneously determined by the stress state. This is consistent with the Zero Stress approximation (Benn and others, 2007; Nick and others, 2010) and van der Veen (1998a)'s basal crevasse LEFM consideration of neglecting time-dependent fracture propagation, but is a possible future extension (Lawn, 1993).

Figure 4. Predicted crevasse depths normalized by ice thickness H for (a) the Zero Stress approximation, (b) Horizontal Force Balance and (c) Linear Elastic Fracture Mechanics given a linear temperature profile in the vertical direction from a basal temperature of −2 °C to surface temperature T _s. The y-axes are non-dimensional height from the ice shelf base, and the x-axes are depth-averaged resistive stress normalized by depth-averaged ice tongue resistive stress. The solid lines are the normalized height of basal crack tips measured from the base upwards, and the dash-dotted lines are the normalized height of surface crevasse tips measured from the surface downwards. Rifting is represented by yellow stars where basal and surface crevasses meet at sea level and by open stars for unstable basal crack depth $d_{\rm b}^\ast$ (e.g. Eqn (5)), respectively. For LEFM, rifting occurs when the basal crack, without the presence of a surface crack, fully penetrates the ice thickness. The isothermal cases presented here at T = −2 °C of the Zero Stress approximation and Horizontal Force Balance are analytical, while all other results are numerical.

2.3.2. Rifting stress

We find that the crevasse to rift transition for cold and warm ice shelf surface temperatures occurs at roughly the same critical stress. Thus, HFB can be approximated by a simple analytical rifting criterion that is independent of the vertical temperature profile,

(24)$${\bar{\! R}^\ast _{xx}\over \bar{\! R}_{xx}^{IT}} = 1.$$

This result from force balance can be explained intuitively. The transition from uncracked to cracked ice causes a force F _x acting on the ice through changes in gravitational potential energy U,

(25)$$F_x = {\partial U\over \partial x} = \int_{-{\rho_{\rm i}\over \rho_{\rm w}}H}^{-{\rho_{\rm i}\over \rho_{\rm w}} H + d_{\rm b}} \left(\rho_{\rm w} - \rho_{\rm i} \right)g z {\rm d}z - \int_{\left(1 - {\rho_{\rm i}\over \rho_{\rm w}} \right)H - d_{\rm s}}^{\left(1 - {\rho_{\rm i}\over \rho_{\rm w}} \right)H} \rho_{\rm i} g z {\rm d}z.$$

The terms on the right-hand side represent the changes in force due to water replacing ice and ice replaced by a free surface for basal and surface crevasse opening, respectively. We reset z = 0 to sea level to highlight the similarity to the compressive buckling problem in equation (B5) of Coffey and others (2022). In the case of rifting with cracks meeting at sea level, this force is

(26)$$F_x = - {1\over 2} \left(1 - {\rho_{\rm i}\over \rho_{\rm w}} \right)\rho_{\rm i} g H^2 = - H \bar{ R}_{xx}^{\ast }.$$

To uphold static HFB ΣF _x = 0, we see that the tensile force in the ice must be the negative of this value. Thus, the resistive stress required for rifting is that of a 1D ice tongue, as in Eqn (24).

2.4. Comparison between the three fracture models

The comparison between the three fracture models is presented in Figure 4. In order of smallest to largest rifting threshold $\bar {\! R}_{xx}^\ast /\bar {\! R}_{xx}^{IT}$, or highest to lowest vulnerability to rifting, we have LEFM, HFB and the Zero Stress approximation for all temperatures analyzed in this study. We note that the influence of temperature is distinct between the three theories; while HFB has negligible temperature dependence, colder surface temperatures lower the rifting threshold $\bar {\!R}_{xx}^\ast /\bar {\! R}_{xx}^{IT}$ for the Zero Stress approximation yet increase the rifting threshold for LEFM. However, for HFB in Figure 4b, colder surface temperatures cause a decrease (increase) in basal (surface) crevasse depth, yet leave the rifting threshold stress unchanged. In summary, the rifting stress ratios that account for a vertically varying temperature profile for each theory and the formula we use to compute the dimensionless stress from ice shelf data are:

$$\eqalign{& \text{Zero Stress}: \quad \frac{\bar{\! R}_{xx}^* }{\bar{\! R}_{xx}^{IT} } = 2 \frac{\rho_w}{\rho_i} \frac{d_b^*}{H} \frac{\overline{\tilde{\! B} \left ( T \right )} }{\tilde{\! B} \left ( T \left ( d_b^*/H \right ) \right )} \cr& \text{LEFM}: \text{Analytical approximation (1%-5%):} \cr& \frac{\bar{R}_{xx}^*}{\bar{R}_{xx}^{IT}} = \frac{\frac{2}{3} \left ( 2 - \frac{\rho_i}{\rho_w} \right )}{2 \!\tilde { \hskip2pt{ z}_0} \left ( 1 - \left (\!\!\! \tilde{ {\hskip3pt z}_0} \left (\exp{\left (\! \frac{1} {\tilde{\hskip2pt z_0}} \right )} - 1 \right ) \right )^{-1} \right )} \cr& \text{HFB}: \frac{\bar{\! R}_{xx}^* }{\bar{\! R}_{xx}^{IT} } = 1 \cr& \text{Data}: \frac{\bar{R}_{xx} }{\bar{\! R}_{xx}^{IT} } = \frac{2 \bar{\! B} \dot{\epsilon}_{xx}^{{1/n}}}{\frac{1}{2} \left ( 1 - \frac{\rho_i}{\rho_w} \right ) \rho_i g H}\,.}$$

Note that when we compute the dimensionless stress with ice shelf data, we first confirm in Figure 8 that regions are approximately 1D to uphold the plane strain assumption of LEFM, and estimate the stress state prior to rifting through estimating the unbroken ice thickness as discussed in detail in Appendix A.

We can understand the temperature effects on crack depth and rifting stress through intuitive explanations of each theory, with extended discussions in Appendices C through E. One way to understand the temperature effects is decomposing the depth-varying resistive stress into a depth-averaged component and a vertically varying component, $R_{xx} ( z ) = \bar {\! R}_{xx} + R_{xx} ' ( z )$. Figure 12 shows that R _xx^′(z) is negative toward the ice base and positive toward the ice surface. The Zero Stress approximation is determined by the net longitudinal stress σ _n in Eqn (1) at the crack tip; thus, due to the sign of R _xx^′(z), we see smaller basal crack depths toward the base and smaller depth-averaged rifting stress thresholds $\bar {\! R}_{xx}^\ast$ for colder T _s in Figure 4a. In contrast, numerical LEFM results in Figure 4c show smaller crack depths and larger depth-averaged rifting stress thresholds $\bar {\! R}_{xx}^\ast$ for colder ice. The analytical explanation is provided in the LEFM section and Appendix D. Finally, our HFB theory inherits the assumption of local force balance at crack tips from the Zero Stress theory, yielding smaller basal (larger surface) stable crack depths for colder T _s with larger vertically varying resistive stress R _xx(z) as shown in Figure 4b. The rifting stress threshold mostly results from stable basal and surface crevasse depths meeting at sea level. The effects of temperature on the ratio of the surface to basal crack depth (19) vanish as the surface and basal cracks approach sea-level, as seen in Figure 4b. Thus, the rifting stress threshold in HFB is well-approximated by the isothermal result of $\ \bar {\! R}_{xx}^\ast /\bar {\! R}_{xx}^{IT} = 1$. Overall, the effect of temperature on basal crevasse depth and rift initiation depends strongly on the chosen theory.

2.5. Robin temperature profile

To check the sensitivity of our conclusions to temperature profiles that are not linear, we run through the analyses of this paper assuming a Robin temperature profile (Robin, 1955). While this solution is strictly valid for an ice divide, the curvature of the profile is closer to observed temperature from borehole data in some cases than a linear temperature profile (Thomas and MacAyeal, 1982; Rist and others, 2002; Craven and others, 2009; Sergienko and others, 2013; Tyler and others, 2013). The goal of this exercise is not to create highly realistic temperature profiles by modeling the computationally expensive temperature evolution and advection from ice divides to ice shelves, but is instead meant as a sensitivity test of the results to the assumed temperature profile. As with the example plotted in Figure 5, the Robin family of temperature profiles take the form

(27)$$T \left(\tilde{z} \right) = T_{\rm s} + {q\over k} \sqrt{{\pi \kappa H_{\rm d}\over 2 \dot{a}}} \left[\hbox{erf} \left(\sqrt{{\dot{a} H_{\rm d}\over 2 \kappa}} \right)- \hbox{erf} \left(\tilde{z} \sqrt{{\dot{a} H_{\rm d}\over 2 \kappa}} \right)\right],\; $$

with surface temperature T _s, thermal diffusivity κ ≡ k/(ρ _ic _p) ≈ 10⁻⁶ m² s⁻¹ defined by thermal conductivity to ice density and specific heat of ice, rescaled vertical coordinate $\tilde {z} = {z}/{H}$ with value 0 at the ice base and 1 at the surface, basal heat flux q, ice divide thickness H _d ≈ 1000 m and snowfall rate $\dot {a}\approx 0.1$$\ {\rm m\, a^{-1}}$ based on Fowler and Ng (2020). We note that the ice divide thickness H _d is set to match the Sandhäger and others (2005) profile at sea level, as demonstrated in Figure 5. We match the temperature profiles near sea level, as this region is important in determining if basal crevasses propagate to form rifts. Considering the ice-ocean temperature at the bottom of ice shelf T _b = −2 °C, we have

(28)$$T_{\rm b} = T_{\rm s} + {q\over k} \sqrt{{\pi \kappa H_{\rm d}\over 2 \dot{a}}} \hbox{erf} \left(\sqrt{{\dot{a} H_{\rm d}\over 2 \kappa}} \right),\; $$

Substituting (28) into (27) gives a simple form of the Robin profile,

(29)$$T \left(\tilde{z} \right) = T_{\rm s} + \left(T_{\rm b} - T_{\rm s} \right)\left(1 - {\hbox{erf} \left(\tilde{z} \sqrt{{\dot{a} H_{\rm d}}\over{2 \kappa}} \right)\over \hbox{erf} \left(\sqrt{{\dot{a} H_{\rm d}}\over{2 \kappa}} \right)} \right).$$

The remainder of the Robin profile analyses for each fracture theory is carried out the same way as presented in the body of this paper. In the next section, we show that the rifting stress thresholds are not significantly impacted by using a Robin temperature profile versus linear temperature profile.

Figure 5. Theoretical ice shelf temperature profiles.

3. Comparison with observations

As the three fracture theories predict distinct critical stresses that drive the basal crevasse to rift transition, we evaluate the applicability of each theory by comparison with observed rift locations. We analyze our results in two ways. First, we plot the predicted rift locations on MODIS MOA (Scambos and others, 2007; Haran and others, 2018) compared with the rifts previously mapped by Walker and others (2013), labeled as ‘true rifts’ and colored in blue on Ross Ice Shelf (RIS) in Figure 6. The goal of these rift formation theories is to maximize the overlap between the predicted and true rifts, colored in green in Figure 6. Because we do not have values of strain rate or surface temperature at the time of rifting, the estimated stress state (Fig. 6e) uses modern surface temperature (van Wessem and others, 2018) and strain rate (Wearing, 2017) values, with limitations discussed in Appendix A. In Figure 6, we see that on the RIS rifts identified by Walker and others (2013), Zero Stress with a vertical temperature profile (Fig. 6b) underpredicts known rifts as shaded in blue. Similarly, LEFM with a depth-averaged resistive stress $\bar {\! R}_{xx}$ (Fig. 6c), with analytical result given by Zarrinderakht and others (2022), overpredicts rifts into areas they were not observed as shaded in red. LEFM with a vertical temperature profile (Fig. 6a) and HFB (Fig. 6d) are the most accurate theories for these RIS rifts. Their differences are small enough to warrant more observations to distinguish which theory is most applicable.

Figure 6. Figures (a) to (d) are map views in kilometers of observed (marked in blue) and theoretically predicted rifts (marked in red) on the Ross Ice Shelf, overlain on MODIS MOA 2014 (Haran and others, 2018; Scambos and others, 2007), with correctly predicted rifts in green. Rift formation theories are (a) LEFM with depth-varying stress due to temperature variation R _xx(z), (b) Zero Stress with R _xx(z), (c) LEFM with depth-averaged stress $\bar {R}_{xx}$ and (d) Horizontal Force Balance with R _xx(z). Figure (e) is the estimated resistive stress ratio used for panel (a–d) with axes in meters. (f) Observed rifts were identified by Walker and others (2013) and shown in blue.

Second, we construct a crack stability plot of the critical stress $\bar {\! R}^\ast _{xx}/\bar {\! R}_{xx}^{IT}$ for rift formation for each theory as a function of the ice surface temperature T _s, as shown by the curves in Figure 7. Noticeably, surface temperature has a negligible effect on the rifting threshold for the HFB; the rifting stress is that of a freely floating ice tongue, $\bar {\! R}_{xx}^\ast = \bar {\! R}_{xx}^{IT} \equiv ( 1 - {\rho _{\rm i}}/{\rho _{\rm w}} ) \rho _{\rm i} g H/2$ (light blue line in Fig. 7). Comparing the rift formation stress criteria, the depth-averaged Zero Stress approximation (dashed blue line) requires a resistive stress $200\%$ of that of a freely floating ice tongue to cause rifts, while the isothermal LEFM theory as presented by Zarrinderakht and others (2022) (dashed magenta line) requires only $74\%$ of that of a freely floating ice tongue to initiate rifts. It is our goal to constrain this substantial uncertainty in the rift initiation stress threshold of these two theories, presented by the dashed lines of Figure 7.

Figure 7. Non-dimensional resistive stress required to initiate rifts versus surface temperature, assuming a linear temperature profile. The solid lines account for vertically varying temperature structure through depth-varying stress R _xx(z), whereas the dashed lines use depth-averaged resistive stress $\bar {\! R}_{xx}$. For a given theory, the regions above and below the curve are predicted as a rift and non-rift, respectively. Ice shelf stresses averaged over regions that are not rifted and obey the 1D fracture assumptions are plotted in orange symbols with mean and standard deviation. These data are shown in map view in Figure 11. We use the average thickness in the non-rift area as H to calculate $\bar {\! R}_{xx}^{IT}$. In this figure, surface temperatures colder than −25 °C are on the Ross Ice Shelf (RIS), whereas those warmer are on the Larsen C Ice Shelf (LCIS). While the Zarrinderakht and others (2022) line in dashed magenta looks compelling with the intact ice shelf data in orange, map view analyses (see Figs 6c, 10c) show that it overpredicts rifts.

We quantitatively compare the rift criteria with the observed rifts on the RIS in Figure 6 and Larsen C Ice Shelf (LCIS) in Figure 10. We identify the extensional, 1D flow regions excluding the high-strain-rate rift locations on the RIS and LCIS in Figure 11, and plot the mean depth-averaged resistive stress and surface temperature across these regions in orange symbols in Figure 7 with one standard deviation of uncertainty due to the variation of resistive stress and surface temperature across these regions. The bulk of the non-rift ice shelf data should lie below the curves in Figure 7. While LEFM with depth-averaged resistive stress may look accurate in Figure 7, the overprediction in red is clear in Figures 6c and 10c. Additionally, the lack of force balance from isothermal Zero Stress (dashed blue line) and Zero Stress with vertical temperature variation (solid green line) invalidates these predictions on theoretical grounds. This lack of force balance manifests as an underprediction of rifting, as discussed in Appendix A and demonstrated in Figure 16. Therefore, within the uncertainty of our data as discussed in the Uncertainties of the data-model comparison subsection and Appendix F, our results on two ice shelves suggest LEFM and HFB considering the vertically varying ice temperatures are the most accurate theories for predicting the rift locations on the RIS and LCIS.

Figures 14, 15, 16 and 17 show that the conclusions drawn from the linear profile are also applicable to the Robin temperature profile. The rift formation stress curves with linear and Robin temperature profiles in Figure 16 are similar. The relative accuracy of each theory in Figures 14 and 15 is comparable to Figures 6 and 10, respectively. Overall, the conclusions of the paper with a linear temperature profile are unchanged given our estimate of a Robin temperature profile.

4. Discussion and conclusions

4.1. Threshold stress for the transition from basal crevasses to rifts

In this paper, we have determined several rift formation stress criteria as functions of surface temperature for linear temperature profiles. We then use remote-sensing and model output data to determine which theory best predicts observed rifts. We find that Zero Stress with either depth-averaged or vertically varying resistive stress underpredicts rifts, due in part to the inconsistency that the formulations do not uphold force balance in a control volume. On the other hand, we find that LEFM with depth-averaged resistive stress overpredicts rifts. Our result shows that on the RIS and LCIS, HFB and LEFM with vertically varying resistive stress are the most accurate theories for correctly predicting rifts and non-rifts. Further distinction between these two theories is inhibited by the number of rifts and uncertainty of current data products used in this study. However, given that buoyancy is expected to stabilize deep basal crevasses (Logan and others, 2013; Zarrinderakht and others, 2022) and is not available analytically for the LEFM models, we expect the rift initiation stress threshold of LEFM models to increase. Thus, for the initialization of rifts in ice-sheet models, we recommend using the simple, analytical, physically consistent theory of HFB, with rift initiation threshold $\bar {\! R}^\ast _{xx}/\bar {\! R}_{xx}^{IT} = 1$ (Eqn (24)). This rifting threshold has been numerically predicted for an isolated basal crevasse (Buck, 2023). Further, this study validates this result against observed rifts on RIS and LCIS. In the Robin temperature profile subsection, we compared the results between linear and Robin temperature profiles (Robin, 1955), both showing a negligible temperature dependence of this rifting threshold.

4.2. Stability of ice tongues

A question may then naturally arise: How can the freely floating ice tongue stress be sufficient to form rifts, yet we see ice tongues exist in nature? It is important here to draw a distinction between idealized 1D, constant thickness, zero yield-strength ice tongues and those found in nature: any perturbations to the stress field, mass balance, spatially varying thickness and crack geometry can influence the stability of real ice tongues. For example, sea ice or ice mélange can stabilize these structures through potentially providing a force buttressing the ice shelf and dampening ocean waves that would otherwise induce ice shelf flexural stress (Vaughan, 1995; Bromirski and others, 2010; Sergienko, 2010, 2013; Hulbe and others, 2016; Massom and others, 2018; Miles and others, 2020; Gomez-Fell and others, 2022), with recent work arguing this idea may be most applicable for thin ice (Bassis and others, 2024). Recent observational work analyzing the Eastern Antarctic Peninsula found that $94\%$ of ice shelf calving occurred during or shortly after the removal of sea ice (Christie and others, 2022), highlighting the importance of external atmospheric and oceanic conditions. The stress threshold and environmental sensitivity of ice tongue calving is crucial for predicting calving and potential subsequent changes in dense water formation, carbon export and biological productivity, as evidenced by the 2010 Mertz Glacier Tongue calving (Kusahara and others, 2011; Shadwick and others, 2013; Ohshima and others, 2016).

4.3. Model limitations

This idealized study is subject to several limitations. We assume 2D (x,  z), incompressible, homogeneous density, elastic modulus and fracture toughness (Rist and others, 2002) ice shelves with zero across-flow strain rate and zero shear strain rate. These theories do not include local thickness variation, creep closure, sub- and super-buoyant flexure (Benn and others, 2017), ice front bending stresses (Reeh, 1968), grain size dependence of yield stress (Ranganathan and others, 2021), crack tip shielding (Clayton and others, 2022), snow accumulation, basal melting (Bassis and Ma, 2015; Kachuck and others, 2022; Buck, 2023) and marine ice accumulation. While we take a constant density, we note recent work that suggests the importance of vertical density variation for surface crevasses (Gao and others, 2023).

To be consistent with the plane strain assumptions of our fracture theories when comparing with observations, we develop a strain rate criterion in Appendix A to validate locations where the flow is approximately 1D and Mode I fracture is applicable. In reality, ice shelves can also have their stress states altered due to 3D effects such as shear fractures (van der Veen, 1999) and torque from ocean currents (Gomez-Fell and others, 2022; Huth and others, 2022). The interactions of the ocean, sea ice, mass balance, the location of additional cracks or other additional stresses may determine the stability of the ice tongues and shelves in nature by modulating the ice shelf stress. Effects such as the hydrographic conditions in basal crevasses, crevasse-, tidal- or tsunami-induced flexural stresses (Walker and others, 2013; Brunt and others, 2011; Bromirski and others, 2017; Gerstoft and others, 2017), viscous creep closure, fatigue failure (Zehnder, 2012), potentially non-negligible process zones (Zehnder, 2012), realistic ice rheology or the time-dependent evolution (Lawn, 1993) of basal crevasses evolving into rifts are neglected. We leave the coupling of these processes to future work.

One inherent source of model uncertainty is the formation mechanism of the rifts on RIS and LCIS that we analyze in this study. As mentioned in the introduction, we can not definitively claim that the rifts were created by basal crevasse propagation. Surface crevasse hydrofracture is a relatively well-studied problem, and ample surface meltwater can form a rift on an ice shelf (Weertman, 1973, 1974; van der Veen, 1998b; Banwell and others, 2013; Duddu and others, 2020; Lai and others, 2020). However, seawater-filled a basal crevasses are predicted to be roughly nine times deeper than dry surface crevasses across several isothermal crevasse theories (Nye, 1955; Nick and others, 2010; Weertman, 1973; van der Veen, 1998b,1998a). This implies that in the absence of strong surface melt (Morris and Vaughan, 2003; van Wessem and others, 2023) basal crevasse is the likely rift formation mechanism.

4.4. Uncertainties of the data-model comparison

Both data products and ice rheology contribute sources of uncertainty to this study. We discuss uncertainty in data products in Appendix F, where we estimate an upper bound in uncertainty and see that the largest source of uncertainty is in the estimation of ice thickness on the Ross Ice Shelf (Morlighem and others, 2020). Another source of uncertainty comes from the ice rheology. We note that there is uncertainty in the flow law exponent n (Goldsby and Kohlstedt, 2001; Qi and others, 2017; Bons and others, 2018; Millstein and others, 2022) and therefore in the parameters of Eqns (2) and (4) (Zeitz and others, 2020). This will affect both our estimation of stress as shown in Figure 6e, and the rift initiation stress threshold curves in Figures 7, 16 for Zero Stress and LEFM. However, one great benefit of the HFB model is an insensitivity to rheology in the limit of rifting. This is exemplified in Figure 4b: while crack depths are different for different surface temperatures, the rift initiation threshold is approximately the same across various vertical temperature profiles. Furthermore, as demonstrated in the Robin temperature profile subsection the rift initiation stress threshold between linear and Robin temperature profiles is the same (see Figs 7, 16). Additionally, as HFB best matches observations (see Figs 6, 10), the leading rift initiation theory of this study demonstrates a unique resilience to rheological uncertainty.

4.5. Future applications

Areas of future research that may build on this study are at least threefold: (1) initializing rifts with the HFB model in numerical ice sheet/shelf models, (2) analyzing the impacts of fracture model on processes governing the ice shelf stability; HFB predicts much deeper crack depths than the Zero Stress approximation and (3) incorporating higher temporal-resolution observations to observationally constrain the rifting stress threshold.

First, this study provides a simple analytical rift formation stress threshold and crack depths that can be coupled with numerical ice-sheet simulations. Damage mechanics has shown great promise (Huth and others, 2021), capturing the path of rift propagation to form iceberg A68 from Larsen C (Huth and others, 2023). However, this model requires the specification of a starter rift (in this case near the Gipps Ice Rise) to predict its horizontal propagation. Coupling rift propagation models with our HFB theory for rift initialization, $\,\bar {\! R}_{xx}^{\ast }/\bar {\! R}_{xx}^{IT} = 1$ can potentially advance simulations’ ability to model rifts from initiation to iceberg detachment.

Second, another interesting extension would be embedding HFB into the analysis of processes that require analytical fracture theory, such as Bassis and Walker (2012); Pollard and others (2015); Bassis and Ma (2015). For example, because the HFB model predicts much deeper cracks compared with the Zero Stress approximation that was used in Pollard and others (2015), implementing our HFB fracture model instead could yield substantially different outcomes, making ice shelves more vulnerable to fracture than predicted by Pollard and others (2015).

Third, future work monitoring the rift initiation process with high temporal resolution can provide stronger constraints on the rifting stress threshold. As shown in Figure 1c and discussed in Appendix B, the best estimates of stress in this study rely on monthly averaged strain rate data. However, the speed of horizontal rift propagation can be much faster than what can be observed from satellite imagery (Olinger and others, 2024). While it is promising that we see a clear increase in surface strain rate and thus stress estimates in Figure 1c, it is clear that higher temporal-resolution measurements of strain rate, cross-rift altimetry and/or field-based data of the basal crevasse to rift transition would help to constrain these rifting stress threshold theories.

This study advances the common fracture theories used in glaciology through the incorporation of vertical stress variations due to the temperature dependence of the ice hardness, and the HFB model that properly accounts for HFB. The different rifting stress theories coupled with ice shelf simulations could lead to distinct calving predictions. Importantly, the Zero Stress approximation would largely under-predict the crevasse depths and stability compared with the HFB model, and underestimate the mass loss of ice shelves over the coming centuries.