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, Reference McGrath2012a; Colgan and others, Reference Colgan2016; Jeong and others, Reference Jeong, Howat and Bassis2016). 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, Reference Thomas and MacAyeal1982; Rott and others, Reference Rott, Rack, Skvarca and Angelis2002; Rignot and others, Reference Rignot2004; Dupont and Alley, Reference Dupont and Alley2005; Pritchard and others, Reference Pritchard2012; Haseloff and Sergienko, Reference Haseloff and Sergienko2018). 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, Reference Pralong and Funk2005; McGrath and others, Reference McGrath2012b; Luckman and others, Reference Luckman2012; Child and others, Reference Child, Stearns, van der Veen and Elósegui2021). 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, Reference Freed-Brown, Amundson, MacAyeal and Zhang2012). The evolution of basal crevasses has been modeled by the balance of ice shelf and hydrostatic ocean stresses (Zarrinderakht and others, Reference Zarrinderakht, Schoof and Peirce2022), as well as idealized ocean dynamics and the mass balance of melting and freezing (Jordan and others, Reference Jordan, Holland, Jenkins, Piggott and Kimura2014). Basal crevasses may play a central role in the flexure-driven calving of marine-terminating glaciers seen in Greenland and Canada (Wagner and others, Reference Wagner2014; Murray and others, Reference Murray2015; Wagner and others, Reference Wagner, James, Murray and Vella2016; Benn and others, Reference Benn2017). In Antarctica, basal crevasses may initiate rifts (Jeong and others, Reference Jeong, Howat and Bassis2016; Joughin and others, Reference Joughin, Shapero, Smith, Dutrieux and Barham2021) that can propagate across the ice shelf and calve icebergs (Lipovsky, Reference Lipovsky2020). These icebergs can transport fresh meltwater equatorwards and threaten biodiversity of islands in the Southern Ocean (Huth and others, Reference Huth, Adcroft, Sergienko and Khan2022).
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, Reference Lai2020). The depth-averaged deviatoric stress formulations of Nye (Nye, Reference Nye1955), Weertman (Weertman, Reference Weertman1973) and the zero toughness or half-space formulations of Mode I Linear Elastic Fracture Mechanics (LEFM) (van der Veen, Reference van der Veen1998a) 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, Reference Morris and Vaughan2003; van Wessem and others, Reference van Wessem, van den Broeke, Wouters and Lhermitte2023).
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, Reference Goelzer2020; Seroussi and others, Reference Seroussi2020). We consider three theories for Mode I basal crevasses to propagate into rifts: (1) the Zero Stress approximation (Jezek, Reference Jezek1984; Benn and others, Reference Benn, Hulton and Mottram2007; Nick and others, Reference Nick, Van Der Veen, Vieli and Benn2010), (2) LEFM for basal crevasses (van der Veen, Reference van der Veen1998a; Tada and others, Reference Tada, Paris and Irwin2000; Lai and others, Reference Lai2020) and (3) the Horizontal Force Balance (HFB) (Buck, Reference Buck2023). 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, Reference van der Veen1998a) 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, Reference Jezek1984; Benn and others, Reference Benn, Hulton and Mottram2007; Nick and others, Reference Nick, Van Der Veen, Vieli and Benn2010) and Mode I LEFM for basal crevasses (van der Veen, Reference van der Veen1998a) 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 (Reference Buck2023) 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, Reference Walker, Bassis, Fricker and Czerwinski2013) 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, Reference MacAyeal1989) 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, Reference Buck2023). 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, Reference Rist, Sammonds, Oerter and Doake2002; Borstad and others, Reference Borstad2012; Lai and others, Reference Lai2020). We analyze several fracture models (van der Veen, Reference van der Veen1998a; Nick and others, Reference Nick, Van Der Veen, Vieli and Benn2010; Buck, Reference Buck2023) 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, Reference van Wessem2018). To conclude this section, we show that having a Robin temperature profile (Robin, Reference Robin1955) does not change the results of this study.
(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, Reference Jeong, Howat and Bassis2016; Joughin and others, Reference Joughin, Shapero, Smith, Dutrieux and Barham2021) 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, Reference Wuite, Hetzenecker, Nagler and Scheiblauer2021). (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, Reference Haran, Bohlander, Scambos, Painter and Fahnestock2013), grounding lines (Mouginot and others, Reference Mouginot, Scheuchl and Rignot2017b) 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, Reference Nye1955) and LEFM (Smith, Reference Smith1976). However, the use of radar confirmed the presence of crevasses extending upwards from the bottom of the ice (Jezek and others, Reference Jezek, Bentley and Clough1979; Jezek and Bentley, Reference Jezek and Bentley1983). The surface crevasses theories were later extended to study basal crevasses (Jezek, Reference Jezek1984; Benn and others, Reference Benn, Hulton and Mottram2007; Nick and others, Reference Nick, Van Der Veen, Vieli and Benn2010; van der Veen, Reference van der Veen1998a). The Zero Stress approximation we refer to in this paper considers these extensions (Jezek, Reference Jezek1984; Benn and others, Reference Benn, Hulton and Mottram2007; Nick and others, Reference Nick, Van Der Veen, Vieli and Benn2010) 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 (Reference van der Veen1998a) at the approximately vertical basal crevasse interface as
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 (Reference Cuffey and Paterson2010), ice lithostatic pressure and hydrostatic water pressure, respectively. As in the Shallow Shelf Approximation of Stokes flow (MacAyeal, Reference MacAyeal1989), 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−2, vertically integrated ice density ρ i = 917 kg m−3 and ice thickness H. The piezometric head $ z_{\rm h} = H \rho_{\rm i} / \rho_{\rm w} $
as defined in Nick and others (Reference Nick, Van Der Veen, Vieli and Benn2010) depends on the water density; for this study, we assume constant saltwater density ρ w = 1028 kg m−3.
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, Reference Glen1958), the effective viscosity is,
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
Based on lab experiments, the Glen's flow law gives n = 3 (Glen, Reference Glen1958), and B(T) (LeB Hooke, Reference LeB Hooke1981) can be expressed as
where B 0 = 2.207 Pa ⋅ yr1/n, T 0 = 3155 K, T r = 273.39 K, k = 1.17 and C = 0.16612 Kk are constants determined from empirical fit (LeB Hooke, Reference LeB Hooke1981). 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, Reference Nye1955) 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, Reference Jezek1984; Benn and others, Reference Benn, Hulton and Mottram2007; Nick and others, Reference Nick, Van Der Veen, Vieli and Benn2010; Duddu and others, Reference Duddu, Jimènez and Bassis2020). 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, Reference Nye1955; Jezek, Reference Jezek1984; Benn and others, Reference Benn, Hulton and Mottram2007; Nick and others, Reference Nick, Van Der Veen, Vieli and Benn2010; Duddu and others, Reference Duddu, Jimènez and Bassis2020), 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
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 (Reference Weertman1957) 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, Reference De Robin1974; Weertman, Reference Weertman1974), where the spacing between crevasses is much less than the individual crevasse depths (Weertman, Reference Weertman1977). Observations indicate basal crevasse spacing to be roughly one to several ice thicknesses (Luckman and others, Reference Luckman2012; McGrath and others, Reference McGrath2012a; Lawrence and others, Reference Lawrence2023), breaking the densely spaced crevasses assumption. Recent extensions of the Zero Stress approximation include non-zero material strength (Benn and others, Reference Benn, Hulton and Mottram2007) and accumulation and melting effects (Bassis and Ma, Reference Bassis and Ma2015; Huth and others, Reference Huth, Duddu and Smith2021). 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, Reference Buck2023). 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 (Reference Smith1976) and later applied to basal crevasses by van der Veen (Reference van der Veen1998a), considers an isolated basal crevasse with stress concentration near the crack tip and assumes small-scale yielding (Zehnder, Reference Zehnder2012). It has been shown that LEFM agrees with the analytical approach of including stress concentration near crack tips by Weertman (Reference Weertman1973) for small crack depths (Buck and Lai, Reference Buck and Lai2021). Unlike Weertman's infinite thickness assumption (Weertman, Reference Weertman1973), LEFM comes with the advantage of accounting for a prescribed finite thickness (van der Veen, Reference van der Veen1998a). 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 (Reference Zarrinderakht, Schoof and Peirce2022) and Lai and others (Reference Lai2020). Building on previous work (van der Veen, Reference van der Veen1998a; Tada and others, Reference Tada, Paris and Irwin2000; Lai and others, Reference Lai2020), 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 (Reference Tada, Paris and Irwin2000); van der Veen (Reference van der Veen1998a); Lai and others (Reference Lai2020), the stress intensity factor for a Mode I basal crack can be written as
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 (Reference Tada, Paris and Irwin2000) (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, Reference Litwin, Zygielbaum, Polito, Sklar and Collins2012). For rifting, we follow the non-dimensionalization and stability criteria in the Supplemental material of Lai and others (Reference Lai2020): 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 (Reference Zarrinderakht, Schoof and Peirce2022) with isothermal ice to account for the temperature structure effects. As in Appendix C of Zarrinderakht and others (Reference Zarrinderakht, Schoof and Peirce2022), the rifting stress threshold can be determined by zero torque (or moment) into the page associated with a deep basal crevasse. Mathematically,
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 (Reference Zarrinderakht, Schoof and Peirce2022) 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),
with $\tilde {z} = {z}/{H} = 0$
at sea level, constants B 0 = 2.207 Pa ⋅ yr1/n, n = 3, T 0 = 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 (Reference Buck2024). 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,
The numerator of the equation is the same as that in Zarrinderakht and others (Reference Zarrinderakht, Schoof and Peirce2022), 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.
(a) Comparison of depth-averaged ice hardness value $\bar {B}$
normalized by the constant B 0 from LeB Hooke (Reference LeB Hooke1981) 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, Reference Jimènez and Duddu2018; Huth and others, Reference Huth, Duddu and Smith2021; Zarrinderakht and others, Reference Zarrinderakht, Schoof and Peirce2022), 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 (Reference Buck and Lai2021), 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 (Reference Buck2023) 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.
(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, Reference Buck2023), 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, Reference Buck2023). 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, Reference Bassis and Kachuck2023). Taking glaciostatic balance as in Lindstrom and MacAyeal (Reference Lindstrom and MacAyeal1987) as the vertical force balance,
we integrate from a level z to the surface H, i.e. σ zz = ρ ig (z − H), and solve for pressure
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
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 (Reference Rudnicki2014)), 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. the zero material strength assumption of the Zero Stress approximation (Nye, Reference Nye1955; Jezek, Reference Jezek1984; Benn and others, Reference Benn, Hulton and Mottram2007; Nick and others, Reference Nick, Van Der Veen, Vieli and Benn2010),
2. continuity of stress at crack tips (Buck and Lai, Reference Buck and Lai2021) and
3. HFB (Buck, Reference Buck2023).
Previous zero stress models (Nye, Reference Nye1955; Benn and others, Reference Benn, Hulton and Mottram2007; Nick and others, Reference Nick, Van Der Veen, Vieli and Benn2010; Duddu and others, Reference Duddu, Jimènez and Bassis2020) 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),
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
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,
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 (Reference Buck and Lai2021), 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
which easily resolves the constant as
At the basal crevasse tip, we have that
which can readily be simplified to a dimensionless relation between basal and surface crack depth,
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
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
The equation has an analytical solution (Buck, Reference Buck2023),
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$
,
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, Reference Benn, Hulton and Mottram2007; Nick and others, Reference Nick, Van Der Veen, Vieli and Benn2010) and van der Veen (Reference van der Veen1998a)'s basal crevasse LEFM consideration of neglecting time-dependent fracture propagation, but is a possible future extension (Lawn, Reference Lawn1993).
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,
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,
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 (Reference Coffey2022). In the case of rifting with cracks meeting at sea level, this force is
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:
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, Reference Robin1955). 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, Reference Thomas and MacAyeal1982; Rist and others, Reference Rist, Sammonds, Oerter and Doake2002; Craven and others, Reference Craven, Allison, Fricker and Warner2009; Sergienko and others, Reference Sergienko, Goldberg and Little2013; Tyler and others, Reference Tyler2013). 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
with surface temperature T s, thermal diffusivity κ ≡ k/(ρ ic p) ≈ 10−6 m2 s−1 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 (Reference Fowler and Ng2020). We note that the ice divide thickness H d is set to match the Sandhäger and others (Reference Sandhäger, Rack and Jansen2005) 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
Substituting (28) into (27) gives a simple form of the Robin profile,
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.
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, Reference Scambos, Haran, Fahnestock, Painter and Bohlander2007; Haran and others, Reference Haran2018) compared with the rifts previously mapped by Walker and others (Reference Walker, Bassis, Fricker and Czerwinski2013), 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, Reference van Wessem2018) and strain rate (Wearing, Reference Wearing2017) values, with limitations discussed in Appendix A. In Figure 6, we see that on the RIS rifts identified by Walker and others (Reference Walker, Bassis, Fricker and Czerwinski2013), 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 (Reference Zarrinderakht, Schoof and Peirce2022), 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.
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, Reference Haran2018; Scambos and others, Reference Scambos, Haran, Fahnestock, Painter and Bohlander2007), 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 (Reference Walker, Bassis, Fricker and Czerwinski2013) 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 (Reference Zarrinderakht, Schoof and Peirce2022) (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.
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 (Reference Zarrinderakht, Schoof and Peirce2022) 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, Reference Logan, Catania, Lavier and Choi2013; Zarrinderakht and others, Reference Zarrinderakht, Schoof and Peirce2022) 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, Reference Buck2023). 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, Reference Robin1955), 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, Reference Vaughan1995; Bromirski and others, Reference Bromirski, Sergienko and MacAyeal2010; Sergienko, Reference Sergienko2010, Reference Sergienko2013; Hulbe and others, Reference Hulbe2016; Massom and others, Reference Massom2018; Miles and others, Reference Miles2020; Gomez-Fell and others, Reference Gomez-Fell, Rack, Purdie and Marsh2022), with recent work arguing this idea may be most applicable for thin ice (Bassis and others, Reference Bassis2024). 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, Reference Christie2022), 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, Reference Kusahara, Hasumi and Williams2011; Shadwick and others, Reference Shadwick2013; Ohshima and others, Reference Ohshima, Nihashi and Iwamoto2016).
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, Reference Rist, Sammonds, Oerter and Doake2002) 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, Reference Benn2017), ice front bending stresses (Reeh, Reference Reeh1968), grain size dependence of yield stress (Ranganathan and others, Reference Ranganathan, Minchew, Meyer and Peč2021), crack tip shielding (Clayton and others, Reference Clayton, Duddu, Siegert and Martínez-Pañeda2022), snow accumulation, basal melting (Bassis and Ma, Reference Bassis and Ma2015; Kachuck and others, Reference Kachuck, Whitcomb, Bassis, Martin and Price2022; Buck, Reference Buck2023) 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, Reference Gao, Ghosh, Jimènez and Duddu2023).
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, Reference van der Veen1999) and torque from ocean currents (Gomez-Fell and others, Reference Gomez-Fell, Rack, Purdie and Marsh2022; Huth and others, Reference Huth, Adcroft, Sergienko and Khan2022). 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, Reference Walker, Bassis, Fricker and Czerwinski2013; Brunt and others, Reference Brunt, Okal and MacAyeal2011; Bromirski and others, Reference Bromirski2017; Gerstoft and others, Reference Gerstoft2017), viscous creep closure, fatigue failure (Zehnder, Reference Zehnder2012), potentially non-negligible process zones (Zehnder, Reference Zehnder2012), realistic ice rheology or the time-dependent evolution (Lawn, Reference Lawn1993) 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, Reference Weertman1973, Reference Weertman1974; van der Veen, Reference van der Veen1998b; Banwell and others, Reference Banwell, MacAyeal and Sergienko2013; Duddu and others, Reference Duddu, Jimènez and Bassis2020; Lai and others, Reference Lai2020). However, seawater-filled a basal crevasses are predicted to be roughly nine times deeper than dry surface crevasses across several isothermal crevasse theories (Nye, Reference Nye1955; Nick and others, Reference Nick, Van Der Veen, Vieli and Benn2010; Weertman, Reference Weertman1973; van der Veen, Reference van der Veen1998b,Reference van der Veen1998a). This implies that in the absence of strong surface melt (Morris and Vaughan, Reference Morris and Vaughan2003; van Wessem and others, Reference van Wessem, van den Broeke, Wouters and Lhermitte2023) 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, Reference Morlighem2020). Another source of uncertainty comes from the ice rheology. We note that there is uncertainty in the flow law exponent n (Goldsby and Kohlstedt, Reference Goldsby and Kohlstedt2001; Qi and others, Reference Qi, Goldsby and Prior2017; Bons and others, Reference Bons2018; Millstein and others, Reference Millstein, Minchew and Pegler2022) and therefore in the parameters of Eqns (2) and (4) (Zeitz and others, Reference Zeitz, Levermann and Winkelmann2020). 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, Reference Huth, Duddu and Smith2021), capturing the path of rift propagation to form iceberg A68 from Larsen C (Huth and others, Reference Huth, Duddu, Smith and Sergienko2023). 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 (Reference Bassis and Walker2012); Pollard and others (Reference Pollard, DeConto and Alley2015); Bassis and Ma (Reference Bassis and Ma2015). For example, because the HFB model predicts much deeper cracks compared with the Zero Stress approximation that was used in Pollard and others (Reference Pollard, DeConto and Alley2015), implementing our HFB fracture model instead could yield substantially different outcomes, making ice shelves more vulnerable to fracture than predicted by Pollard and others (Reference Pollard, DeConto and Alley2015).
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, Reference Olinger, Lipovsky and Denolle2024). 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.
Acknowledgements
We acknowledge funding from NSF's Office of Polar Programs through OPP-2235051. C.-Y.L. acknowledges NASA for partial support via Cryosphere Award 80NSSC21K1003.
Appendix A. Methods for the data-model comparison
To create rift prediction maps, several important steps must be made to ensure sensible and causal predictions. First, we must validate that the areas containing rifts are in regions that largely obey the 1-dimensional (1D) extensional background flow assumptions of the fracture theories assessed in this study. To do so, we use automatic differentiation to construct strain rate fields based on MEaSUREs ice shelf velocity data (Rignot and others, Reference Rignot, Mouginot and Scheuchl2011; Mouginot and others, Reference Mouginot, Scheuchl and Rignot2012, Reference Mouginot, Rignot, Scheuchl and Millan2017a; Rignot and others, Reference Rignot, Mouginot and Scheuchl2017).
To determine if the 1D fracture theory assumptions are upheld, we utilize the criterion that the normalized resistive stress difference between that of the SSA solution and that assuming 1D flow is within $10\%$
,
where $\bar {R}_{xx}^{{\rm SSA}} = \bar {B} \dot {\epsilon }_{\rm e}^{-1 + ( {1}/{n}) } ( 2 \dot {\epsilon }_{xx} + \dot {\epsilon }_{yy} )$
and $\bar {R}_{xx}^{1{\rm D}} = 2 \bar {B} \dot {\epsilon }_{xx}^{{1}/{n}}$
. This equation can be simplified and written in terms of dimensionless measures of strain rate,
where the first term with $\alpha = {\dot {\epsilon }_{yy}}/{\dot {\epsilon }_{xx}}$
and $\xi = {\dot {\epsilon }_{xy}}/{\dot {\epsilon }_{xx}}$
comes from the second invariant of the strain rate tensor under the 3D incompressible condition and assumption of negligible vertical shear stresses. Thus, while one could calculate the depth-averaged ice hardness $\bar {B}( T)$
assuming some vertical temperature profile, this is unnecessary as this term cancels and does not appear in Eqn (A2). Figure 8 shows the region which satisfies the criterion (A2), from which we may select our rifts.
Regions where the background viscous ice flow is approximately 1D based on Eqn (A2). Rifts in these regions are regarded as having formed due to 1D tension, called Mode I failure. The color scale is dimensionless strain rate deviation from 1D flow as in Eqn (A2), and the axes are in units of meters.
Second, the rifts identified by Walker and others (Reference Walker, Bassis, Fricker and Czerwinski2013) in these regions must have their thickness values padded to an estimate of the unbroken state thickness. These rifts are currently filled with a conglomerate ice material composed of sea ice, snow and ice shelf fragments, collectively termed as ice mélange in Rignot and MacAyeal (Reference Rignot and MacAyeal1998); MacAyeal and others (Reference MacAyeal, Rignot and Hulbe1998); Hulbe and others (Reference Hulbe, Rignot and Macayeal1998). Since our goal is to predict if the observed rift could have formed, rather than the stability of the mélange in the rift, we need to estimate the state of the unbroken ice. Here, we generate bounding boxes around rifts of interest and infill the mélange thickness with an average of the local unbroken ice thickness in BedMachine Version 2 (Morlighem, Reference Morlighem2020; Morlighem and others, Reference Morlighem2020), as shown in Figure 9. Combined with the map of regions that uphold the 1D extensional flow assumption of our fracture theories in Figure 8, we generate rift prediction maps in Figures 6 and 10.
Mélange padding in kilometers with local, unbroken ice thickness on the RIS and LCIS. Subfigures (a) and (c) are the original data products from BedMachine2, while subfigures (b) and (d) are the mélange-padded results used in this study to see if fracture theories can correctly predict rifts in areas where they are known to have occurred. LCIS thickness has an upper bound of 300 m in the color bars of subfigures (c) and (d) to enhance the visibility of padded rifts.
Same as Figure 6, except on LCIS, assuming a linear temperature profile.
The caveat to this method is that in the absence of remotely sensed data on the time-dependent evolution of basal crevasses into rifts, strain rates and temperature profiles at the time of the basal crack to rift transition remain unknown. We present an analysis of time series data in Appendix B, but these data are largely unavailable for many existing rifts. As such, we utilize the modern values of surface temperatures and strain rates from van Wessem and others (Reference van Wessem2018) and Wearing (Reference Wearing2017), acknowledging that these may have changed since rift formation. We note that the 1/n exponent dependence of strain rate given Glen's flow law in Eqn (3) shields the resistive stress from strain rate changes, thus decreasing the sensitivity to precise strain rate estimates.
Due to the uncertainty associated with the temperature and strain rate evolution since the time of rift formation, we construct datasets of stress estimates on the non-rift ice shelf regions, as shown in Figure 11. The advantage of these datasets is that they do not have the time-evolution problem of rift datasets. Therefore, given our methods, we have more confidence in a theory that correctly predicts many rifts and minimally overpredict rifts than a theory which correctly predicts most or all rifts while overpredicting rifts in non-rift regions. In our work, this emphasizes that LEFM with depth-averaged resistive stress result in Zarrinderakht and others (Reference Zarrinderakht, Schoof and Peirce2022) has too low of a threshold for rift initiation.
Extensional, approximately 1D regions of RIS and LCIS, assuming a linear temperature profile, that exclude both observed rifts and rifts predicted by a resistive stress greater than twice the freely floating resistive stress. The color shows the magnitude of $\bar {R}_{xx} / \bar {R}_{xx}^{IT}$
, and axes are in meters.
We note that since Walker and others (Reference Walker, Bassis, Fricker and Czerwinski2013) do not necessarily identify all rifts in their rift catalog, there could be rifts included in our non-rift region datasets. To help remove rifts not classified by Walker and others (Reference Walker, Bassis, Fricker and Czerwinski2013) in the non-rift regions of Figure 11, we exclude data that would unanimously be predicted a rift in every fracture theory presented in this study. As such, we take the highest stress threshold for rifting in the Zero Stress approximation with isothermal depth-averaged resistive stress, $\bar {R}_{xx}^\ast /\bar {R}_{xx}^{IT} = 2$
. This is reflected in the magnitude of the color bar of Figure 11.
Appendix B. Temporal observation of a basal crevasse to rift transition
Figures 1c,e show the ratio of depth-averaged resistive stress to isothermal ice tongue resistive stress $\bar {R}_{xx}/R_{xx}^{IT}$
over Pine Island Ice Shelf in January and May 2019 respectively. The resistive stress was found using Eqn (7) with along-flow strain rates estimated from 200 × 200 m resolution, monthly averaged ice velocity observations made using feature tracking applied to Sentinel-1 image pairs (Wuite and others, Reference Wuite, Hetzenecker, Nagler and Scheiblauer2021) (https://cryoportal.enveo.at/data/). The strain rate components were calculated via numerical differentiation of the easting and northing velocities using methods developed in Chartrand (Reference Chartrand2017). Estimates of ice shelf thickness were according to BedMachine version 2 by Morlighem (Reference Morlighem2020), and the temperature profile was assumed to be linear between −2 °C at the ice shelf base and temperature given by RACMO (van Wessem and others, Reference van Wessem2018) at the surface. The accompanying satellite images shown in (d) and (f) are geocoded, multi-looked and radiometrically terrain-corrected Single-Look Complex backscatter data from the European Space Agency and European Commission Copernicus’ Sentinel-1 satellites – shown at 50 m resolution.
These show the concurrent evolution of the ratio of stresses near the terminus of the ice shelf alongside the evolution of a rift, likely from a central basal crevasse as argued by Jeong and others (Reference Jeong, Howat and Bassis2016), that eventually led to the calving of the B49 iceberg in February 2020. We see clear changes to the along-flow strain rates over the rift as it widens and propagates laterally. This provides a motivating example for determining stress conditions under which basal crevasses transition into rifts, such as those discussed in this article. However, Figures 1c,e show that we measure stress ratios below those required for full-thickness rifts according to each of the theories discussed in this article. That is, the stress ratio over the rift that is clearly visible by May 2019 (Fig. 1f) is below the value of 1 predicted by the Zero Stress approximation, modified to maintain horizontal force balance (Eqn (24)). However, Figure 1 should not be seen as an example aimed at validation of one of the theories considered here, merely as motivation for the work. In part, this is because the central part of Pine Island Ice Shelf does not conform to the assumption of one-dimensional flow (Eqn (A2)) (though the speed of the ice tongue varies little laterally, the flow is dominated by advection and across-flow strain rates are similar in magnitude to along-flow strain rates), nor is the ice shelf cavity necessarily hydrostatic with temperature T = −2 °C. Additionally, by using satellite-derived measures of ice velocity averaged over monthly intervals, we cannot hope to capture the maximum strain rates over a crevasse of this scale as the data are too limited in both spatial and temporal resolution. In the future, a validation of the theories discussed in this article using time series of strain rate or stress data should be carried out with the use of higher-resolution satellite data or data collected on the ground, e.g. with the use of an ApRES system as in Nicholls and others (Reference Nicholls2015). Similarly, it is not possible to accurately determine when the crevasse transitioned into a rift from satellite images alone, or whether it was ever a basal crevasse at all. Future work that aims to use time series data should do so in conjunction with other datasets that provide further information on the type of crevasse under consideration.
Appendix C. Zero stress for rift formation via surface and basal crevasses
The Zero Stress approximation argues that a vertical crack will propagate so long as there is no net compression of the net longitudinal stress σ n at the crack tip defined in Eqn (1). Written mathematically, the Zero Stress condition (σ c = 0) claims that a crack propagates when
Under the Zero Stress approximation, a basal and surface crevasse will form a rift when the criterion (C1) holds for all depths.
The depth of surface and basal crevasses that result from this Zero Stress condition (C1) are available in the literature (Nye, Reference Nye1955; Jezek, Reference Jezek1984; Benn and others, Reference Benn, Hulton and Mottram2007; Nick and others, Reference Nick, Van Der Veen, Vieli and Benn2010; Duddu and others, Reference Duddu, Jimènez and Bassis2020),
The criterion (C1) can also be re-written in terms of a dimensionless resistive stress, with z = 0 at the bottom of the ice,
as visualized by the red dotted lines in Figure 12, which shows basal crevasse depth and rift formation with a linear temperature profile modifying the resistive stress R xx(z).
Visualizing the Zero Stress condition. The dotted red curve, defined by the piecewise distribution in Eqn (C3), is the magnitude of the lithostatic pressure minus water pressure that opposes fracture; the intersection below sea level of a given resistive stress with the dotted red curve determines basal crevasse depth. Solid curves utilize R xx(z), whereas dashed lines utilize $\bar {R}_{xx}$
. Green curves with surface temperature T s = −22 °C are the reference for cyan and magenta curves; cyan curves utilize T s = −32 °C, and magenta curves have along-flow strain rate doubled. Depth-averaged resistive stresses have unstable basal crevassing occur only at sea level, and may be solved for analytically; vertical temperature profiles may have crevasses unstably propagate before sea level, and require numerical treatment.
Since the Zero Stress approximation has zero material strength, the minimum required resistive stress to form a rift is defined by Eqn (C3). If the resistive stress is not in net tension at the ice shelf base, no basal crevasse is predicted. If the resistive stress is in net tension at the base but becomes less than the dotted red curve in Figure 12 at a larger height z > 0, the point of equality below sea level is the basal crevasse depth. For example, if the resistive stress takes the value shown by the dashed green line of Figure 12, a basal crevasse would propagate up to a depth about $60\%$
of the unbroken ice thickness. However, resistive stresses that are greater than or equal to the dotted red curve for all heights will form rifts because basal crevasses can propagate all the way to the surface. Figure 12 clearly demonstrates the underestimation of rifts when depth-averaged resistive stress theories, the dashed lines, are used instead of their depth-dependent counterparts, the solid curves.
Next, we develop the mathematical expression for the rift initiation stress threshold of the Zero Stress approximation. These expressions are plotted in dashed blue and solid green in Figures 7 and 16. Taking the assumption that the second invariant of strain rate is approximately the along-flow strain rate for consistency with LEFM, the rift formation criterion given isothermal, depth-averaged resistive stress can be written as a dimensionless stress ratio,
with $\bar {R}_{xx}^{IT} = ( 1- {\rho _{\rm i}}/{\rho _{\rm w}} ) \rho _{\rm i} g H/2 $
the depth-averaged ice tongue resistive stress. This equation can be understood visually from Figure 12 as the corner of the dotted red curve located on the x-axis at $R_{xx} / \left ( \left ( \rho_w - \rho_i \right ) g H \right ) = \rho_i /\! \rho_w \approx 0.89 $
, which is equivalent to (C4). Similarly, in the depth-dependent case, we have that
Here $\tilde {B}( T) = B( T) /B( T = -2$
°C) is the dimensionless ice hardness, $\bar {\tilde {B}}( T )$
is the depth-averaged dimensionless ice hardness, and $d_{\rm b}^\ast$
is the unstable basal crevasse depth at which a basal crevasse will propagate to form a rift. The unstable basal crevasse depth $d_{\rm b}^\ast$
depends upon temperature profile and the prescribed stresses. For isothermal ice, the unstable basal crevasse height is sea level without tides, $ d_b^* = H \rho_i / \rho_w $
, and we also have that $\bar {\tilde {B}}( T ) = \tilde {B}( T( {d_{\rm b}^\ast }/{H}) )$
, so the above Eqn (C5) reduces to the depth-averaged case in Eqn (C4). For vertical temperature profiles that become colder toward the ice shelf surface, the unstable basal crevasse depth $d_{\rm b}^\ast$
can decrease. Given the ice hardness function of LeB Hooke (Reference LeB Hooke1981) and a linear temperature profile from T b = −2 °C at the base, the unstable basal crevasse depth $d_{\rm b}^\ast$
falls below sea level for surface temperatures at least as cold as T s = −25 °C in the Zero Stress theory as well as HFB, as shown with the blue curves of Figure 4a with T s = −32 °C.
Appendix D. Temperature dependence of linear elastic fracture mechanics for ice shelf basal crevasses
In this Appendix section, we first explain the stability of LEFM basal crevasses given (1) and (6), and then the effects of ice-shelf temperature on crack depth and rifting stress, as shown in Figure 4c.
D.1. LEFM basal crack stability
Figure 13 indicates the possible solutions for LEFM through increasing the resistive stress, similar to presentation in Lai and others (Reference Lai2020). In the blue curve, we see that the stress intensity factor is less than the fracture toughness for all crack depths, implying no fracture could exist for that given resistive stress. For a larger resistive stress (red curve), a tangent solution exists where the stress intensity factor is equal to the fracture toughness, so a stable crack exists, but it cannot propagate. This crack depth is the maximum initial flaw size required to form a stable crack. Following (Lai and others, Reference Lai2020), the maximum initial flaw size is κ(K Ic/(ρ ig))2/3 ≈ 15 m for a basal crevasse, with κ ≈ 2.27, K Ic = 0.15 MPa⋅m1/2 from Litwin and others (Reference Litwin, Zygielbaum, Polito, Sklar and Collins2012), ρ i = 917 kg m−3 and g = 9.8 m s−2. For an even larger resistive stress, an initial flaw (open circle) would propagate as long as K I ≥ K Ic, potentially finding a stable crack depth, as shown by the solid circle on the yellow curve. Finally, there exists a rifting stress $\bar {R}_{xx}^\ast$
above which initial flaws would unstably propagate through the entire ice thickness, as shown by the purple curve.
Understanding possible LEFM solutions for isothermal ice with a constant thickness H = 300 m and fracture toughness K Ic = 150 kPa$\sqrt {\hbox{m}}$
(dashed line). The solid curves are the stress intensity factors from (6) with varying resistive stress. Descriptions of the possible solutions are in the text of this Appendix D.
Below, we discuss the dependence of LEFM results on vertical temperature profile in terms of (1) stable crevasse depths and (2) the rifting stress threshold $\bar {R}_{xx}^\ast$
that corresponds to the first unstable solution.
D.2. Temperature effects on basal crevasse depths and the rifting stress threshold
As shown in our Figure 4c, for the same depth-averaged resistive stress, colder ice has shallower crack depths. This may seem at odds with Figure 3 of van der Veen (Reference van der Veen1998a) that shows colder ice temperatures promote deeper crack depths. However, it is actually consistent; Figure 3 of van der Veen (Reference van der Veen1998a) is for a constant strain rate, whereas our Figure 4 is for constant normalized resistive stress. Since the resistive stress depends exponentially on the temperature through the effective viscosity as shown in (2), (3) and (4), colder temperatures will have higher viscosity and thus decreased strain rate for the same value of depth-averaged resistive stress. Thus, for the same depth-averaged resistive stress, colder ice has shallower crack depths.
As shown by the LEFM solutions in Figure 4c that turn toward vertical as d b goes to H; there is an unstable basal crevasse depth $d_{\rm b}^\ast$
at which rifts are formed, similar to the other theories in this study. Since $d_{\rm b}^\ast < H$
is when rifting is determined, and $\int _0^{d_{\rm b}} ( \bar {R}_{xx}-R_{xx} ( z) ) \rm{d}z$
increases for colder surface temperatures given a monotonically decreasing temperature profile toward the ice surface, the diverging rifting stress between the depth-averaged resistive stress $\bar {R}_{xx}$
and vertically varying resistive stress R xx(z) is expected, as shown in Figure 7.
Appendix E. Derivation of horizontal force balance
Here we demonstrate how the Zero Stress approximation does not uphold horizontal force balance on an isothermal ice shelf through an Eulerian control volume argument based on Buck (Reference Buck2023), also see Section 2.3. The main argument of the control volume approach, as has been applied by Weertman (Reference Weertman1957) and Jezek (Reference Jezek1984) to solve for the net tension we call $\bar {R}_{xx}^{IT}$
at ice fronts, is Newton's Second Law. The sum of the forces acting on the control volume are equal to the product of mass and acceleration of fluid entering the control volume. In our case, there is no net acceleration of fluid into or out of the control volume, and the shear stresses on surface and bottom boundaries are negligible. Thus, we can write the horizontal force balance for a control volume between a crevassed location x = x c and an uncrevassed downstream location x = x c + Δx as
The horizontal force balance model for an isothermal ice shelf was developed in Buck (Reference Buck2023) and is summarized below. At the downstream location x = x c + Δx that is sufficiently far away from the bending stresses near the ice front (Reeh, Reference Reeh1968; Wagner and others, Reference Wagner, James, Murray and Vella2016), we have
At the crevassed location, we will follow the Zero Stress assumption and have dual surface and basal crevasses with depths d s and d b,
Note that the stresses cannot be the same at both locations, or σ xx(x = x c) ≠ σ xx(x = x c + Δx), because the intact ice is effectively thinner at the crevassed location. If we were to evaluate the force balance of Eqn (E1) with the incorrect assumption of σ xx(x = x c) = σ xx(x = x c + Δx), the crack depths would be twice as deep as that of the Zero Stress theory,
and the stress distribution at the surface crevasse tip would not be continuous, σ xx(z = H − d s) = −ρ ig d s + R xx ≠ 0. Thus, to satisfy a continuous stress at the surface crack tip, we have
Similarly, using stress continuity at the basal crack tip gives a relation between surface and basal crack depths,
With the crack depth relation in Eqn (E6), plugging the stress definitions in Eqns (E2), (E3), (E5) into the force balance condition of Eqn (E1) yields the analytical crack depth predictions of Buck (Reference Buck2023) for an isothermal ice shelf,
For more insight into the role of temperature dependence, we now specify the form of Eqns (19) and (20) for a simplified, approximate ice hardness function and linear vertical temperature profile. In Eqn (4), the second term in the brackets, −C/(T r − T)k, is two to three orders of magnitude smaller than the first term, T 0/T. Similarly, with temperature $T( \tilde {z}) = T_{\rm b} [ 1 - ( 1 - T_{\rm s}/T_{\rm b}) \tilde {z} ]$
in Kelvin, the term $( 1 - T_{\rm s}/T_{\rm b} ) \tilde {z}$
is at least an order of magnitude smaller than unity, so we may Taylor expand the exponent to first order in $( 1 - T_{\rm s}/T_{\rm b} ) \tilde {z}$
. We define the approximated ice hardness function B a with these two simplifications,
with the dimensionless e-folding length scale $\tilde {z}_0 \equiv ( ( {T_0}/{T_{\rm b}}) ( 1- {T_{\rm s}}/{T_{\rm b}}) ) ^{-1}$
. Therefore, the crevasse depth relation of Eqn (19) may be written as
and the horizontal force balance of Eqn (20) may be written as
The same map view plot in kilometers as Figure 6, except that a Robin temperature profile is taken to generate the rift formation stresses and assumed in calculating the value of stress in the data.
Even with the simplified ice hardness, these equations (E9) and (E10) are not algebraically solvable due to the nature of the Arrhenius equation. Although the result including vertically varying temperature requires numerical treatment, the rift initiation stress threshold produced using LeB Hooke (Reference LeB Hooke1981)'s ice hardness function is within $0.1\%$
of the analytical isothermal solution ${R_{xx}^\ast }/{\bar {R}_{xx}^{IT}} = 1$
for all surface temperatures used for both linear and Robin temperature profiles. Therefore, we can well-approximate the rift initiation stress threshold as that of a freely floating ice shelf without buttressing, i.e.
Appendix F. Result robustness: uncertainty estimation
A discussion of result robustness is incomplete without considering the uncertainty in data products. The largest data uncertainty comes from the measurements of ice thickness (Morlighem, Reference Morlighem2020; Morlighem and others, Reference Morlighem2020), where the uncertainties in our regions of interest are 100 m for the majority of the RIS, or around a third of the ice thickness, and around 30 m for the LCIS, or about a tenth of the ice thickness. Alone, one standard deviation of this uncertainty would shift the data points up or down by about a third for RIS data or about a tenth for LCIS data on Figures 7 and 16. As such, we look at LCIS for result robustness. Importantly, if the ice shelf data of interest is governed by the 1D SSA momentum equation (MacAyeal, Reference MacAyeal1989),
then the depth-averaged resistive stress scales linearly with H. To compute uncertainty accurately for dependent variables, we would have to use covariance (Taylor, Reference Taylor1982); however, we cannot meaningfully compute the covariance for each pixel of ice shelf data, and so we estimate the upper bound on uncertainty σ β with
Here, our variable of interest is the dimensionless resistive stress $\beta = \bar {R}_{xx}/\bar {R}_{xx}^{IT}$
. The uncertainties in thickness, strain rate and equivalent temperature are $\sigma _{H},\; \, \sigma _{\dot {\epsilon }_{xx}}$
and σ T*, with equivalent temperature T* defined as the temperature at which $\bar {B} = B( T^\ast )$
(Sergienko, Reference Sergienko2014). We take the strain rate uncertainty associated with 20 km from the ice front from Table C.1 of Wearing (Reference Wearing2017) and apply this to the whole ice shelf. Given that we do not have defined uncertainties associated with equivalent temperature, we estimate σ T* = 3 K from the uncertainty range associated with modeled and observed RACMO surface temperature data in Figure 3a of van den Broeke (Reference van den Broeke2008). In this calculation, we assume the Robin temperature profile in our ice hardness and equivalent temperature calculations, as we do not expect profiles warmer than linear, but this choice is negligible in the final results.
The same idea as Figure 7, except that we include approximate rift stress data and evaluate linear and Robin temperature profiles for generating the depth-dependent rift formation stresses (green, cyan and red curves). Relative to the linear temperature profile, the Robin profile raises the value of the resistive stresses in the rifts and non-rift ice shelf datasets (black and orange data, respectively), but has a negligible effect on the rifting stress threshold curves. Horizontal Force Balance is a much more accurate rift initiation theory compared with the Zero Stress approximation, and is largely or fully insensitive to surface temperature assuming a linear or Robin profile. This provides robustness, as the same conclusions are drawn with either a linear or Robin temperature profile.
The same as Figure 11, except that the Robin temperature profile raises the value of the depth-averaged resistive stress relative to the linear temperature profile. Axes are in meters.
We plot the upper bound of dimensionless resistive stress uncertainty σ β in Figure 18. Given the distributions of these datasets have some large outliers that skew the mean, we report the estimated median uncertainties for RIS and LCIS are σ β = 0.27 and σ β = 0.14, respectively. The RIS median uncertainty is large as anticipated, and the LCIS median uncertainty is comparable to the difference between LEFM with R xx(z) and Horizontal Force Balance given the temperatures on LCIS (see red and cyan curves in Figs 7, 16). Therefore, more precise measurements of ice thickness, strain rate and temperature are needed to further observationally constrain the optimal theory for tensile rift initiation from basal crevasses.
Estimates of dimensionless resistive stress uncertainty σ β defined in Eqn (F2) on RIS and LCIS, with thickness padded for known rifts. The small uncertainty on LCIS provides more confidence in our work. The color scale is capped at σ β = 1, as there are a few outliers due to thickness uncertainty being comparable or larger than ice thickness. Axes are in meters.
Open access
