1. Introduction
Recent reviews of ice calving and stability (Alley and others, Reference Alley2023; Bassis and others, Reference Bassis2024) elucidate how glacial retreat can transition across different calving regimes, e.g. from ice shelves to deep water glaciers to shallow water glaciers (Fig. 1). Assessing the impact of the marine ice cliff instability (MICI) (Bassis and Walker, Reference Bassis and Walker2012; Pollard and others, Reference Pollard, DeConto and Alley2015; DeConto and Pollard, Reference DeConto and Pollard2016) on sea level rise via ice sheet models is challenging due to the incomplete parametrization of calving rates (Morlighem and others, Reference Morlighem2024) that are poorly constrained by observations. The recent Intergovernmental Panel on Climate Change assessment shows that MICI contributes to a highly uncertain high-end scenario, which can result in about one meter of global sea-level rise by 2100 (Pollard and others, Reference Pollard, DeConto and Alley2015; DeConto and Pollard, Reference DeConto and Pollard2016; Fox-Kemper and others, Reference Fox-Kemper2021).
Figure 1. Schematic of the three cases considered in this paper: crevasses on an ice shelf, a marine-terminating glacier and a land-terminating glacier. Surface crevasses are considered in all three cases. Basal crevasse depth depends on whether it is filled with ocean saltwater or subglacial meltwater. Calving occurs when the surface crack and basal crack depths (
$d_s, d_b$) fully occupy ice’s thickness, 
$d_s+d_b=H$ which gives various calving laws derived in this paper.
Despite its importance, some of the most fundamental questions surrounding ice crevassing remain unanswered, including a predictive calving criterion. Various calving laws have been developed to model either the retreat rate (e.g. the von Mises law (Morlighem and others, Reference Morlighem2016), eigencalving law (Levermann and others, Reference Levermann, Albrecht, Winkelmann, Martin, Haseloff and Joughin2012)) or position (e.g. the crevasse-depth law (Nye, Reference Nye1955; Benn and others, Reference Benn, Hulton and Mottram2007a; Nick and others, Reference Nick, van der Veen, Vieli and Benn2010)) of the calving front. Despite its wide use, it has been shown (Buck, Reference Buck2023; Coffey and others, Reference Coffey, Lai, Wang, Buck, Surawy-Stepney and Hogg2024) that in the Zero-Stress approximation (Nye, Reference Nye1955; Benn and others, Reference Benn, Hulton and Mottram2007a; Nick and others, Reference Nick, van der Veen, Vieli and Benn2010), the depth-integrated force at the crevassed and non-crevassed location are unbalanced. This has led to a modified crevasse-depth model for constant-thickness ice shelves (ISs) that satisfies Horizontal Force-Balance (HFB) (Buck, Reference Buck2023). Importantly, HFB analytically predicts that the tensile stress required for calving is only half of that in the Zero-Stress approximation, which may substantially underestimate a glacier’s vulnerability to calving. HFB provides reasonable agreement with observed ice shelf rift locations and yields a calving threshold that is insensitive to the vertical temperature profile (Coffey and others, Reference Coffey, Lai, Wang, Buck, Surawy-Stepney and Hogg2024). Since analytical theories serve as foundational cases for developing physical understanding and benchmarking numerical methods that can simulate more complicated phenomena, it is crucial to develop physically self-consistent fracture models that can apply across diverse glacial environments, from ISs and marine-terminating glaciers (MTGs) to land-terminating glaciers (LTGs). This is the goal of the paper.
In this paper, we generalize the HFB approach of Buck Reference Buck(2023) and Coffey and others Reference Coffey, Lai, Wang, Buck, Surawy-Stepney and Hogg(2024) from constant-thickness ISs to land- and marine-terminating glaciers (Fig. 1). We consider six “crack configurations” involving different combinations of dry or meltwater surface and saltwater- or meltwater-filled basal crevasses (Fig. 3). Section 2 lays out the general formulation of HFB and the driving and resisting mechanisms for calving. Applying the general HFB formulation gives the analytical crack depths and calving criteria predictions for ISs, MTGs and LTGs in sections 2.1, 2.2 and 2.3, respectively. We hypothesize that the variety of observed calving styles may arise from different dominant balances in our force-balance equation, between sources of buttressing and drivers of calving, as discussed in section 3.
2. A force-balance framework (HFB) to predict crevasse depths
We begin this section by outlining the general steps to use the HFB approach to determine crevasse depths, depicted in Fig. 2. This approach is then applied to a variety of cases, from ISs (section 2.1, Fig. 5) to LTGs (section 2.2, Fig. 7) and MTGs (section 2.3, Fig. 2).
Figure 2. Schematic of the forces that drive calving (red) or inhibit calving (green), with force balance conceptualized atop the cartoon. See Table 1 for descriptions of symbols. Throughout this paper, saltwater and meltwater are shown in blue and light green, respectively.
Assuming no acceleration, all of the forces F acting on a body must balance such that their sum is zero, 
$\Sigma F = 0$. This paper focuses on using the 
$\Sigma F = 0$ constraint to determine the crack depths and calving criteria.
The force balance 
$\Sigma F = 0$ can be expressed with a volume integral of the Stokes equation within a control volume V,
\begin{equation}
\iiint \left[\underline{\nabla} \cdot \underline{\underline{\sigma}} + \rho_i \underline{g} \right]dV = \underline{0}.
\end{equation} In (1), 
$\underline{\nabla}$ is the vector differential operator, 
$\underline{\underline{\sigma}}$ is the ice shelf Cauchy stress tensor, ρi is the ice density and 
$\underline{g} $ is the gravitational acceleration vector. Applying Gauss’s theorem to (1), we obtain
\begin{equation}
\unicode{x222F} \underline{\underline{\sigma}} \cdot \underline{n} d S + \iiint \rho_i \underline{g} dV = \underline{0}.
\end{equation} The horizontal component (perpendicular to gravity) of (2) per unit width (into the page) can be written in terms of the surface forces or traction 
$\underline{t} = \underline{\underline{\sigma}} \cdot \underline{n}$ (e.g. Malvern, Reference Malvern1969; Dahlen and Tromp, Reference Dahlen and Tromp1998; Rudnicki, Reference Rudnicki2014),
\begin{equation}
\left ( \underline{\hat{x}} \right ) : \unicode{x222F} \underline{\hat{x}} \cdot \underline{\underline{\sigma}} \cdot \underline{n} d S = \unicode{x222F} \underline{\hat{x}} \cdot\underline{t} d S = \Sigma F_x = 0.
\end{equation} This equation is the HFB per unit width. All that we need to know is the traction along the boundaries of our domain, or the sum of the forces acting on our body. This idea of integrating the momentum equation over a control volume is commonly used in many problems, such as with the Betti reciprocal relation (e.g. Dahlen and Tromp, Reference Dahlen and Tromp1998) and the Boundary Element Method (e.g. Crouch, Reference Crouch1976; Crouch and Starfield, Reference Crouch and Starfield1983; Zarrinderakht and others, Reference Zarrinderakht, Schoof and Peirce2022). We use the glaciology convention (e.g. van der Veen and Whillans, Reference van der Veen and Whillans1989; Cuffey and Paterson, Reference Cuffey and Paterson2010) and define the resistive stress 
$\underline{\underline{R}}$ such that
\begin{equation}
\underline{\underline{\sigma}} = -p_l \underline{\underline{I}}+\underline{\underline{R}},
\end{equation} or written in terms of the 
$j, k$ components of the stress tensor, 
$\sigma_{jk} = -\delta_{jk}p_l \ + \ R_{jk}$, where 
$p_l\equiv \rho_i g (s-z)$ is the lithostatic pressure (scalar), δjk is the Kronecker delta function, s is the surface elevation and the vertical coordiante z is zero at the sea level, increasing upward.
We now consider two states: unfractured and fractured states (e.g. before and after crevassing). In the unfractured state, we may write (3) for a control volume delineated by the yellow dashed lines in Fig. 2 as
\begin{align}
&\int_{y_L}^{y_R}\int_{b}^s \underline{\hat{x}} \cdot \underline{\underline{\sigma}}^0 \left(x=x_c\right)\cdot \underline{n} dz dy= \iint\limits^{\text{Surface}} \underline{\hat{x}} \cdot \underline{\underline{\sigma}}^0 \cdot \underline{n} d S+ \nonumber \\ & \iint\limits^{\text{Base}} \underline{\hat{x}} \cdot \underline{\underline{\sigma}}^0 \cdot \underline{n} d S+\iint\limits^{\text{Front}} \underline{\hat{x}} \cdot \underline{\underline{\sigma}}^0 \cdot \underline{n} d S + \iint\limits^{\text{Left}} \underline{\hat{x}} \cdot \underline{\underline{\sigma}}^0 \cdot \underline{n} d S + \nonumber \\ & \iint\limits^{\text{Right}} \underline{\hat{x}} \cdot \underline{\underline{\sigma}}^0 \cdot \underline{n} d S ,
\end{align} where the superscript 0 denote the stress without crack formation. In the fractured state at 
$x=x_c$, the stresses at the cracked location differ from the unfractured state and (3) becomes
\begin{align}
&\int_{y_L}^{y_R}\int_{b}^s \underline{\hat{x}} \cdot \underline{\underline{\sigma}} \left(x=x_c\right)\cdot \underline{n} dz dy= \iint\limits^{\text{Surface}} \underline{\hat{x}} \cdot \underline{\underline{\sigma}} \cdot \underline{n} d S+ \nonumber \\ & \iint\limits^{\text{Base}} \underline{\hat{x}} \cdot \underline{\underline{\sigma}} \cdot \underline{n} d S+\iint\limits^{\text{Front}} \underline{\hat{x}} \cdot \underline{\underline{\sigma}} \cdot \underline{n} d S +\iint\limits^{\text{Left}} \underline{\hat{x}} \cdot \underline{\underline{\sigma}} \cdot \underline{n} d S +\nonumber \\ & \iint\limits^{\text{Right}} \underline{\hat{x}} \cdot \underline{\underline{\sigma}} \cdot \underline{n} d S.
\end{align} Since the terms on the right-hand side of both equations are associated with horizontal forces acting on the control volume boundaries excluding the crevassing plane, which do not change between the fractured and the unfractured states, the right-hand sides of both equations are the same. Thus, with negligible vertical shear stress or approximately vertical cracks at 
$x=x_c$, the HFB for a control volume can be reduced to a local equation
\begin{equation}
\int_{y_L}^{y_R}\int_{b}^s \sigma_{xx}^0 \left(x=x_c\right) dz dy = \int_{y_L}^{y_R}\int_{b}^s \sigma_{xx}\left(x=x_c\right) dz dy.
\end{equation}For convenience, the stresses are width-averaged in the rest of the paper. The local HFB equation thus becomes
\begin{equation}
\int_{b}^s \sigma_{xx}^0 \left(x=x_c\right) dz = \int_{b}^s \sigma_{xx}\left(x=x_c\right) dz.
\end{equation}The HFBs in (3) and (8) applied to glacier crevassing have been extensively discussed in Buck Reference Buck(2023) and Coffey and others Reference Coffey, Lai, Wang, Buck, Surawy-Stepney and Hogg(2024). In this paper, we will use the local force balance (8) for each case considered in Fig. 1, from ISs through MTGs to LTGs.
Figure 2 shows an example of the forces acting on a control volume aligned with tensile, vertical crevasses on a MTG that can be described with (3). The driver of calving (red arrows) in our static framework is the hydrostatic water pressure in surface and basal crevasses. The forces that inhibit calving (green arrows), via reducing the tensile glaciological stresses around the crevasses and introducing buttressing, are basal drag (Vallot and others, Reference Vallot2018), lateral drag (Dupont and Alley, Reference Dupont and Alley2005), sea ice or ice mélange force (Amundson and others, Reference Amundson, Fahnestock, Truffer, Brown, Lüthi and Motyka2010; Meng and others, Reference Meng2025), the ice ligament force beneath the surface crevasse (Buck, Reference Buck2023; Coffey and others, Reference Coffey, Lai, Wang, Buck, Surawy-Stepney and Hogg2024) and the hydrostatic force from the ocean imposed at the calving front.
It will be shown that the outcome of the crack depths and calving criteria will depend on the net buttressing
\begin{equation}
B\left ( x, t \right ) = B_{M{\acute{e}}lange}\left(t\right)
+ B_{Drags}\left ( x \right ).
\end{equation} 
$B_{M\acute{e}lange}$ and BDrags result from mélange and lateral or basal drag forces, respectively, which contribute positively to buttressing. In practice, it could be challenging to quantify each term in (9) due to the poorly known functional forms of these buttressing forces (see section 2.4). However, buttressing can be calculated if the glacial stress states are known. Ice shelf buttressing can be quantified using a similar line integral to (3) (Sergienko, Reference Sergienko2025). For our modeling convenience we define the local dimensionless buttressing B(x) (see Appendix D for more details) as
\begin{equation}
B \equiv 1 - \frac{\overline{R}_{xx}^0}{\overline{R}_{xx}^{B=0}},
\end{equation} where 
$\overline{R}_{xx}^0$ is the depth- and width-averaged resistive stress (defined in (4)) in the unfractured state (Buck and Lai, Reference Buck and Lai2021) and 
$\overline{R}_{xx}^{B=0}$ is the depth- and width-averaged resistive force of a glacier with no buttressing (
$B_{M\acute{e}lange}=B_{Drags}=0$), defined as
\begin{equation}
\overline{R}_{xx}^{B=0} \equiv \frac{1}{2} \left ( 1 - \frac{\rho_i}{\rho_w}\lambda^2 \right ) \rho_i g H,
\end{equation} where 
$\lambda \equiv -\frac{\rho_w}{\rho_i} \frac{b}{H}$ is the dimensionless water level relative to flotation 
$\rho_i H / \rho_w$, and b is the bed elevation relative to the sea level. For ISs λ = 1, and the buttressing definition in (10) converges to that of Gudmundsson Reference Gudmundsson(2013) at the grounding line. Each symbol in (11) is defined in Table 1 and Fig. 2. Thus, the horizontal force per unit width at 
$x=x_c$ in the unfractured state can be written in terms of B as
\begin{equation}
\int_{b}^s \sigma_{xx}^0\left(x=x_c\right) dz \equiv H \overline{R}_{xx}^0 - \frac{1}{2} \rho_i g H^2 = \left ( 1 - B \right ) H \overline{R}_{xx}^{B=0} - \frac{1}{2} \rho_i g H^2.
\end{equation}Table 1. Mathematical Symbols Glossary
We will find after solving for crack depths that there is a range of B where cracks can exist, from crack formation B F to calving 
$B^{*}$, such that crack-depth solutions are valid within these bounds 
$B^{*} \leq B \leq B^{F}$.
To solve for only the surface crack depth, the force balance (8) is sufficient, since it yields an explicit relationship between crack depth, buttressing B and dimensionless water level λ, e.g. (57). To solve for dual crack depths (surface cracks atop basal cracks), both the force balance (8) and crevasse-depth relation (13), i.e. the explicit dependence of the basal crack depth db on the surface crack depth ds, are needed; otherwise, the system is underdetermined. In this paper, we will solve for both surface and dual cracks, as we generalize the previous HFB model (Buck, Reference Buck2023; Coffey and others, Reference Coffey, Lai, Wang, Buck, Surawy-Stepney and Hogg2024; Slater and Wagner, Reference Slater and Wagner2025) to six “crack configurations” illustrated in Fig. 3, each with their relevant setting, i.e. IS, MTG and LTG. As discussed in Buck (Reference Buck2023), Coffey and others (Reference Coffey, Lai, Wang, Buck, Surawy-Stepney and Hogg2024), two conditions are used to determine the surface and basal crack-depth relation: first, stress is continuous at crack tips (Buck and Lai, Reference Buck and Lai2021); second, the ice has zero material strength. In the case of a surface crevasse and a basal crevasse with constant ice and saltwater densities ρi and ρw, this yields a crack depth relation of the form
\begin{equation}
\tilde{d}_b = f\left(\tilde{d}_s, \rho_i, \rho_w, ... \right )
\end{equation}
Figure 3. The six tensile crack configurations considered in this paper. The boxes in the top left corners of each case denote which of the three scenarios—IS, MTG, or LTG—is being considered. Throughout this paper, saltwater and meltwater are shown in blue and light green, respectively. The parameters used to generate these crack depths are B = 0.1, λ = 0.75, 
$\tilde{h}_w = 0.1$, 
$\tilde{z}_h = 0.7$.
where the dimensionless basal and surface crack depths 
$\tilde{d}_b, \tilde{d}_s$ are defined as the basal and surface crevasse depths normalized by ice thickness, i.e. 
$\tilde{d}_{b,s} \equiv d_{b,s}/H$. An example of (13) is (17). Solving (8) and (13) analytically yields simple expressions for tensile crevasse-induced calving laws that can be used to predict calving in numerical ice sheet models. In summary, our recipe for solving for crack depths in this paper is to solve
(1) (8) for HFB, and
(2) (13), if there are dual crevasses, which is obtained in this paper with the assumptions of zero material strength (Nye, Reference Nye1955; Jezek, Reference Jezek1984; Benn and others, Reference Benn, Hulton and Mottram2007a; Nick and others, Reference Nick, van der Veen, Vieli and Benn2010; Buck, Reference Buck2023; Coffey and others, Reference Coffey, Lai, Wang, Buck, Surawy-Stepney and Hogg2024) and continuity of stress at crack tips (Buck and Lai, Reference Buck and Lai2021). Finite ice strength can be added, but it is not within the scope of this paper.
In the Zero-Stress approximation, the HFB (8) is not satisfied (Buck, Reference Buck2023; Coffey and others, Reference Coffey, Lai, Wang, Buck, Surawy-Stepney and Hogg2024). (13) alone is used to solve for crack depths with a background stress state defined in the absence of fractures (Nye, Reference Nye1955; Benn and others, Reference Benn, Hulton and Mottram2007a). Figure 4 demonstrates the issue through an example with a prescribed buttressing 
$B = 1 - x/L$: the Zero-Stress approximation underpredicts the crevasse depths, resulting in the calving stress threshold under the Zero-Stress approximation being twice as large as that predicted by the HFB framework (8) (Buck, Reference Buck2023; Coffey and others, Reference Coffey, Lai, Wang, Buck, Surawy-Stepney and Hogg2024), for a dry surface crevasse and a saltwater basal crevasse. In this paper, we derive the HFB-crack depths of six plausible crack configurations in Fig. 3 for ISs through MTGs to LTGs.
Figure 4. Comparing the crack depths predicted using the HFB and the Zero-Stress approximation for dry surface crevasses and saltwater basal crevasses (DS+SB) given an idealized dimensionless buttressing number of 
$B=1-x/L$. In HFB, crack depths are deeper than that predicted from the Zero-Stress approximation. Importantly, all HFB cases have calving occur at the ice front where B = 0, while the Zero-Stress approximation does not predict calving. According to HFB, instead of a critical stress criteria, zero buttressing B = 0 is the common calving criteria among the ice-shelf, marine-terminating and LTG cases (for a dry surface crevasse and potentially saltwater-filled basal crevasse in the absence of basal melting and material strength). The crack-depth envelopes are plotted as smooth curves, while the jaggedness is plotted to convey that these envelopes represent crack tip depth. Crack spacing is arbitrary in these plots.
2.1. Ice shelf
2.1.1. Dry surface crevasse atop a saltwater basal crevasse (DS+SB)
The simplest dual crack solution exists for a freely floating ice shelf, and was originally derived in Buck Reference Buck(2023). Here, we restate the derivation for completeness, and see in Appendix D that the solution with variable thickness and isostasy is the same. We begin by determining the stress state of the unbroken ice ligament, l, defined in Fig. 5. In an unfractured state, the stress in the unbroken ligament is 
$\sigma_{xx}=-\rho_i g (s-z) + R_{xx}$. In a fractured state, the longitudinal stress at the crevassed location 
$x=x_c$ satisfying the conditions of continuity of stress at the crack tips and zero material strength can be written in a piecewise expression (see derivation in Buck and Lai (Reference Buck and Lai2021); Buck (Reference Buck2023); Coffey and others (Reference Coffey, Lai, Wang, Buck, Surawy-Stepney and Hogg2024))
\begin{align}
&\sigma_{xx} \left (x_c, z \right ) = \nonumber \\
& \left\{
\begin{array}{lr}
0, & s-d_s \leq z \leq s \quad \text{(surface crevasse)} \\
- \rho_i g \left ( s - z \right ) + c R_{xx}, & b+d_b \leq z \leq s-d_s \quad \text{(unbroken ice ligament)} \quad \ \\
\rho_w g z, & b \leq z \leq b+d_b \quad \text{(basal crevasse)} \ \ \,
\end{array}
\right.
\end{align}
Figure 5. Equivalent version of Figure 2 for an ice shelf (IS) with meltwater in a surface crevasse and saltwater in a basal crevasse.
Here, 
$cR_{xx}(x_c,z)$ parameterizes the sum of the crevasse-induced compressive stress in the unbroken ice ligament and 
$R_{xx} = 2 \tau_{xx} + \tau_{yy}$ is the background resistive stress in the unfractured state (see Table 1 for definitions of resistive and deviatoric stresses). The constant c allows the stress in the unbroken ice ligament to update between the fractured and unfractured states, and will be determined shortly. Note that the Zero-Stress approximation corresponds to c = 1, i.e. not allowing the unbroken ice ligament stress state to differ between the unfractured and fractured ice states. For further discussion of 
$c R_{xx}$ and the inconsistency of the Zero-Stress approximation, see section 2.3 and Appendix E of Coffey and others Reference Coffey, Lai, Wang, Buck, Surawy-Stepney and Hogg(2024). The continuity of stress between the crevasse and the unbroken ligament in (14) gives
\begin{equation}
-\rho_i g \left ( s - \left ( s - d_s \right ) \right ) + c R_{xx} = 0 \quad \textrm{at} \quad z\left ( x_c \right ) = s\left ( x_c \right ) - d_s
\end{equation}at the surface crack tip and
\begin{equation}
-\rho_i g \left ( s -\left ( b + d_b \right ) \right ) + c R_{xx} = \rho_w g \left ( b + d_b \right )\ \textrm{at}\ z\left ( x_c \right ) = b\left ( x_c \right ) + d_b
\end{equation}at the basal crack tip. The pair of (15) and (16) removes the unknown c and gives the following crack-depth relation (13), written in dimensionless form as
\begin{align}
\tilde{d}_b = \tilde{d}_s \frac{\rho_i}{\rho_w - \rho_i},
\end{align} where 
$\tilde{d} \equiv d/H$ is the dimensionless crevasse depth and H is the ice thickness at 
$x=x_c$ in the unfractured state (Fig. 5). Note that this crevasse-depth relation (17) assumes isothermal ice. A vertically varying temperature profile would modify the crevasse relation as shown in Coffey and others Reference Coffey, Lai, Wang, Buck, Surawy-Stepney and Hogg(2024).
For an ice shelf, we write the HFB of (8) in its dimensionless form, normalized by 
$H \overline{R}_{xx}^{0}$, as
\begin{equation}
\frac{- \left (1 - \tilde{d}_s - \tilde{d}_b \right )^2 + \tilde{d}_b \frac{\rho_w}{\rho_i} \left (\tilde{d}_b - 2 \frac{\rho_i}{\rho_w} \right )}{1 - \frac{\rho_i}{\rho_w}} = -B - \frac{\rho_i/\rho_w}{1 - \frac{\rho_i}{\rho_w}}.
\end{equation} As in (8), the forces on the left-hand side and right-hand side represent the horizontal forces at 
$x=x_c$ in the fractured and unfractured ice states, respectively. Analytically solving (17) with (18), the crack-depth solutions for an ice shelf that is floating under Archimedean buoyancy everywhere are
\begin{equation}
\textrm{Surface crack depth:} \quad \tilde{d}_s = \left (1 - \frac{\rho_i}{\rho_w} \right ) \left (1 - \sqrt{B} \right ),
\end{equation}
\begin{equation}
\textrm{Basal crack depth:} \quad \tilde{d}_b = \frac{\rho_i}{\rho_w} \left (1 - \sqrt{B} \right ).
\end{equation} We plot these crack depths in Fig. 6a,b. The total fraction of ice that is fractured, 
$\tilde{D}$, can be written in terms of B,
\begin{equation}
\tilde{D} \equiv \tilde{d}_s + \tilde{d}_b = 1 - \sqrt{B}.
\end{equation}
Figure 6. Ice tongues do not form with HFB unless there is a non-zero material strength, positive mass balance, or non-zero buttressing. We demonstrate the case of buttressing with the HFB solutions of (19) and (20) with zero buttressing in panel a and small buttressing in panel b. The ice thickness profile is the analytical solution of Van der Veen Reference van der Veen1986. Crack spacing is arbitrary in these plots. Panel c shows the EPSG:3031 projection of the Drygalski Ice Tongue, Scott Coast, East Antarctica from Sentinel-2 on 7 March 2020 with Highlight Optimized Natural Color from the European Space Agency. Long, bright and dark shadow surface features perpendicular to flow may represent surface depressions atop basal crevasses (Luckman and others, Reference Luckman, Jansen, Kulessa, King, Sammonds and Benn2012).
These algebraic equations provide several insights. First, calving of a varying thickness ice shelf occurs (
$\tilde{D}=1$) if there is no buttressing, setting the lower bound on buttressing 
$B^{*}$ as
\begin{equation}
B^* = 0,
\end{equation}which is the same calving criterion as that for constant thickness ISs as reported in (Buck, Reference Buck2023; Coffey and others, Reference Coffey, Lai, Wang, Buck, Surawy-Stepney and Hogg2024).
Second, the upper bound of buttressing B F for there to be no fractures (fractures of zero depth) and thus no tension at the ice surface is
\begin{equation}
B^{F} = 1.
\end{equation} Thus, for isostatic ISs with varying thicknesses, the nondimensional buttressing (9) must be 
$0 \leq B \leq 1$ to permit the formation of dry surface crevasses and saltwater-filled basal crevasses. The calving criterion is B = 0, as shown by the crack depths in Figs. 6a and 4. Throughout this paper, the lower bound of buttressing 
$B^{*}$ occurs when crack(s) penetrate the entire ice thickness (
$\tilde{D}=1$), and thus is our calving criteria. In general, there will be bounds on the amount of buttressing from full thickness fracture to no fracture, 
$B^{*} \leq B \leq B^{F}$, summarized along with crack-depth formulas in Tables 2, 3 and 4.
Table 2. Crack depths 
$(\tilde{d}_s, \tilde{d}_b)$, calving criteria 
$B^{*}$, buttressing required for dual crack formation B F and the corresponding range of 
$\tilde{h}_w$ for an ice shelf, derived in section 2.1 and illustrated in Figure 3 middle panels. The MS+SB column converges to DS+SB when 
$\tilde{h}_w=0$
Table 3. Crack depths 
$(\tilde{d}_s, \tilde{d}_b)$, calving criteria 
$B^{*}$, buttressing required for surface crack formation B F and the corresponding range of 
$\tilde{h}_w$ for a surface crevasse on a marine-terminating glacier (MTG), derived in sections 2.3.1 and 2.3.2 and illustrated in Figure 3 left panels. The results converge to an LTG when λ = 0. The MS column converges to DS when 
$\tilde{h}_w=0$.
Table 4. Crack depths 
$(\tilde{d}_s, \tilde{d}_b)$, calving criteria 
$B^{*}$, buttressing required for dual crack formation B F and the corresponding range of 
$\tilde{h}_w$ for dual cracks on a marine-terminating glacier (MTG), derived in sections 2.3.3 and 2.3.4 and illustrated in Figure 3 middle and right panels. The results converge to DS+MB/SB when 
$\tilde{h}_w=0$. The MS+MB column converges to an LTG when λ = 0. The MS+SB column converges to an IS when ice is at flotation λ = 1, 
$b/H = - \rho_i/\rho_w$.
According to our HFB model, unbuttressed ISs with 
$B\left(x\right)=0$, i.e. ice tongues, are vulnerable to calving as the predicted surface and basal crevasses meet at the sea level everywhere (Fig. 6a). The existence of un-calved ice tongues can be attributable to non-zero material strength (e.g. Wells-Moran and others Reference Wells-Moran, Ranganathan and Minchew(2025)), non-zero buttressing such as sea ice (Christie and others, Reference Christie, Benham, Batchelor, Rack, Montelli and Dowdeswell2022; Gomez-Fell and others, Reference Gomez-Fell, Rack, Purdie and Marsh2022), and a positive mass balance (Bassis and Ma, Reference Bassis and Ma2015; Lawrence and others, Reference Lawrence2023). This version of HFB does not consider more complicated 3D effects such as basal melt channels and their associated ice shelf flexure (Indrigo and others, Reference Indrigo, Dow, Greenbaum and Morlighem2021) or suture zones (Khazendar and others, Reference Khazendar, Rignot and Larour2011; Jansen and others, Reference Jansen, Luckman, Kulessa, Holland and King2013; Kulessa and others, Reference Kulessa, Jansen, Luckman, King and Sammonds2014; McGrath and others, Reference McGrath2014).
In the following sections, we apply the same procedure to obtain crevasse depths and calving criteria for various crack configurations.
2.1.2. Meltwater-containing surface crevasse atop a saltwater basal crevasse (MS+SB)
In warmer regions, meltwater on the ice shelf surface can enter into surface crevasses and has been implicated in the breakup of ISs. Here, we extend the HFB framework to consider surface crevasse deepening via hydrofracture (Scambos and others, Reference Scambos2009); the effects of surface loading of a lake on the ice shelf surface as in Banwell and others (Reference Banwell, MacAyeal and Sergienko2013), MacAyeal and Sergienko (Reference MacAyeal and Sergienko2013) are left for future research.
The first difference between a meltwater-containing surface crevasse from the previous dry surface crevasse model is the crack-depth relation of (17), since this must now depend on the meltwater depth hw in the surface crevasse. The stress profile at xc in the fractured state also differs from (14) due to meltwater in the surface crevasse,
\begin{align}
&\sigma_{xx} \left (x_c, z \right ) = \nonumber \\
& \left\{
\begin{array}{lr}
0, & s-d_s+h_w \leq z \leq s \quad \text{(air in surface crevasse)} \\
-\rho_m g \left ( s-d_s+h_w - z \right ), & s - d_s \leq z \leq s-d_s+h_w \quad \text{(meltwater in surface crevasse)} \\
- \rho_i g \left ( s - z \right ) + c R_{xx}, & b+d_b \leq z \leq s-d_s \quad \text{(unbroken ice ligament)} \quad \ \\
\rho_w g z, & b \leq z \leq b+d_b \quad \text{(basal crevasse)} \ \ \,
\end{array}
\right.
\end{align} At the surface crevasse tip 
$z = s - d_s$, we have
\begin{equation}
-\rho_m g \left ( s - d_s + h_w - \left ( s - d_s \right ) \right ) = - \rho_i g \left ( s - \left ( s - d_s \right ) \right ) + c R_{xx},
\end{equation} where ρm is the density of meltwater. Similarly, at the basal crevasse tip 
$z = b + d_b$, we have that
\begin{equation}
\rho_w g \left ( b + d_b \right ) = - \rho_i g \left ( s - \left ( b + d_b \right ) \right ) + cR_{xx}.
\end{equation} Combining (25) and (26) while noting that 
$s = \left ( 1 - \frac{\rho_i}{\rho_w} \right ) H$ and 
$b = - \frac{\rho_i}{\rho_w} H$ for a freely floating ice shelf, we determine the crack-depth relation
\begin{equation}
\tilde{d}_s = \frac{\rho_m}{\rho_i} \tilde{h}_w + \left ( \frac{\rho_w}{\rho_i} - 1 \right ) \tilde{d}_b, \quad \textrm{for} \quad 0\leq \tilde{h}_w \leq \frac{\rho_i}{\rho_m},
\end{equation} where 
$\tilde{h}_w\equiv h_w/H$ and its upper bound that permits a dual crack solution will be explained in (37). Note that when 
$\tilde{h}_w = 0$, (17) and (27) are the same.
Next, the force balance in (18) must change to account for the both the hydrostatic meltwater pressure in the surface crevasse and the subsequent change in the stress in the unbroken ligament, resulting in
\begin{align}
& \frac{- \frac{\rho_m}{\rho_i} \tilde{h}_w^2 - \left (1 - \tilde{d}_s - \tilde{d}_b \right )^2 - 2 \frac{\rho_m}{\rho_i} \tilde{h}_w \left (1 - \tilde{d}_s - \tilde{d}_b \right ) + \tilde{d}_b \frac{\rho_w}{\rho_i} \left (\tilde{d}_b - 2 \frac{\rho_i}{\rho_w} \right )}{1 - \frac{\rho_i}{\rho_w}}\nonumber \\ & \qquad = -B - \frac{\rho_i/\rho_w}{1 - \frac{\rho_i}{\rho_w}}.
\end{align} In the limiting case of 
$\tilde{h}_w = 0$, (18) and (28) would be identical.
Solving (27) and (28), we find that
\begin{align}
&\textrm{Surface crack depth:} \quad \tilde{d}_s = \frac{\rho_m}{\rho_i} \tilde{h}_w + \left (1 - \frac{\rho_i}{\rho_w} \right ) \nonumber \\ & \quad \left (1 - \sqrt{B + \frac{\rho_m - \rho_i}{\rho_w - \rho_i} \frac{\rho_m \rho_w}{\rho_i^2} \tilde{h}_w^2} \right ) \nonumber \\ & \quad \textrm{for} \quad 0\leq \tilde{h}_w \leq \frac{\rho_i}{\rho_m},
\end{align}
\begin{align}
& \textrm{Basal crack depth:} \quad \tilde{d}_b = \frac{\rho_i}{\rho_w} \left (1 - \sqrt{B + \frac{\rho_m - \rho_i}{\rho_w - \rho_i} \frac{\rho_m \rho_w}{\rho_i^2} \tilde{h}_w^2 } \right ) \nonumber\\ & \quad \textrm{for} \quad 0\leq \tilde{h}_w \leq \frac{\rho_i}{\rho_m}.
\end{align} We now construct bounds for buttressing 
$B^{*} \leq B \leq B^{F}$ as before. The minimum stress or maximum buttressing to permit basal crack formation is
\begin{equation}
B^{F} = 1 - \frac{\rho_m - \rho_i}{\rho_w - \rho_i} \frac{\rho_m \rho_w}{\rho_i^2} \tilde{h}_w^2 \quad \textrm{for} \quad 0\leq \tilde{h}_w \leq \frac{\rho_i}{\rho_m}.
\end{equation}
\begin{equation}
\tilde{D} = \frac{\rho_m}{\rho_i} \tilde{h}_w + 1 - \sqrt{B + \frac{\rho_m - \rho_i}{\rho_w - \rho_i} \frac{\rho_m \rho_w}{\rho_i^2}\tilde{h}_w^2},
\end{equation} thus calving or 
$\tilde{D} = 1$ occurs when
\begin{equation}
B^* = \frac{\rho_w - \rho_m}{\rho_w - \rho_i} \frac{\rho_m}{\rho_i} \tilde{h}_w^2 \quad \textrm{for} \quad 0\leq \tilde{h}_w \leq \frac{\rho_i}{\rho_m}.
\end{equation} Comparing (33) with (22), we see that adding meltwater allows calving to occur in buttressed regions of the ice shelf. Comparing (31) with (33), we see that B F decreases with 
$\tilde{h}_w^2$, whereas 
$B^{*}$ increases with 
$\tilde{h}_w^2$. This buttressing range 
$B^{*} \leq B \leq B^{F}$ permits the dual crack configuration to still obey the force balance constraint. Continuously adding surface meltwater 
$\tilde{h}_w$ will either cause calving or provide enough pressure to close the basal crevasse. In the latter case, 
$B^F \leq B^*$, and a transition to a new crack configuration, i.e. a meltwater-containing surface crevasse without a basal crevasse (MS), must occur. By adding meltwater to highly buttressed ice shelf regions, calving is reached by a meltwater-containing surface crevasse alone.
2.1.3. Meltwater-containing surface crevasse without a basal crevasse (MS)
As in the previous section, the stress balance at the surface crevasse tip (25) is identical, but there is no basal crevasse and thus no need for (26). The force-balance equation
\begin{equation}
\frac{- \frac{\rho_m}{\rho_i} \tilde{h}_w^2 - \left ( 1 - \tilde{d}_s \right )^2 - 2 \frac{\rho_m}{\rho_i} \tilde{h}_w \left ( 1 - \tilde{d}_s \right )}{1-\frac{\rho_i}{\rho_w}} = -B - \frac{\rho_i/\rho_w}{1 - \frac{\rho_i}{\rho_w}}
\end{equation} differs from (28) in that the basal crevasse does not exist 
$\tilde{d}_b = 0$. The analytical solution to (34) is
\begin{align}
&\textrm{Surface crack depth:} \quad \tilde{d}_s = 1 + \frac{\rho_m}{\rho_i} \tilde{h}_w \nonumber\\ & \quad - \sqrt{B \left(1-\frac{\rho_i}{\rho_w} \right ) + \frac{\rho_i}{\rho_w} + \frac{\rho_m}{\rho_i} \left ( \frac{\rho_m}{\rho_i} - 1 \right ) \tilde{h}_w^2} \quad \textrm{for} \nonumber \\ & \quad \frac{\rho_i}{\rho_m} \leq \tilde{h}_w \leq 1.
\end{align} Calving occurs when 
$\tilde{d}_s = 1$, and we may substitute this into (35) to solve for the calving threshold,
\begin{equation}
B^* = \frac{- \frac{\rho_i}{\rho_w} + \frac{\rho_m}{\rho_i} \tilde{h}_w^2 }{1 - \frac{\rho_i}{\rho_w}} \quad \textrm{for} \quad \frac{\rho_i}{\rho_m} \leq \tilde{h}_w \leq 1.
\end{equation} Calving is possible at higher levels of buttressing due to larger amounts of meltwater 
$\tilde{h}_w$. To determine the transition from dual crevasses to a meltwater-containing surface crevasse without a basal crevasse (MS), we must set (31) and (33) equal to solve for meltwater depth 
$\tilde{h}_w^T$ at the limit of 
$B^* = B^F$. This gives the buttressing BT and water depth 
$\tilde{h}_w^T$ at the transition from calving by dual crevasses to calving by a meltwater-containing surface crevasse.
\begin{equation}
B^{T} = \frac{\frac{\rho_w}{\rho_m} - 1}{\frac{\rho_w}{\rho_i} - 1}, \quad \textrm{and} \quad \tilde{h}_w^{T} = \frac{\rho_i}{\rho_m}.
\end{equation} Thus, in HFB, the meltwater-containing surface crevasse without a basal crevasse (MS) case on an ice shelf exists for 
$B^{T} \leq B^{*} \leq B$, defined in (37) and (36). For calving to occur, the lowest amount of meltwater in this configuration is 
$\tilde{h}_w^{T} = \rho_i/\rho_m\approx0.917$, showing that more than 
$90\%$ of the ice thickness must be filled with meltwater in a surface crevasse to cause calving from a meltwater-containing surface crevasse without a basal crevasse (MS) on an ice shelf. This is consistent with the previously reported results using Linear Elastic Fracture Mechanics (LEFM) (Lai and others, Reference Lai2020).
2.2. Land-terminating glacier
LTGs are defined here in the absence of water at the ice front (Figs. 1 and 7). As such, no saltwater basal crevasses can form. Note that in reality, LTGs often terminate in a toe or snout instead of an idealized vertical cliff as in Fig. 7. Additionally, basal crevasses on LTGs forced by basal meltwater is a hypothetical scenario and may be challenging to observe. However, the derivation of these hypothetical cases in this subsection offers a useful theoretical comparison with MTGs, which we explore in the next subsection. Specifically, assuming a flat bed and constant thickness, the LTG is equivalent to the MTG with λ = 0, and (11) becomes
\begin{equation}
\overline{R}_{xx}^{B=0}\left(\lambda=0\right) = \frac{1}{2} \rho_i g H.
\end{equation}
Figure 7. Equivalent version of Figure 2 for an LTG with meltwater in crevasses. The ocean force and floating ice mélange are absent for an LTG.
Distinct from the previous ice shelf cases, the scenarios considered for LTGs are dry surface crevasses (DS), meltwater-containing surface crevasses (MS) and a combination of either of those surface crevasses with a meltwater basal crevasse (DS+MB/MS+MB). We now elaborate on each case.
2.2.1. Dry surface crevasse without a basal crevasse (DS)
The case of a dry tensile surface crevasse on an LTG represents a simple starting point for considering stability to fracture. Since we do not have a basal crevasse, there is no relation between crack depths (13). Instead of having two equations, (13) and (8), with two unknowns, 
$\tilde{d}_b$ and 
$\tilde{d}_s$ given B, we have one equation (8) with one unknown 
$\tilde{d}_s$ given B.
As before, we use continuity of stress at the surface crevasse tip given the zero-strength assumption to solve for the stress state in the unbroken ice. At the surface crack tip 
$x=x_c, z\left(x_c\right)=s\left(x_c\right)-d_s$, the stress balance is
\begin{equation}
-\rho_i g \left ( s - \left ( s - d_s \right ) \right ) + c R_{xx} = 0.
\end{equation}Solving for the constant c gives the stress profile in the ice
\begin{equation}
\sigma_{xx} \left ( x=x_c, z \leq s-d_s \right ) = - p_l + cR_{xx} = - \rho_i g \left ( s - d_s - z \right ).
\end{equation}Then we may write the force balance (8) in dimensional form as
\begin{equation}
- \frac{\rho_i g}{2} \left ( H - d_s \right )^2 = \left ( 1 - B \right ) H \overline{R}_{xx}^{B=0} - \frac{\rho_i g H^2}{2}.
\end{equation}(41) gives the nondimensional surface crevasse depth,
\begin{equation}
\textrm{Surface crack depth:} \quad \tilde{d}_s = 1 - \sqrt{B }.
\end{equation}Thus, for the surface crack to form,
\begin{equation}
B^F = 1,
\end{equation}and for calving, the buttressing is
\begin{equation}
B^* = 0.
\end{equation} This equation shows that for an LTG with zero material strength, a dry surface crevasse will reach the base of the glacier in the absence of any buttressing, e.g. no basal or lateral drag, 
$B=B^*=0$ (shown in purple in Fig. 4 and in Fig. 9i). For there to be no cracks, i.e. 
$\tilde{d}_s = 0$, the buttressing must be 
$B =B^F= 1$. These 
$B^*, B^F$ results are the same as the dual crack configuration for the IS in (21).
2.2.2. Meltwater-containing surface crevasse without a basal crevasse (MS)
The addition of meltwater in a surface crevasse changes both the stress in the unbroken ice ligament (40) and adds a water pressure term to the force balance in (41). To satisfy continuity of stress at the crack tip, the analogous equation to the dry surface crevasse case (39) is
\begin{equation}
- \rho_i g \left ( s - \left ( s - d_s \right ) \right ) + c R_{xx} = - \rho_m g h_w.
\end{equation}Thus, the stress in the ice below the surface crevasse, analogous to (40), is
\begin{equation}
\sigma_{xx} \left ( x = x_c, z \right ) = - \rho_i g \left ( s - d_s - z \right ) - \rho_m g h_w.
\end{equation}Similarly, the force-balance equation analogous to (41) now has two new terms from the water pressure and updated stress in the unbroken ligament,
\begin{align}
& - \frac{\rho_i g}{2} \left ( H - d_s \right )^2 - \rho_m g h_w \left ( H - d_s \right ) - \frac{\rho_m g h_w^2}{2} \nonumber \\ & \qquad = \left ( 1 - B \right ) H \overline{R}_{xx}^{B=0} - \frac{\rho_i g H^2}{2}.
\end{align} This quadratic equation can be simplified and solved for the dimensionless crack depth 
$\tilde{d}_s$,
\begin{align}
&\textrm{Surface crack depth:} \quad \tilde{d}_s = 1 + \frac{\rho_m}{\rho_i} \tilde{h}_w \\ & \nonumber - \sqrt{B + \frac{\rho_m}{\rho_i} \left (\frac{\rho_m}{\rho_i} -1 \right ) \tilde{h}_w^2} \quad \text{for} \quad 0 \leq \tilde{h}_w \leq 1.
\end{align} In the absence of meltwater or 
$\tilde{h}_w = 0$, this equation converges to the dry surface crevasse case (41). We now discuss the range of 
$B^* \leq B \leq B^F$ that permits this crack-depth solution.
First, the smallest crack depth in (48) is 
$\tilde{d}_s = \tilde{h}_w$. Using this equality, one can solve for the maximum buttressing that permits a meltwater hydrofracture,
\begin{equation}
B^{F} = 1 + \left (\frac{\rho_m}{\rho_i} -1 \right ) \tilde{h}_w\left (2 - \tilde{h}_w \right ) \quad \text{for} \quad 0 \leq \tilde{h}_w \leq 1.
\end{equation} When 
$\tilde{h}_w = 0$ this equation converges to the dry surface crevasse case. However, when there is meltwater 
$\tilde{h}_w \gt 0$, the buttressing upper bound is larger, indicating that hydrofracture is less stable than dry fracture.
Second, the case of calving is determined by setting 
$\tilde{d}_s = 1$ in (48),
\begin{equation}
B^{*} = \frac{\rho_m}{\rho_i} \tilde{h}_w^2 \quad \text{for} \quad 0 \leq \tilde{h}_w \leq 1.
\end{equation} Unlike the dry surface crevasse case, calving due to a meltwater surface crevasse can now occur in buttressed regions B > 0. Ice becomes less stable with a larger density ratio 
$\rho_m/\rho_i$ and with more meltwater 
$\tilde{h}_w \equiv h_w / H$.
2.2.3. Surface crevasse atop a meltwater basal crevasse (DS+MB/MS+MB)
As seen in the previous section, the solution for the dry surface crevasse on an LTG is a limiting case of the meltwater-containing surface crevasse when there is no meltwater, or 
$\tilde{h}_w = 0$. As such, we now derive the dual solution for a surface crevasse, dry or with meltwater and a basal crevasse filled with subglacial meltwater.
The stress in the unbroken ice ligament has the same form as (46). However, since we have two crevasses, we now seek a crack-depth relation as outlined in (13). Continuity of stress, as applied in (39) and (45) at the surface crevasse tip, is applied at the basal crack tip and gives the crack-depth relation
\begin{equation}
\tilde{d}_s = 1 - \tilde{d}_b - \frac{\rho_m}{\rho_i} \left (\tilde{z}_h - \tilde{h}_w - \tilde{d}_b \right ),
\end{equation} where 
$\tilde{z}_h\equiv z_h/H$ and zh is the piezometric head height, i.e. the height to which the water would rise relative to the ice bed 
$z=b \lt 0$ in a borehole. Note that in dimensional form, the hydrostatic water pressure in the basal crevasse is set to be 
$- \rho_m g \left ( b + z_h - z \right )$.
Next, we specify the HFB of (8). Relative to (47), there is an additional force due to the water pressure in the basal crevasse, and an adjustment to the ice pressure forces at 
$x=x_c$ because the ligament l is now smaller with dual cracks (see Fig. 7). Written in a nondimensional form, the HFB equation becomes
\begin{align}
& - \left ( 1 - \tilde{d}_s - \tilde{d}_b \right )^2 - 2 \frac{\rho_m}{\rho_i} \tilde{h}_w \left ( 1 - \tilde{d}_s - \tilde{d}_b \right ) - \frac{\rho_m}{\rho_i} \tilde{h}_w^2 \nonumber \\ & \qquad - \frac{\rho_m}{\rho_i} \tilde{d}_b \left ( 2 \tilde{z}_h - \tilde{d}_b \right ) = -B.
\end{align}With the crack-depth relation (51) and HFB (52), we analytically solve for crack-depth solutions,
\begin{align}
&\textrm{Surface crack depth:} \quad \tilde{d}_s = 1 + \frac{\rho_m}{\rho_i}\tilde{h}_w - \tilde{z}_h \nonumber \\ & - \sqrt{\left (1 - \frac{\rho_i}{\rho_m} \right ) \left[ B - \frac{\rho_m}{\rho_i} \left( \tilde{z}_h^2 - \left (\frac{\rho_m-\rho_i}{\rho_i}\right) \tilde{h}_w^2 \right)\right]}\\
&\quad \text{for} \quad 0 \leq \tilde{h}_w \leq \tilde{z}_h, \nonumber
\end{align}
\begin{align}
&\textrm{Basal crack depth:} \quad \tilde{d}_b = \tilde{z}_h \nonumber \\ & - \frac{\rho_i}{\rho_m}\sqrt{\frac{\rho_m}{\rho_m - \rho_i} \left[ B - \frac{\rho_m}{\rho_i} \left( \tilde{z}_h^2 - \left (\frac{\rho_m-\rho_i}{\rho_i}\right) \tilde{h}_w^2 \right)\right] } \quad \text{for} \nonumber \\ & \quad 0 \leq \tilde{h}_w \leq \tilde{z}_h.
\end{align} We next seek the buttressing range that permits crack formation, 
$B^{*} \leq B \leq B^{F}$. Unlike the case of MS without a basal crevasse (49), the smallest crack depth in this dual crack solution occurs when there is no basal crevasse, 
$\tilde{d}_b = 0$. Plugging this in to (54) gives
\begin{equation}
B^F = \frac{\rho_m}{\rho_i} \left ( \frac{\rho_m}{\rho_i} \tilde{z}_h^2 - \frac{\rho_m - \rho_i}{\rho_i} \tilde{h}_w^2 \right ) \quad \text{for} \quad 0 \leq \tilde{h}_w \leq \tilde{z}_h.
\end{equation} The calving criterion can be determined from (53) and (54) with 
$\tilde{D} \equiv \tilde{d}_s+\tilde{d}_b= 1$,
\begin{equation}
B^* = \frac{\rho_m}{\rho_i} \tilde{z}_h^2 \quad \text{for} \quad 0 \leq \tilde{h}_w \leq \tilde{z}_h.
\end{equation} Interestingly, the calving threshold for a meltwater-containing surface crevasse over a meltwater basal crevasse does not depend on the amount of water in the surface crevasse. However, the buttressing bounds 
$B^{*} \leq B \leq B^{F}$ require that 
$B^* \leq B^F$: according to (56) and (55), this inequality holds on LTGs when 
$\tilde{z}_h \geqslant \tilde{h}_w$. Thus, calving is determined by the basal crevasse for the case of a meltwater-containing surface crevasse atop a meltwater basal crevasse, and this dual crack configuration is only valid when the head height of the meltwater is at least larger than the meltwater depth in the surface crevasse, 
$\tilde{z}_h \geqslant \tilde{h}_w$. However, if the meltwater basal crevasse closes because 
$\tilde{h}_w \geqslant \tilde{z}_h$, the calving threshold would be set by the meltwater-crevasse without a basal crevasse (MS) case defined in the previous section (50).
2.3. Marine-terminating glacier
MTGs cover the widest range of scenarios and have crack solutions converge to LTGs when the water at the ice front is zero. MTGs also have crack solutions converge to ISs when the basal crack is saltwater-filled and the ice is at flotation everywhere. Thus, this section contains the most generalized crack solutions that, under some limits, converge to the previously presented cases. Specifically, to maintain a consistent definition of buttressing between ISs, MTGs and LTGs, the MTG case considers constant thickness and flat bed slope (see Appendix D).
2.3.1. Dry surface crevasse without a basal crevasse (DS)
The force-balance equation is the similar to that of LTGs (41), but is now dependent on the dimensionless water level 
$\lambda\equiv \frac{\rho_w}{\rho_i} \frac{-b}{H}\geqslant 0$, which will be important throughout the MTG section,
\begin{equation}
- \frac{\rho_i g}{2} \left ( H - d_s \right )^2 = \left ( 1 - B \right ) H \overline{R}_{xx}^{B=0} - \frac{\rho_i g H^2}{2}.
\end{equation}As before, we use ρw as the saltwater density, but it can be defined for alternative applications as the lake water density. (57) gives the surface crevasse depth,
\begin{equation}\textrm{Surface crack depth:} \quad \tilde{d}_s = 1 - \sqrt{ B \left ( 1 - \frac{\rho_i}{\rho_w}\lambda^2 \right ) + \frac{\rho_i}{\rho_w}\lambda^2 }.\end{equation}This solution is shown by the red curve in Fig. A1 in Appendix A. A proglacial body of water alone acts to reduce the depth of cracks. For there to be no crevasse, we would have that
\begin{equation}B^F=1,\end{equation}while calving would occur when
\begin{equation}B^\ast=\frac{-\frac{\rho_i}{\rho_w}\lambda^2}{1-\frac{\rho_i}{\rho_w}\lambda^2}.\end{equation}Thus, negative buttressing is needed for calving to occur via a dry surface crevasse in a MTG. According to (9), negative buttressing would not occur in our width-averaged framework.
2.3.2. Meltwater-containing surface crevasse without a basal crevasse (MS)
Following the corresponding section for an LTG, we extend the derivation to include proglacial water at the calving front. The stress in the ice ligament from the base to the surface crevasse tip has the same form as (46). The force balance (47) now has the slightly modified form,
\begin{align}
&- \frac{\rho_i g}{2} \left ( H - d_s \right )^2 - \rho_m g h_w \left ( H - d_s \right ) - \frac{\rho_m g h_w^2}{2} \nonumber \\& = \left ( 1 - B \right ) H \overline{R}_{xx}^{B=0} - \frac{\rho_i g H^2}{2}.
\end{align}The solution for surface crevasse depth then becomes
\begin{align}
&\textrm{Surface crack depth:} \quad \tilde{d}_s = 1 + \frac{\rho_m}{\rho_i} \tilde{h}_w \nonumber \\ & - \sqrt{B \left ( 1 - \frac{\rho_i}{\rho_w}\lambda^2 \right ) + \frac{\rho_i}{\rho_w}\lambda^2 + \frac{\rho_m}{\rho_i} \left (\frac{\rho_m}{\rho_i} -1 \right ) \tilde{h}_w^2} \quad \text{for} \nonumber \\ & \quad 0 \leq \tilde{h}_w \leq 1.
\end{align} To have the minimal meltwater-filled crack depth 
$\tilde{d}_s=\tilde{h}_w$, the buttressing must be
\begin{equation}
B^{F} = 1 + \frac{\left (\frac{\rho_m}{\rho_i} -1 \right ) \tilde{h}_w\left (2 - \tilde{h}_w \right ) }{1 - \frac{\rho_i}{\rho_w}\lambda^2} \quad \text{for} \quad 0 \leq \tilde{h}_w \leq 1.
\end{equation}To have calving, the buttressing must be
\begin{equation}
B^{*} = \frac{- \frac{\rho_i}{\rho_w}\lambda^2 + \frac{\rho_m}{\rho_i} \tilde{h}_w^2}{1 - \frac{\rho_i}{\rho_w}\lambda^2} \quad \text{for} \quad 0 \leq \tilde{h}_w \leq 1.
\end{equation} It is clear that the crack depths, 
$B^*$, and BF all converge to the results in the LTG cases when λ = 0.
2.3.3. Surface crevasse atop a meltwater basal crevasse (DS+MB/MS+MB)
We now extend the previous case to consider a basal crevasse filled with meltwater, e.g. from the subglacial water, underneath a surface crevasse. We follow a similar format to that for the LTG in section 2.2.3. The stress in the ice again has the form of (46), and the crack-depth relation is defined by (51). For the force balance equation, the only change is the water pressure at the ice front in the force balance of (52),
\begin{align}
&\frac{- \left ( 1 - \tilde{d}_s - \tilde{d}_b \right )^2 - 2 \frac{\rho_m}{\rho_i} \tilde{h}_w \left ( 1 - \tilde{d}_s - \tilde{d}_b \right ) - \frac{\rho_m}{\rho_i} \tilde{h}_w^2 - \frac{\rho_m}{\rho_i} \tilde{d}_b \left ( 2 \tilde{z}_h - \tilde{d}_b \right )}{1 - \frac{\rho_i}{\rho_w}\lambda^2} \nonumber \\ & = -B - \frac{\frac{\rho_i}{\rho_w}\lambda^2}{1 - \frac{\rho_i}{\rho_w}\lambda^2}.
\end{align}Solving (65) with (51) for dimensionless crevasse depths, we find that
\begin{align}
&\textrm{Surface crack depth:} \quad \tilde{d}_s = 1 + \frac{\rho_m}{\rho_i}\tilde{h}_w - \tilde{z}_h\\
&- \sqrt{\left (1 - \frac{\rho_i}{\rho_m} \right ) \left[ B \left ( 1 - \frac{\rho_i}{\rho_w}\lambda^2 \right ) + \frac{\rho_i}{\rho_w}\lambda^2 - \frac{\rho_m}{\rho_i} \left( \tilde{z}_h^2 - \left (\frac{\rho_m-\rho_i}{\rho_i}\right) \tilde{h}_w^2 \right)\right]} \nonumber \\
&\quad \text{for} \quad 0 \leq \tilde{h}_w \leq \tilde{z}_h, \nonumber
\end{align}
\begin{align}
&\textrm{Basal crack depth:} \quad \tilde{d}_b = \tilde{z}_h - \frac{\rho_i}{\rho_m} \nonumber \\ &\sqrt{\frac{\rho_m}{\rho_m - \rho_i} \left[ B \left ( 1 - \frac{\rho_i}{\rho_w}\lambda^2 \right ) + \frac{\rho_i}{\rho_w}\lambda^2 - \frac{\rho_m}{\rho_i} \left( \tilde{z}_h^2 - \left (\frac{\rho_m-\rho_i}{\rho_i}\right) \tilde{h}_w^2 \right)\right] }\\
&\quad \text{for} \quad 0 \leq \tilde{h}_w \leq \tilde{z}_h. \nonumber
\end{align}When the basal crevasse depth is zero, the corresponding buttressing is
\begin{equation}
B^F = \frac{- \frac{\rho_i}{\rho_w}\lambda^2 + \frac{\rho_m}{\rho_i} \left ( \frac{\rho_m}{\rho_i} \tilde{z}_h^2 - \frac{\rho_m - \rho_i}{\rho_i} \tilde{h}_w^2 \right ) }{1 - \frac{\rho_i}{\rho_w}\lambda^2} \quad \text{for} \quad 0 \leq \tilde{h}_w \leq \tilde{z}_h.
\end{equation} Calving, i.e. 
$\tilde{D}\equiv \tilde{d}_s+\tilde{d}_b=1$, occurs when the buttressing number reaches
\begin{equation}
B^* = \frac{- \frac{\rho_i}{\rho_w}\lambda^2 + \frac{\rho_m}{\rho_i} \tilde{z}_h^2 }{1 - \frac{\rho_i}{\rho_w}\lambda^2} \quad \text{for} \quad 0 \leq \tilde{h}_w \leq \tilde{z}_h.
\end{equation} The crack depths, 
$B^*$, and BF converge to the results in the LTG cases when λ = 0. Note the lack of dependence of the calving stress threshold on the amount of meltwater in the surface crevasse, 
$\tilde{h}_w$. This is the same situation as the LTG; as in section 2.2.3, for 
$B^* \leq B \leq B^F$ to hold, this dual crevasse case can only exist when 
$\tilde{z}_h \geqslant \tilde{h}_w$. Thus, a stable meltwater basal crevasse will not form beneath a meltwater-containing surface crevasse unless 
$\tilde{z}_h \geqslant \tilde{h}_w$. However, if the meltwater basal crevasse closes because 
$\tilde{h}_w \geqslant \tilde{z}_h$, the calving threshold would be set by the MS case without a basal crevasse defined in the previous section (64).
2.3.4. Surface crevasse atop a saltwater basal crevasse (DS+SB/MS+SB)
In this section, we will study a simplified version of a saltwater-filled basal crevasse beneath a surface crevasse. Saltwater intrusions have been studied theoretically (Wilson and others, Reference Wilson, Wells, Hewitt and Cenedese2020; Robel and others, Reference Robel, Wilson and Seroussi2022; Gadi and others, Reference Gadi, Rignot and Menemenlis2023; Bradley and Hewitt, Reference Bradley and Hewitt2024; Ehrenfeucht and others, Reference Ehrenfeucht, Rignot and Morlighem2024) and observationally (e.g. Kim and others, Reference Kim, Rignot and Holland2024; Rignot and others, Reference Rignot, Ciracì, Scheuchl, Tolpekin, Wollersheim and Dow2024). As the precise form of water pressure and density will not have a simple analytical form and likely evolve with entrainment, tides and grounding line migration, we develop an end-member model for fully saltwater-filled basal crevasses.
Similar to the case of a meltwater-containing surface crevasse on an ice shelf, which combines (25) and (26), we have the crack-depth relation
\begin{equation}
\tilde{d}_s = 1+\frac{\rho_m}{\rho_i} \tilde{h}_w + \frac{\rho_w-\rho_i}{\rho_i} \tilde{d}_b -\lambda \quad \text{for} \quad 0 \leq \tilde{h}_w \leq -\frac{\rho_w}{\rho_m} \frac{b}{H}.
\end{equation}The dimensionless force-balance equation is
\begin{align}
&\left [ - \frac{\rho_m}{\rho_i} \tilde{h}_w^2 - \left ( 1 - \tilde{d}_s - \tilde{d}_b \right )^2 - 2 \frac{\rho_m}{\rho_i} \tilde{h}_w \left ( 1 - \tilde{d}_s - \tilde{d}_b \right )\right. \nonumber \\ & \quad \left.+ \ \tilde{d}_b \frac{\rho_w}{\rho_i} \left ( 2 \frac{b}{H} + \tilde{d}_b \right ) \right ] / \left [ 1 - \frac{\rho_i}{\rho_w}\lambda^2 \right ] = -B - \frac{\frac{\rho_i}{\rho_w}\lambda^2}{1 - \frac{\rho_i}{\rho_w}\lambda^2}.
\end{align}Combined, we find that
\begin{align}
&\textrm{Surface crack depth:} \quad \tilde{d}_s = 1 + \frac{b}{H} + \frac{\rho_m}{\rho_i} \tilde{h}_w\nonumber \\ & - \sqrt{\left(1 - \frac{\rho_i}{\rho_w}\right) \left[B \left(1 - \frac{\rho_i}{\rho_w}\lambda^2 \right) + \frac{\rho_m}{\rho_i} \frac{\rho_m - \rho_i}{\rho_i} \tilde{h}_w^2 \right ]} \nonumber \\
&\text{for} \quad 0 \leq \tilde{h}_w \leq -\frac{\rho_w}{\rho_m} \frac{b}{H},
\end{align}
\begin{align}
&\textrm{Basal crack depth:} \quad \tilde{d}_b = - \frac{b}{H} \nonumber \\ &- \frac{\rho_i}{\rho_w} \sqrt{\frac{1}{1 - \frac{\rho_i}{\rho_w}} \left[B \left(1 - \frac{\rho_i}{\rho_w}\lambda^2 \right) + \frac{\rho_m}{\rho_i} \frac{\rho_m - \rho_i}{\rho_i} \tilde{h}_w^2 \right ]} \quad \text{for} \nonumber \\ & \quad 0 \leq \tilde{h}_w \leq -\frac{\rho_w}{\rho_m} \frac{b}{H}.
\end{align} The dual crack-depth solutions as a function of buttressing are presented in Fig. 8 without surface meltwater 
$\tilde{h}_w=0$. The range of buttressing that permits dual crack formation is 
$B^* \leq B \leq B^F$, where 
$B=B^F$ when 
$\tilde{d}_b=0$ (denoted by the blue stars in Fig. 8),
\begin{equation}
B^F = \frac{\left ( 1 - \frac{\rho_i}{\rho_w} \right ) \lambda^2 - \frac{\rho_m}{\rho_i} \frac{\rho_m - \rho_i}{\rho_i} \tilde{h}_w^2}{1 - \frac{\rho_i}{\rho_w}\lambda^2} \quad \text{for} \quad 0 \leq \tilde{h}_w \leq -\frac{\rho_w}{\rho_m} \frac{b}{H},
\end{equation}
Figure 8. Crevasse-depth solutions for Marine-Terminating Glaciers (MTGs) (58), (72) and (73) as a function of dimensionless buttressing B for dry surface and saltwater basal crevasses (DS+SB). Colors correspond to the dimensionless water level, 
$\lambda \equiv-\frac{\rho_w}{\rho_i} \frac{b}{H}$. Solid lines are basal crack depths 
$d_b/H$ measured from the ice base at 0. Dash-dotted lines are surface crack depths 
$d_s/H$ measured downward from the surface. For all cases with dry surface and saltwater basal crevasses (DS+SB), the calving criterion (
$d_s+d_b=H$) of MTGs is 
$B^{*} = 0$ (marked with yellow stars) as seen in (75). The unlikely super-buoyant scenario, λ > 1 (Benn and others, Reference Benn2017) is represented here with the red curves. Intruding saltwater under grounded MTGs do not form basal crevasses unless 
$B\leq B^{F}$, defined by (74) and shown by the blue stars.
while if 
$B=B^*$ then calving occurs (
$\tilde{D}\equiv\tilde{d}_s+\tilde{d}_b=1$; denoted by the yellow stars in Fig. 8),
\begin{equation}
B^* = \frac{\frac{\rho_m}{\rho_i} \left ( 1 - \frac{\rho_m}{\rho_w} \right ) \tilde{h}_w^2}{1 - \frac{\rho_i}{\rho_w}\lambda^2} \quad \text{for} \quad 0 \leq \tilde{h}_w \leq -\frac{\rho_w}{\rho_m} \frac{b}{H}.
\end{equation} There are several interesting limits for this calving threshold. First, if there is no water in the surface crevasse or if 
$\rho_w=\rho_m$, then the calving threshold is 
$B^{*}=0$ regardless of the value of λ (denoted by the yellow stars in Fig. 8). Thus, if there is no net source of buttressing, the glacier would calve. Second, more meltwater in the surface crack will increase the magnitude of 
$B^{*}$, making ice more vulnerable to calving. Lastly, in the case of 
$\rho_w=\rho_m$ where the density of the water in the surface and basal crevasse is the same, the calving threshold does not depend on the amount of water in the surface crevasse, hw. Similar to the meltwater basal crevasse case, this holds as long as this dual crack configuration can stably exist, or 
$B^* \leq B \leq B^F$. By evaluating 
$B^* \leq B^F$ with (74) and (75), the stable dual crevasse configuration can exist so long as 
$\tilde{h}_w \leq -\frac{\rho_w}{\rho_m} \frac{b}{H}$. However, if the saltwater basal crevasse closes because 
$\tilde{h}_w \geqslant -\frac{\rho_w}{\rho_m} \frac{b}{H}$, the calving threshold would be set by the MS case without a basal crevasse defined in the previous section 2.3.2.
2.4. Buttressing and the causes of calving
Much of the results of crevasse depths and calving criteria depend on the buttressing number B (10). In Fig. 4, we used an idealized buttressing distribution 
$B\left ( \frac{x}{L} \right )= 1 - \frac{x}{L}$ to illustrate how buttressing affects crack depths. To understand how buttressing can evolve with time for an MTG on a flat bed and constant thickness, as in (9) the buttressing number (Appendix D) may be written as
\begin{align}
B (x,t)= & \phi \left ( \frac{H_M\left(t\right)}{H} \right )^2 \frac{1 - \frac{\rho_i}{\rho_w}}{1 - \frac{\rho_i}{\rho_w} \lambda^2}\nonumber \\ &+ \frac{\int_{x}^{x_t} \tau_{bx}\left(x', t\right) dx' + \int_{x}^{x_t} {\tau}_w\left(x', t\right) dx'}{\frac{1}{2} \left ( 1 - \frac{\rho_i}{\rho_w} \lambda^2 \right ) \rho_i g H^2}.
\end{align} The first term is the positive buttressing provided by the floating ice mélange (Amundson and others, Reference Amundson, Robel, Burton and Nissanka2025; Meng and others, Reference Meng2025) with porosity ϕ and thickness 
$H_M\left(t\right)$. When mélange thickness is very small, as most times in the summer, the mélange buttressing is negligible. The second term comes from the lateral and basal drag forces exerted on ice, consisting of the force per unit width (into the page) 
$\int_{x}^{x_t} \tau_{bx} dx'$ due to the basal drag along the bedrock and the horizontal force per unit width due to the depth-averaged lateral drag on both lateral walls 
$\int_{x}^{x_t} {\tau}_w dx'$. For real glaciers, buttressing may depend on spatially varying drag in 3-dimensions with 
$\tau_{bx}(x, y)$ and 
$\tau_w(x,z)$, and bed topography 
$\rho_i g H \partial_x b$. These buttressing extensions are discussed in Appendix D assuming negligible variation in glacier width along x.
Importantly, different dominant balances between the terms in (76) can lead to a diverse set of calving styles. Seasonal calving behavior observed for Greenland can occur if any of the components altering buttressing (mélange, drag forces, or other potential contributors) vary seasonally. Thus, quantifying the relative magnitude of each term in (76) can help in understanding the calving styles of glaciers.
3. Discussion
In this section, we summarize and interpret our results derived in previous sections, for an ice shelf (Table 2), MTG (Tables 3 and 4) and LTG (setting λ = 0 in Table 3 and the MS+MB column of Table 4).
Generally, the calving criteria depend on a set of two key dimensionless numbers 
$\{B, \lambda\}$ and the densities of ice and saltwater. If meltwater is present in either the surface or basal cracks, the calving criteria can depend on two more parameters 
$\tilde{h}_w, \tilde{z}_h$ and the density of meltwater. The calving regime diagram, Fig. 9a, showcases MTGs, LTGs (λ = 0) and ISs (λ = 1), such that the calving criteria can be plotted as a function of B and the dimensionless water level λ. In all cases explored, decreasing the buttressing B eventually leads to calving (dashed curves in Fig. 9a, also illustrated by the cartoons i through iv in b). The effect of the dimensionless water level λ is not as simple, as seen in the bounds for buttressing 
$B^{*}$ and B F in Tables 2–4 and shown by the dashed and solid curves in Fig. 9a, respectively.
Figure 9. Panel a displays the calving regime diagram as a function of the dimensionless buttressing B and the dimensionless water level 
$\lambda \equiv -\rho_w b / (\rho_i H)$ (= 0 for a land-terminating glacier with a flat bed; = 1 for an ice shelf). The onset of calving and crack initiation is shown by the dashed curves and bold curves, respectively, with plausible crack(s) existence living within the shaded regions. We use head height values of 
$\tilde{z}_h = \rho_i / ( 2 \rho_m )$ for DS+LMB in blue and 
$\tilde{z}_h = 3 \rho_i / ( 4 \rho_m)$ for DS+HMB in black. Panel b shows four cases of glaciers reaching the calving threshold corresponding to different locations in panel a, labeled as i to iv, with ocean saltwater shown in blue and freshwater shown in green. DS = Dry Surface, MS = Meltwater Surface, DS+SB = Dry Surface and Saltwater Basal, DS+L(H)MB = Dry Surface and Low (High)-pressure Meltwater Basal.
3.1. Dry surface crevasse atop a saltwater basal crevasse (DS+SB)
In the case of a dry surface crevasse atop a saltwater basal crevasse (DS+SB) with constant ice thickness and a flat bed, the calving criterion is 
$B^{*} = 0$ regardless of the dimensionless water level λ (dashed purple line in Fig. 9a and cartoon b-iii). Since the saltwater basal crevasse has its pressure driven by the proglacial water height, lowering the dimensionless water level λ will also decrease the pressure in the basal crevasse. Hence, decreasing the water level will decrease the basal crack depth, as crossing the solid purple line in Fig. 9a represents transitioning from dual crevasses to a dry surface crevasse without a basal crevasse (DS).
3.2. Dry surface crevasse atop a meltwater basal crevasse (DS+MB)
For a basal crevasse filled with subglacial meltwater, the water pressure in the basal crevasse can be set by the subglacial hydrology, and thus is decoupled from the water level λ. But lowering the water level λ will promote meltwater basal crevasse-driven calving. For example, Figs. 9a and A1 in Appendix A show the buttressing bounds for a dry surface crevasse atop a meltwater basal crevasse (DS+MB) in blue and black: decreasing the dimensionless water level λ while keeping buttressing B fixed can result in basal crevasse formation and calving. Thus, increasing the water level toward flotation (λ = 1) will increase the supporting force from the ocean acting on the ice front and stabilize the glacier.
The set of two curves, blue and black, in Fig. 9a shows the sensitivity of meltwater basal crevasses to the piezometric head height 
$\tilde{z}_h$. The head height in the basal crevasse analyzed in this paper depends on the subglacial hydrology, which varies spatiotemporally (e.g. Figure S1 of Harper and others Reference Harper, Bradford, Humphrey and Meierbachtol(2010)) and is beyond the scope of this paper. The lower-pressure (LMB) blue curves in Fig. 9a correspond to 
$\tilde{z}_h = \frac{1}{2} \frac{\rho_i}{\rho_m}$, while the higher-pressure (HMB) black curves assume an arbitrarily high pressure 
$\tilde{z}_h = \frac{3}{4} \frac{\rho_i}{\rho_m}$, which substantially reduces the critical stresses (increase the critical buttresssing) required for calving (black dashed curve in Fig. 9a), weakening the glacier. Thus, the sensitivity of calving to subglacial water pressure is strong. Additionally, the basal crevasse formation threshold B F strongly depends on the head height 
$\tilde{z}_h$. Basal crevasses on LTGs have been observed in Ensminger and others (Reference Ensminger, Alley, Evenson, Lawson and Larson2001), Fountain and others (Reference Fountain, Jacobel, Schlichting and Jansson2005), Harper and others (Reference Harper, Bradford, Humphrey and Meierbachtol2010) and surge-type glaciers (Rea and Evans, Reference Rea and Evans2011).
3.3. Dry (DS) or meltwater-containing surface crevasse (MS) without a basal crevasse
Finally, we consider surface crevasses, dry (DS) and with meltwater (MS), without basal crevasses. The simplest case is the LTG (λ = 0), shown by Fig. 9b-i and b-ii: a dry surface crevasse (DS) will calve at B = 0, as shown by i, while a meltwater-containing surface crevasse (MS) with 
$\tilde{h}_w = \frac{1}{2}$ will calve at 
$B^* = \frac{\rho_m}{\rho_i} \tilde{h}_w^2$, as shown by (ii). Note that the HFB model predicts glaciers to be much more vulnerable to calving compared with the Zero-Stress approximation. The Zero-Stress approximation predicts calving at 
$B^*=-1$ and 
$B^* = 2 \frac{\rho_m}{\rho_i} \tilde{h}_w-1$ for the same dry and meltwater cases described above. In terms of critical calving stresses for dry surface crevasses on a constant thickness LTG, the amount of tensile resistive stress 
$\overline{R}_{xx}$ required for calving for the HFB model is only half that of the Zero-Stress approximation, similar to previous findings for floating ISs (Buck, Reference Buck2023, Coffey and others, Reference Coffey, Lai, Wang, Buck, Surawy-Stepney and Hogg2024).
For meltwater hydrofracture-induced calving, we analyze if more tension (less buttressing) is required for calving for the same amount of meltwater 
$\tilde{h}_w$ in HFB than the Zero-Stress approximation. The two calving criteria, 
$B^{*} = \frac{\rho_m}{\rho_i}\tilde{h}_w^2$ for HFB and 
$B^{*} = 2 \frac{\rho_m}{\rho_i} \tilde{h}_w - 1$ for Zero-Stress, predict the same calving stress threshold at 
$\tilde{h}_w = 1 - \sqrt{1-\frac{\rho_i}{\rho_m}} \approx 0.72$. For calving to occur with the meltwater in a surface crevasse less than 
$72\%$ of the ice thickness, less tensile stress (more buttressing) can lead to calving in HFB than the Zero-Stress approximation.
3.4. Model limitations
Our HFB calving models have a list of assumptions, which are all also used by the original Zero-Stress approximation (Nye, Reference Nye1955), including the (1) zero material strength (Appendix C.1), (2) no elastic deformation associated with lake-induced flexure and tidal or wave perturbations (Appendix C.2), (3) constant density (Appendix C.3), (4) neglected thermomechanical effects (Appendix C.4) and (5) neglected ice rheological effects (Appendix C.5). Additionally, as mentioned in Appendix D, the analysis for MTGs in this paper requires the idealistic assumptions of constant thickness and flat bed slope. The consideration of these aforementioned effects is beyond the scope of the current study. Note that assumptions (2)-(5) are also commonly assumed by existing fracture models like LEFM (van der Veen, Reference van der Veen1998a, Reference van der Veen1998b; Jiménez and Duddu, Reference Jiménez and Duddu2018; Lai and others, Reference Lai2020), the Zero-Stress approximation (Nye, Reference Nye1955; Jezek, Reference Jezek1984; Benn and others, Reference Benn, Hulton and Mottram2007a; Nick and others, Reference Nick, van der Veen, Vieli and Benn2010), and HFB for constant-thickness ISs (Buck, Reference Buck2023; Coffey and others, Reference Coffey, Lai, Wang, Buck, Surawy-Stepney and Hogg2024). While important for crack depths, temperature dependence has been found to be negligible for the stress criteria for calving when 
$\tilde{D}=1$ in Coffey and others Reference Coffey, Lai, Wang, Buck, Surawy-Stepney and Hogg(2024), and is not considered in this paper.
3.5. Comparison with existing ideas
3.5.1. Diverse calving styles
The wide-ranging calving styles can be conceptually attributed to the fact that different terms in the HFB equation dominate in varying scenarios. The method of dominant balances is the idea that the equations may be described by the balance of the two (or more) most important terms. For example, the dominant balances of different terms in the Navier-Stokes equation can describe wide-ranging phenomena from glacial flow to hurricanes. Similarly, different combinations of the dominant terms in the HFB (3) across various glacial settings can exhibit a wide range of calving styles via different dominant balances between the horizontal hydrostatic forces acting on the crevasse wall and calving front, and various buttressing forces (9) such as the basal drag, lateral drag and mélange buttressing. This can qualitatively explain a diverse range of calving styles, as listed in the following.
HFB can capture the seasonal signature of calving through the dependence on a seasonal buttressing force. For Greenland MTG that experience a loss of buttressing in summer, through thinner mélange (Xie and others, Reference Xie, Dixon, Holland, Voytenko and Vaňková2019; Meng and others, Reference Meng2025) and reduced drag, calving occurs more frequently in the summer and thus exhibits seasonality (Zhang and others, Reference Zhang, Catania and Trugman2023; Greene and others, Reference Greene, Gardner, Wood and Cuzzone2024). Similarly, if a grounded glacier experiences saltwater intrusions or begins to float, the basal drag may be substantially reduced and push the system toward calving. While not modeled fully in this paper, large geometric effects at the ice front, such as a buoyant foot, water line melt, or undercutting, may lead to a mixed-mode, flexurally-driven or shear-driven calving style (Wagner and others, Reference Wagner2014; Reference Wagner, James, Murray and Vella2016; Slater and others, Reference Slater, Nienow, Goldberg, Cowton and Sole2017; How and others, Reference How2019; Sartore and others, Reference Sartore, Wagner, Siegfried, Pujara and Zoet2025).
On ISs, if the basal drag and mélange or sea ice buttressing are negligible, calving occurs when the lateral drag, the major source of buttressing, is lost as ice passes by and loses contact with land, such as an ice rise. This can trigger rifting near the ice rise location, causing tabular icebergs, as seen around the Roosevelt Island on the Ross Ice Shelf and the Gipps Ice Rise on the Larsen C Ice Shelf. A similar effect has been modeled at Sermeq Kujalleq (Store Glacier), Greenland (Benn and others, Reference Benn2023). However, we note that our HFB crevasse-depth model in its current form is more appropriate to describe rift initiation from vertical crevasses (Coffey and others, Reference Coffey, Lai, Wang, Buck, Surawy-Stepney and Hogg2024), as it does not describe horizontal rift propagation (Lipovsky, Reference Lipovsky2020) and its timescales (Bassis and others, Reference Bassis2007). Lastly, grounding zone saltwater intrusions by tides (Rignot and others, Reference Rignot, Ciracì, Scheuchl, Tolpekin, Wollersheim and Dow2024), causing reduced basal drag and varying water density (Wilson and others, Reference Wilson, Wells, Hewitt and Cenedese2020; Robel and others, Reference Robel, Wilson and Seroussi2022; Bradley and Hewitt, Reference Bradley and Hewitt2024), may impact the crevasse formation and calving threshold.
While each of these examples has been previously studied, the key point is that the force-balance equation can serve as a unified fundamental equation to describe a range of diverse calving styles.
3.5.2. Dependence on thickness H
One important feature of HFB, the Zero-Stress approximation, and LEFM is that the calving thresholds scale linearly with the ice thickness H (Coffey and others, Reference Coffey, Lai, Wang, Buck, Surawy-Stepney and Hogg2024). Physically speaking, a greater ice thickness results in higher lithostatic stresses, meaning the ice needs correspondingly larger tensile stresses to calve. However, this is not the case when a constant critical stress threshold for different ice thicknesses is assumed for calving. When calving is set to occur above a critical stress or thickness, we would expect runaway behavior because the ice is generally thicker upstream and the depth-averaged deviatoric stress scales with H (Haseloff and Sergienko, Reference Haseloff and Sergienko2022) (e.g. Marine Ice Cliff Instability (MICI) (Bassis and Walker, Reference Bassis and Walker2012; DeConto and Pollard, Reference DeConto and Pollard2016)). Future theoretical work implementing HFB with MICI is important for understanding ice cliff stability.
3.5.3. Comparison with existing calving laws
HFB, Zero-Stress and LEFM can all be considered as “crevasse-depth” calving laws (Benn and others, Reference Benn, Warren and Mottram2007b; Reference Benn2023; Choi and others, Reference Choi, Morlighem, Wood and Bondzio2018; Wilner and others, Reference Wilner, Morlighem and Cheng2023) in that calving occurs when the crevasse (surface plus basal) depth equals the ice thickness. Most calving laws involve tuning unmeasurable parameters, e.g. the parameter σmax in the von Mises law (Morlighem and others, Reference Morlighem2016; Choi and others, Reference Choi, Morlighem, Wood and Bondzio2018; Downs and others, Reference Downs, Brinkerhoff and Morlighem2023; Wilner and others, Reference Wilner, Morlighem and Cheng2023). Our results involve measurable and known parameters (stresses and water level), except for the meltwater-filled crevasse cases. Our analytical HFB model shows that the fundamental parameters that govern calving are the in-situ buttressing stresses B, the dimensionless water level λ and the densities of ice and water in the crevasses. If meltwater is present in the crevasses, both the surface meltwater depth 
$\tilde{h}_w$ and the meltwater head height 
$\tilde{z}_h$ will affect the calving criteria. Thus, our HFB formulation can yield diverse calving criteria, suggesting that using one universal threshold value of stress for all glacier calving is not appropriate. Multiple thresholds must be used to effectively capture the various scenarios, including but not limited to those mentioned in this paper. Further work comparing these criteria with observations and integration in numerical models with additional effects (Bassis and Ma, Reference Bassis and Ma2015; Berg and Bassis, Reference Berg and Bassis2022) will be informative. An initial comparison between the HFB theory and the observed ice-shelf rift locations is available in Coffey and others Reference Coffey, Lai, Wang, Buck, Surawy-Stepney and Hogg(2024).
Recent work (Zarrinderakht and others, Reference Zarrinderakht, Schoof and Peirce2022; Coffey and others, Reference Coffey, Lai, Wang, Buck, Surawy-Stepney and Hogg2024) shows that LEFM for a surface or basal crevasse gives stress thresholds for ice shelf calving that can be understood through torque balance, 
$\Sigma \tau = 0$. It has also been shown and argued (Yu and others, Reference Yu, Rignot, Morlighem and Seroussi2017; 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; Coffey and others, Reference Coffey, Lai, Wang, Buck, Surawy-Stepney and Hogg2024) that the hydrostatic ocean restoring force can not be neglected. In Appendix B, we show that the calving stress from the conventional LEFM solution used for isolated ice shelf basal crevasses (van der Veen, Reference van der Veen1998a) written with torque balance (Zarrinderakht and others, Reference Zarrinderakht, Schoof and Peirce2022) can be modified to consider the ocean restoring force and allow for ice shelf deflection, thus predicting the same calving stress threshold as HFB.
4. Conclusion
In this paper, we generalize the HFB fracture model (Buck, Reference Buck2023; Coffey and others, Reference Coffey, Lai, Wang, Buck, Surawy-Stepney and Hogg2024) across ISs, LTGs and MTGs. We examine six tensile crack configurations in HFB defined in Fig. 3, which lead to different calving stress thresholds (see Tables 2, 3 and 4). Our generalized HFB model yields analytical solutions for crack depths and the calving criteria. We show that in the absence of meltwater in the surface or basal cracks, the calving criteria fundamentally depend on the dimensionless buttressing force B, the dimensionless water level 
$\lambda \equiv \frac{- \rho_w b}{\rho_i H}$, as well as the ice to saltwater density ratio (constant). These parameters can either be measured or calculated. In the cases when meltwater is present in crevasses, a few more parameters play a role: dimensionless meltwater depth 
$\tilde{h}_w \equiv h_w / H$ in the surface crevasse, dimensionless head height 
$\tilde{z}_h \equiv z_h / H$ in the meltwater basal crevasse, and the ice to meltwater density ratio (constant). In other words, with a specified 
$\{B, \lambda, \tilde{h}_w, \tilde{z}_h\}$ and crack configuration, an HFB calving criteria can be obtained analytically. Our result is summarized in Tables 2, 3 and 4.
The generalized crevasse-depth calving laws using HFB have great explanatory power. In general, lower buttressing B can lead to calving in every crack configuration explored. Similarly, lowering the dimensionless water level λ has the same effect, except in the case of calving driven by saltwater basal crevasses. For LTGs with dry surface crevasses (DS), and MTGs and ISs with dry surface crevasses and saltwater-filled basal crevasses (DS+SB), the calving criterion is simply no buttressing, 
$B^{*} = 0$, regardless of the dimensionless water level λ. Thus, this result of MTGs converges to ISs and LTGs with a flat bed when λ = 1 and 0, respectively.
Our result indicates that there is generally no reason to expect a universal threshold value of calving stress for all glaciers due to each crack configuration’s (Fig. 3) varying dependence on the parameters 
$B, \lambda, \tilde{h}_w, \tilde{z}_h$. HFB can be used to model a range of calving styles in a unified framework. For example, the seasonality of calving in Greenland can be caused by a decrease of buttressing in the summer via 
$B=B^{*} = 0$. This loss of summer buttressing can be attributed to various sources, e.g. reduced basal and lateral drag, thinner mélange and their combinations. We note that the HFB calving threshold only depends on the net buttressing B and is agnostic to the buttressing loss mechanisms. Climate warming is a threat as the calving criteria is very sensitive to the surface meltwater depth 
$\tilde{h}_w$ and subglacial meltwater head height 
$\tilde{z}_h$.
Modeling of the diverse calving styles (Alley and others, Reference Alley2023; Bassis and others, Reference Bassis2024) is a challenge. Our HFB model can analytically predict transitions between six crack configurations across ISs, MTGs and LTGs. The dynamical coupling between HFB and an ice flow model to assess calving behavior during glacial retreat is a topic for future study.
Acknowledgements
Ching-Yao Lai acknowledges partial funding from NSF’s Office of Polar Programs through OPP-2344690.
Competing interests
The authors declare that they have no conflict of interest.
Appendix A. Crack Depth 
$d/H$ versus Dimensionless Water Level λ
In this Appendix, we seek to further develop intuition for the HFB theory, specifically for crack depths as functions of the dimensionless water level λ, buttressing B and the dimensionless surface crevasse meltwater depth 
$\tilde{h}_w$. Figure 9a represents a calving regime diagram as a function of dimensionless water level versus buttressing. Choosing one crack configuration, one may select a point on this calving regime diagram and have a corresponding measure of crack depths. In Fig. 8, we examine horizontal slices, or fixed values of dimensionless water level λ, of Fig. 9a for a dry surface crevasse atop a saltwater basal crevasse (DS+SB). In this appendix, we consider vertical slices, or fixed dimensionless buttressing B, of Fig. 9a with variable surface meltwater 
$\tilde{h}_w\equiv h_w / H$ as plotted in Fig. A1.
Figure A1. Transects of the dimensionless water level and buttressing calving regime diagram of Figure 9a, with dimensionless crack depths versus water level given buttressing. The left column a, c and e use a buttressing value of B = 0.15, while the right column b, d and f use B = 0.3. The first row, a and b, have dry surface crevasses 
$\tilde{h}_w\equiv h_w / H=0$, while the second and third rows, c, d and e, f, have surface crevasses with meltwater to fill 
$10\%$ or 
$50\%$ of the ice thickness, respectively. Shaded contours represent where crack configuration solutions (DS, MS, DS/MS+SB, DS/MS+LMB, DS/MS+HMB) exist. We use a similar crack configuration color key as Figure 9a.
In Fig. A1, a surface crevasse without a basal crevasse on a MTG would have depth 
$d_s^{\natural}$. If basal crevasses form, the surface crevasse associated with this dual fracture setup ds branches slightly off from this surface-crevasse-only curve, and may enable calving.
The first row, a and b, represent vertical transects of Fig. 9a at B = 0.15 and B = 0.3, respectively. We see in b that by increasing the buttressing, we shut off the low-pressure meltwater basal crevasses (DS+LMB in blue shading in a) from forming, as there is not enough tension. Furthermore, the dimensionless water level required for calving from a high pressure meltwater basal crevasse (DS+HMB in dash-dotted black) decreases, while the calving threshold for a saltwater basal crevasse (DS+SB in dash-dotted purple) does not change. Thus, lowering ocean water level stabilizes a saltwater basal crevasse yet destabilizes a meltwater basal crevasse.
Comparing each row, we see minimal change between dry surface crevasses (DS) from a to b and surface crevasses with meltwater (MS) filling 
$10\%$ of the ice thickness from c to d. However, comparing e to f where meltwater fills 
$50\%$ of the ice thickness, we see a large change for surface crevasses with meltwater (MS). By increasing the amount of water relative to the ice thickness 
$\tilde{h}_w \equiv h_w / H$, we see calving around 
$\lambda \sim 0.4$. Furthermore, the meltwater-surface-crevasse-alone crack depth 
$d_s^{\natural}$ is typically greater in e than in f, and all other crack configurations branch from this solution. Thus, whereas in b there is not enough tension for a low-pressure meltwater basal crevasse (DS+LMB in blue) to form, in e this configuration cannot form because the meltwater-containing surface crevasse has already calved at a higher dimensionless water level.
We end with two subtle points related to crack formation (dashed lines) and calving thresholds (dashed-dotted lines) for meltwater and saltwater basal crevasses. First, akin to how smaller head height zh results in a smaller window for meltwater basal crevasses (DS+MB) to exist in Fig. 9a, we see the same behavior in (MS+MB) by increasing surface meltwater 
$\tilde{h}_w$ comparing each row of Fig. A1. However, this is purely due to a change in the crack formation threshold BF; the calving threshold 
$B^*$ does not change. On the other hand, for surface crevasses atop saltwater basal crevasses (MS+SB), both crack formation BF and calving 
$B^*$ thresholds change with 
$\tilde{h}_w$, requiring higher dimensionless water level λ and thus higher saltwater pressure for the crack to form BF, and less pressure for calving 
$B^*$. Importantly, with a dry surface crevasse and saltwater basal crevasse (DS+SB), the critical dimensionless water level λ for calving remains above 1, which represents flotation. Scenarios with λ > 1, and thus calving by saltwater basal crevasse in these panels, represent an unlikely super-buoyant case (Benn and others, Reference Benn2017).
Appendix B. Torque Equilibrium Calving Argument
Recent work has established that the stress required for calving from a one horizontal dimension ice shelf basal crevasse may be described through a torque-balance argument (Zarrinderakht and others, Reference Zarrinderakht, Schoof and Peirce2022). This torque-balance argument matches with the result of LEFM (van der Veen, Reference van der Veen1998a, Tada and others, Reference Tada, Paris and Irwin2000, Lai and others, Reference Lai2020). Additionally, this torque balance has been used to approximately describe the calving stress in Mode I basal crevasse LEFM including a vertical temperature profile (Coffey and others, Reference Coffey, Lai, Wang, Buck, Surawy-Stepney and Hogg2024). However, as noted in Yu and others (Reference Yu, Rignot, Morlighem and Seroussi2017), Jiménez and Duddu (Reference Jiménez and Duddu2018), Zarrinderakht and others (Reference Zarrinderakht, Schoof and Peirce2022), Coffey and others (Reference Coffey, Lai, Wang, Buck, Surawy-Stepney and Hogg2024), an issue with this formulation is that the ice shelf base is treated as stress-free. Instead, a hydrostatic ocean would exert pressure on the ice base as it deforms and would affect the stress required for calving. Including this pressure in our torque balance, we show that the only rifting stress that can satisfy both force balance and torque balance with beam flexure explicitly modeled is 
$\overline{R}_{xx} = \frac{1}{2} \left ( 1 - \rho_i / \rho_w \right ) \rho_i g H$ or 
$B^* = 0$, the HFB result.
In static equilibrium, where forces and torques sum to zero, the relevant component of torque τy at the crevassed location 
$x=x_c$ may be written as
\begin{equation}
\tau_y = \left ( \underline{r} \wedge \underline{F} \right ) \cdot \underline{\hat{y}} = r_z F_x - r_x F_z = 0,
\end{equation}where following Zarrinderakht and others Reference Zarrinderakht, Schoof and Peirce(2022), we take 
$\underline{r}$ to be the distance vector from the crevassed location at the surface 
$\left ( x=x_c, z=s \right )$, 
$\underline{F}$ is the force vector, 
$\underline{\hat{y}}$ is the unit vector into the page, and 
$r_x, r_z$ and 
$F_x, F_z$ are the distances and forces in the horizontal and vertical, respectively. The first term, 
$r_z F_x$, is included in Zarrinderakht and others (Reference Zarrinderakht, Schoof and Peirce2022), Coffey and others (Reference Coffey, Lai, Wang, Buck, Surawy-Stepney and Hogg2024), and accounts for the forces acting along the newly-formed rift walls. The second term, 
$-r_x F_z$, is missing from previous work, and accounts for the hydrostatic restoring force of the ocean due to ice shelf base vertical displacement. Thus, the torque is not balanced purely at the rift location 
$x=x_c$, but instead is balanced some length in x of the ice shelf or newly-formed iceberg.
In the case where the ice shelf can be modeled as incompressible, (B.1) may be written as
\begin{align}
&2 \left [ \int_b^s \left ( R_{xx} - \rho_i g \left ( s - z \right ) \right ) \left ( s - z \right ) dz + \int_b^0 \rho_w g z \left ( s - z \right ) dz \right ] \nonumber \\ & - 2 \rho_w g \int_0^L w\left(x\right) x dx = 0.
\end{align}The first two terms come from 
$r_z F_x$ and are the torques acting along the newly-formed rift walls, while the last term comes from 
$-r_x F_z$ and is the force applied along the deflecting ice shelf base, with 
$w \left (x \right ) \gt 0$ indicating downward ice shelf deflection (Turcotte and Schubert, Reference Turcotte and Schubert2014). This configuration considers the basal crack xc to be in the center of a symmetric ice shelf of length 2L.
The rifting stress can be obtained by solving (B.2) for 
$\overline{R}_{xx}$. This has been done in Zarrinderakht and others (Reference Zarrinderakht, Schoof and Peirce2022), Coffey and others (Reference Coffey, Lai, Wang, Buck, Surawy-Stepney and Hogg2024), neglecting the last term in (B.2) due to the ocean restoring force,
\begin{equation}
\frac{\overline{R}_{xx}}{\frac{1}{2} \left ( 1 - \frac{\rho_i}{\rho_w} \right ) \rho_i g H} = \frac{2}{3} \left ( 2 - \frac{\rho_i}{\rho_w} \right ).
\end{equation}However, this rifting stress ignoring ocean restoring force is inconsistent with the rifting stress obtained from HFB 
$\overline{R}_{xx} = \overline{R}_{xx}^{IT} \equiv \frac{1}{2} \left ( 1 - \frac{\rho_i}{\rho_w} \right ) \rho_i g H$. Below we show how by considering the ocean restoring force (the last term in (B.2)), the torque balance would not be inconsistent with the HFB.
If we model ice as a thin elastic beam (Hetényi and Hetbenyi, Reference Hetényi and Hetbenyi1946, Turcotte and Schubert, Reference Turcotte and Schubert2014), which inherently assumes force balance in the calculation of deflection 
$w\left(x\right)$ with the bending moment at the crevassed location 
$M_c \equiv M \left ( x=x_c \right )$, we will converge to the same answer as HFB theory but with the addition of explicitly modeling bending. From Hetényi and Hetbenyi (Reference Hetényi and Hetbenyi1946), Turcotte and Schubert (Reference Turcotte and Schubert2014), the deflection profile is akin to the bending of the elastic lithosphere at an ocean trench (Turcotte and Schubert, Reference Turcotte and Schubert2014). With x = 0 placed at rift location, the deflection may be written as
\begin{equation}
w \left ( x \right ) = \frac{\alpha^2 M_c}{2D} \exp{\left ( \frac{-x}{\alpha} \right )} \left ( \left ( \frac{V_c \alpha}{M_c} + 1 \right ) \cos \frac{x}{\alpha} - \sin \frac{x}{\alpha} \right ),
\end{equation}with the flexural wavelength
\begin{equation}
\alpha \equiv \left [ \frac{4D}{\rho_w g} \right ]^{1/4},
\end{equation} where Vc is the shear force and 
$D = E H^3 / \left ( 12 \left ( 1 - \nu^2 \right ) \right )$ is the rigidity. Unlike the lithosphere problem in Turcotte and Schubert Reference Turcotte and Schubert2014, there is no shear force applied to the beam surface that pushes the beam downward/upward so 
$V_c = 0$. Thus, evaluating (B.2) with (B.4) and (B.5), and taking the limit as the length of the ice shelf L approaches infinity,
\begin{equation}
\lim_{L \rightarrow \infty} \int_0^L w \left ( x \right ) x dx = - \frac{\alpha^2}{2} \frac{\alpha^2 M_c}{2D} = - \frac{M_c}{\rho_w g},
\end{equation}and we may write the torque along the ice shelf base as
\begin{equation}
-r_x F_z = 2 M_c.
\end{equation}Substituting (B.7) into (B.2) we have
\begin{equation}
2 \left [ \int_b^s \left ( R_{xx} - \rho_i g \left ( s - z \right ) \right ) \left ( s - z \right ) dz + \int_b^0 \rho_w g z \left ( s - z \right ) dz \right ] - 2 M_c = 0,
\end{equation}with bending moment
\begin{equation}
M_c = \int_b^s \left ( R_{xx} - \rho_i g \left ( s - z \right ) \right ) \left ( s - z \right ) dz + \int_b^0 \rho_w g z \left ( s - z \right ) dz.
\end{equation}This bending moment has been defined in previous literature (e.g. Weertman (Reference Weertman1957), Buck (Reference Buck2024)) with a different z = 0; however, this does not change the final expression. What is important is the moment arm (the vertical distance from the ice surface in this problem 
$\left ( s-z \right )$) is the same as that used in Zarrinderakht and others Reference Zarrinderakht, Schoof and Peirce(2022) and in the bracket term in (B.8). Therefore, the torque-balance constraint cannot be used to solve for the depth-averaged rifting stress 
$\overline{R}_{xx}$ when the ice shelf base can deflect, as used (Zarrinderakht and others, Reference Zarrinderakht, Schoof and Peirce2022, Coffey and others, Reference Coffey, Lai, Wang, Buck, Surawy-Stepney and Hogg2024) to describe rifting from basal crevasse LEFM (van der Veen, Reference van der Veen1998a).
Instead, Euler-Bernoulli beam theory inherently assumes static equilibrium, or force and torque balance (Turcotte and Schubert, Reference Turcotte and Schubert2014). As such, the HFB rifting stress, 
$\overline{R}_{xx} = \frac{1}{2} \left ( 1- \rho_i/\rho_w \right ) \rho_i g H$ or 
$B^*=0$, is the unique rifting stress that can satisfy force and torque balance. This helps explain the agreement of HFB with the numerical simulation of basal crevasse rifting with flexure in the blue curves of Fig. 6 in Buck Reference Buck(2023). Additionally, it is clear that a vertical temperature profile would not modify this result, providing a simple explanation for the result of Coffey and others Reference Coffey, Lai, Wang, Buck, Surawy-Stepney and Hogg(2024) that the rifting stress threshold with HFB is independent of the vertical temperature profile.
Appendix C. Model Limitations
C.1. Material strength and energy
One potential limitation of our HFB model may be the lack of a material strength or fracture toughness (Lawn, Reference Lawn1993, Litwin and others, Reference Litwin, Zygielbaum, Polito, Sklar and Collins2012). We develop a scaling argument that suggests that the added force required to overcome the material strength through the entire ice thickness is negligible compared to the force required to calve a glacier with zero fracture toughness.
One way to envision the zero material strength assumption in HFB is to consider a row of touching domino tiles. There is no cohesive force between each tile and its neighbor; thus, at each tile border, there is zero strength. The force required for full separation of tiles, i.e. calving, may be written considering conservative forces only as
\begin{equation}
F_x = \frac{\partial U}{\partial x},
\end{equation} where in the absence of a fracture toughness or other sources of potential energy, the potential energy is just gravitational, 
$U = U_{G}$. This has already been shown to describe the HFB theory (Coffey and others, Reference Coffey, Lai, Wang, Buck, Surawy-Stepney and Hogg2024) and compressive buckling (Coffey and others, Reference Coffey2022). In fracture mechanics (Lawn, Reference Lawn1993, Anderson, Reference Anderson2017), the Griffith criterion for crack growth is
\begin{equation}
\frac{\partial \mathcal{U}}{\partial C} = 0,
\end{equation} where 
$\mathcal{U}$ is the total energy of the system, and C is crack length times length into the page (Lawn, Reference Lawn1993, Anderson, Reference Anderson2017). Following Lawn Reference Lawn1993, a nonzero material strength (thus nonzero fracture toughness) may be defined with a surface energy US. The surface energy release rate may be written as
\begin{equation}
\frac{\partial U_S}{\partial C} = \frac{K_{Ic}^2}{E'},
\end{equation} where KIc is the mode I fracture toughness of ice, measured as 0.15 MPa · m
$^{\frac{1}{2}}$ (Litwin and others, Reference Litwin, Zygielbaum, Polito, Sklar and Collins2012), and E ʹ is the Young’s Modulus E in plane stress, or 
$E / \left ( 1 - \nu^2 \right )$ in plane strain with Poisson’s ratio ν. To have calving, one would need to have 
$\frac{\partial }{\partial C} \left ( U_G + U_I+W_{ext} \right ) \geqslant - \frac{\partial U_S}{\partial C}$ for C increasing from initial flaw size through the full thickness, with internal energy UI, such as elastic strain energy, heat, or chemical effects and external work Wext. We have tacitly assumed kinetic energy is negligible. If there was a nonzero fracture toughness 
$K_{Ic} \gt 0$, what would be the force required to overcome this surface energy?
\begin{equation}
F_x = \frac{\partial}{\partial x} \int\frac{\partial U_S}{\partial C}dC = \frac{\partial}{\partial x} \int \frac{K_{Ic}^2}{E'} dC \sim \frac{K_{Ic}^2 H \Delta y}{E' \Delta x}.
\end{equation}For a saltwater basal crevasse and a dry surface crevasse, the force required to calve with zero fracture toughness is
\begin{equation}
F_x = \frac{\partial U_G}{\partial x} = \frac{\Delta y}{2} \left ( 1 - \frac{\rho_i}{\rho_w}\lambda^2 \right ) \rho_i g H^2.
\end{equation}The ratio of these forces gives,
\begin{equation}
\frac{2 K_{Ic}^2}{E' \Delta x \left ( 1 - \frac{\rho_i}{\rho_w}\lambda^2 \right ) \rho_i g H} \sim O \left (10^{-6} \right ),
\end{equation} where we have taken fracture toughness as above, 
$E' \sim 10^9$ Pa, 
$\Delta x = 100$ m, 
$\rho_i = 917$ kg m−3, 
$\rho_w = 1028$ kg m−3, g = 9.8 m s−2, H = 300 m, and λ = 1. The surface energy US term (related to material strength) in Griffith’s energy balance 
$\frac{\partial }{\partial C} \left ( U_G + U_I+W_{ext} \right ) \geqslant - \frac{\partial U_S}{\partial C}$ is negligible compared with the gravitational potential term UG. Thus, our simple scaling argument suggests that the force required to calve a glacier with zero fracture toughness, in the absence of buttressing, is approximately the same as that of a glacier with the fracture toughness of ice (Litwin and others, Reference Litwin, Zygielbaum, Polito, Sklar and Collins2012).
C.2. Alternate mechanisms
In this section, we summarize the alternate mechanisms that we do not consider in this paper. Our force balance predicts tensile hydrofracture-induced calving involves full-thickness fracture. However, as suggested by Weertman and others (Weertman, Reference Weertman1971, Zarrinderakht and others, Reference Zarrinderakht, Schoof and Peirce2022), the closure of surface crevasse tips may also enable drainage and hydrofracture without necessarily causing calving.
Second, we do not model the non-isostatic effects of lake loading, drainage, flexural unloading, nor calving aftereffects (MacAyeal and others, Reference MacAyeal, Scambos, Hulbe and Fahnestock2003, Scambos and others, Reference Scambos2009, Amundson and others, Reference Amundson, Fahnestock, Truffer, Brown, Lüthi and Motyka2010, Banwell and others, Reference Banwell, MacAyeal and Sergienko2013, MacAyeal and Sergienko, Reference MacAyeal and Sergienko2013). Dolines, moulins and blisters underneath the ice sheet (Moore, Reference Moore1993, Chase and others, Reference Chase, Lai and Stone2021, Lai and others, Reference Lai2021, Banwell and others, Reference Banwell, Willis, Stevens, Dell and MacAyeal2024, Hageman and others, Reference Hageman, Mejía, Duddu and Martínez-Pañeda2024) are consequences of lake loading and drainage on ISs and ice sheets. The nearly axisymmetric nature of these drainage features suggests that a full 3D stress field may be required to model drainage that can lead to calving (Banwell and others, Reference Banwell, MacAyeal and Sergienko2013, MacAyeal and Sergienko, Reference MacAyeal and Sergienko2013).
Tides and ocean waves can also elastically deform the ice and induce stresses that interact with crevasses (Freed-Brown and others, Reference Freed-Brown, Amundson, MacAyeal and Zhang2012, Nekrasov and MacAyeal, Reference Nekrasov and MacAyeal2023). While we do not model these time-dependent processes, if we allow for thin beam flexure, we show in Appendix B that a torque-balance argument will converge to our 
$B^*=0$ solution.
C.3. Constant density
We assume constant densities for both the ice and the water source. We note that firn has been implicated as important for surface crevasses (Gao and others, Reference Gao, Ghosh, Jiménez and Duddu2023, Clayton and others, Reference Clayton, Duddu, Hageman and Martínez-Pañeda2024, Meng and others, Reference Meng, Culberg and Lai2024). We leave 
$\rho \left ( z \right )$ effects as an area of possible future study.
C.4. Thermomechanical effects
If coupling force balance with the heat equation, one could potentially consider crevasse wall refreezing or melting. Future work fully accounting for these time-dependent thermomechanical processes may yield insightful results.
C.5. Complicated or unknown fracture processes
There are many possible complications with fracture mechanics and rheology (Zarrinderakht and others, Reference Zarrinderakht, Schoof and Zwinger2023), such as stress concentration, the fracture process zone length (Pralong and Funk, Reference Pralong and Funk2005), material behavior (brittle, quasibrittle, or ductile) and a potential size effect (Bažant and others, Reference Bažant, Le and Salviato2021), crack tip shielding (Zarrinderakht and others, Reference Zarrinderakht, Schoof and Peirce2024), crack nucleation given ice fabric and deformation mechanism (Frost, Reference Frost2001), flexure, shear (Clerc and others, Reference Clerc, Minchew and Behn2019, Bassis and others, Reference Bassis, Berg, Crawford and Benn2021, Needell and Holschuh, Reference Needell and Holschuh2023) and densely spaced fractures.
An assumption of HFB is that changes in traction applied to our ice geometry (Fig. 2) must be balanced by crack-depth changes. In reality, the rheology of ice permits internal viscous deformation that can accommodate changing boundary conditions.
Appendix D. Deriving Buttressing from van der Veen and Whillans’ Force Balance
The goal of this Appendix is twofold. First, we find the common form of buttressing in (10), as shown for ISs in Gudmundsson Reference Gudmundsson(2013). We show the assumptions that are required for ISs, marine- and land-terminating glaciers to arrive at (10). We assume constant thickness and flat bed slope for MTGs, and flat bed slope for LTGs. Second, we discuss extending the force balance to 3D. We conclude by discussing a challenge of applying HFB in 3D.
Depth-integrating the Stokes equations and using a traction-free surface boundary condition, the 
$\hat{x}$-component of the momentum equation was derived in Van Der Veen and Whillans Reference van der Veen and Whillans1989 (see equation 14):
\begin{align}
&\partial_x \int_b^s R_{xx} dz + \partial_y \int_b^s R_{xy} dz - \rho_i g H \partial_x s - R_{xz} \vert_{\left ( z = b \right )} \nonumber \\ &+ R_{xx} \vert_{\left ( z = b \right )} \partial_x b + R_{yx} \vert_{\left ( z = b \right )} \partial_y b = 0.
\end{align}
Following the convention set in Cuffey and Paterson Reference Cuffey and Paterson2010, we will denote the second term as −τw, and group the last 3 terms as the total basal resistance −τbx. To turn this equation into a force-balance equation per unit width, we integrate the above equation from xc to xt and use 
$H = s - b$,
\begin{equation}
H \overline{R}_{xx} \vert^{x_t}_{x_c} - \int_{x_c}^{x_t} \left ( \tau_w+\tau_{bx} \right ) dx - \rho_i g \left ( \frac{H^2}{2} \vert_{x_c}^{x_t} + \int_{x_c}^{x_t} H \partial_x b dx \right ) = 0.
\end{equation} For our local force-balance argument, the quantity of interest at xc is 
$\int_b^s \sigma_{xx} dz$. We may rearrange to solve for this,
\begin{align}
&\int_b^s \sigma_{xx} \vert^{x_c} dz = H \overline{R}_{xx} \vert^{x_c} - \frac{\rho_i g }{2} H^2 \vert^{x_c} = H \overline{R}_{xx} \vert^{x_t} - \frac{\rho_i g }{2} H^2 \vert^{x_t} \nonumber \\ & - \int_{x_c}^{x_t} \left ( \tau_w + \tau_{bx} \right ) dx - \rho_i g \int_{x_c}^{x_t} H \partial_x b dx.
\end{align}We focus here on the last two terms on the right-hand side: 1) gradients in the depth-integrated horizontal shear stresses due to the lateral wall τw and the basal drag τbx, and 2) the basal topography-induced stress. Given constant thickness, retrograde bed slope provides a force in the opposite direction of the flow, and vice versa. Typically, we expect τw and τbx to greatly contribute to normal buttressing, as these forces originate from drag applied at the lateral and bottom boundaries of the glacier.
D.1. Ice shelf buttressing
For ISs, the horizontal component of normal stress along the base 
$\rho_i g H \partial_x b$ can be calculated using the isostatic assumption 
$b = - \frac{\rho_i}{\rho_w} H$. The horizontal force per unit width at xc becomes
\begin{align}
\int_b^s \sigma_{xx} \big|^{x_c} dz
&= \int_b^s \sigma_{xx} \big|^{x_t} dz
- \int_{x_c}^{x_t} \left( \tau_w + \tau_{bx} \right) dx
+ \rho_i g \int_{x_c}^{x_t} \frac{\rho_i}{\rho_w} H \partial_x H dx
\end{align}
\begin{align}
&= \int_{b}^{s} \left [ R_{xx} - \rho_i g \left ( s - z \right ) \right ]^{x_t} dz - \int_{x_c}^{x_t} \left ( \tau_w + \tau_{bx} \right ) dx + \rho_i g \frac{\rho_i}{\rho_w} \frac{H^2}{2} \vert_{x_c}^{x_t}
\end{align}
\begin{align}
&= H_t \overline{R}_{xx} - \frac{1}{2} \rho_i g H_t^2 - \int_{x_c}^{x_t} \left ( \tau_w + \tau_{bx} \right ) dx + \rho_i g \frac{\rho_i}{\rho_w} \frac{H^2}{2} \vert_{x_c}^{x_t}
\end{align}
\begin{align}
&= H_t \overline{R}_{xx}^{B=0} - \left ( B_{M\acute{e}lange} \right ) \frac{1}{2} \left ( 1 - \frac{\rho_i}{\rho_w} \right ) \rho_i g H_c^2 - \frac{1}{2} \rho_i g H_t^2 \nonumber \\
&- \int_{x_c}^{x_t} \left ( \tau_w + \tau_{bx} \right ) dx + \rho_i g \frac{\rho_i}{\rho_w} \frac{H^2}{2} \vert_{x_c}^{x_t}
\end{align}
\begin{align}
&= \frac{1}{2} \left ( 1 - \frac{\rho_i}{\rho_w} \right ) \rho_i g H_t^2 - \left ( B_{M\acute{e}lange} \right ) \frac{1}{2} \left ( 1 - \frac{\rho_i}{\rho_w} \right ) \rho_i g H_c^2 - \frac{1}{2} \rho_i g H_t^2\nonumber \\
&- \int_{x_c}^{x_t} \left ( \tau_w + \tau_{bx} \right ) dx + \rho_i g \frac{\rho_i}{\rho_w} \frac{H^2}{2} \vert_{x_c}^{x_t}
\end{align}
\begin{align}
&= \frac{1}{2} \left ( - \frac{\rho_i}{\rho_w} \right ) \rho_i g H^2_c - \left ( B_{M\acute{e}lange} \right ) \frac{1}{2} \left ( 1 - \frac{\rho_i}{\rho_w} \right ) \rho_i g H_c^2 - \int_{x_c}^{x_t} \left ( \tau_w + \tau_{bx} \right ) dx,
\end{align} where Ht and Hc are the ice thickness at xt and xc, respectively. In the transition from (D.4c) to (D.4d), if there is no ice mélange, then 
$ B_{M\acute{e}lange}=0$. We can write the left-hand side of the equation in terms of 
$\sigma_{xx}=R_{xx}-\rho_i g (s-z)$ and rearrange (D.4f),
\begin{align}
H_c \overline{R}_{xx} \vert^{x_c} = & \frac{1}{2} \left ( 1 - \frac{\rho_i}{\rho_w} \right ) \rho_i g H^2_c - \left ( B_{M\acute{e}lange} \right ) \frac{1}{2} \left ( 1 - \frac{\rho_i}{\rho_w} \right ) \rho_i g H_c^2 \nonumber \\ & - \int_{x_c}^{x_t} \left ( \tau_w + \tau_{bx} \right ) dx.
\end{align}We define the dimensionless buttressing force as
\begin{equation}
B \equiv B_{M\acute{e}lange} + \frac{\int_{x_c}^{x_t} \left ( \tau_w + \tau_{bx} \right ) dx}{\frac{1}{2} \left ( 1 - \frac{\rho_i}{\rho_w} \right )\rho_i g H^2_c},
\end{equation}such that (D.5) can be written as
\begin{equation}
\frac{\overline{R}_{xx} \vert^{x_c}}{\frac{1}{2} \left ( 1 - \frac{\rho_i}{\rho_w} \right )\rho_i g H_c} = 1 - B.
\end{equation}If there is no drag applied to the ice shelf base or lateral margins, then B = 0; if there is buttressing then B > 0.
D.2. Marine-terminating glacier buttressing
For MTGs and LTGs, the assumption of local isostasy is modified due to the bedrock. Rearranging (D.3), we have that
\begin{align}
&\frac{\overline{R}_{xx}\vert^{x_c}}{\frac{1}{2} \left ( 1 - \lambda^2 \frac{\rho_i}{\rho_w} \right ) \rho_i g H_c} \nonumber \\ & = \frac{H_t \overline{R}_{xx} \vert^{x_t} - \frac{\rho_i g }{2} H^2 \vert^{x_t}_{x_c} - \int_{x_c}^{x_t} \left ( \tau_w + \tau_{bx} \right ) dx - \rho_i g \int_{x_c}^{x_t} H \partial_x b dx}{\frac{1}{2} \left ( 1 - \lambda^2 \frac{\rho_i}{\rho_w} \right ) \rho_i g H^2_c},
\end{align} where 
$\lambda = - \frac{\rho_w}{\rho_i} \frac{b}{H} \vert^{x_c}$ is the dimensionless water level. To achieve consistency with the ice shelf case, we can simplify the MTG case to have a constant thickness and a flat bed. This allows (D.8) to be written as
\begin{align}
&\frac{\overline{R}_{xx}\vert^{x_c}}{\frac{1}{2} \left ( 1 - \lambda^2 \frac{\rho_i}{\rho_w} \right ) \rho_i g H_c} \nonumber \\&= - B_{M\acute{e}lange} + \frac{\frac{1}{2} \left ( 1 - \lambda^2 \frac{\rho_i}{\rho_w} \right ) \rho_i g H^2_c -0 - \int_{x_c}^{x_t} \left ( \tau_w + \tau_{bx} \right ) dx - 0}{\frac{1}{2} \left ( 1 - \lambda^2 \frac{\rho_i}{\rho_w} \right ) \rho_i g H^2_c}.
\end{align}Thus, the buttressing number for the constant thickness and flat bed MTG is defined as
\begin{equation}
B \equiv 1 - \frac{\overline{R}_{xx}\vert^{x_c}}{\frac{1}{2} \left ( 1 - \lambda^2 \frac{\rho_i}{\rho_w} \right ) \rho_i g H_c}=B_{M\acute{e}lange}+ \frac{\int_{x_c}^{x_t} \left ( \tau_w + \tau_{bx} \right ) dx}{\frac{1}{2} \left ( 1 - \lambda^2\frac{\rho_i}{\rho_w} \right )\rho_i g H^2_c},
\end{equation}Following Meng and others Reference Meng(2025), a simple form for the dimensionless mélange buttressing force per unit width is
\begin{align}
B_{M\acute{e}lange}\left ( t \right ) &= \left [ \frac{1}{2} \left ( 1 - \frac{\rho_i}{\rho_w} \right ) \phi \rho_i g H_M^2 \left ( t \right ) \right ] / \left [ \frac{1}{2} \left ( 1 - \lambda^2 \frac{\rho_i}{\rho_w} \right )\rho_i g H_c \right ] \nonumber \\ & = \phi \left ( \frac{H_M\left(t\right)}{H_c} \right )^2 \frac{1 - \frac{\rho_i}{\rho_w}}{1 - \frac{\rho_i}{\rho_w} \lambda^2},
\end{align} with mélange porosity ϕ and seasonally-varying mélange thickness 
$H_M\left(t\right)$.
D.3. Land-terminating glacier buttressing
For LTGs, there is no mélange at the ice front. Force balance at the ice front guarantees that 
$H_t \overline{R}_{xx} \vert^{x_t} - \frac{\rho_i g }{2} H^2_t = 0$. Thus, (D.8) simplifies to
\begin{equation}
\frac{\overline{R}_{xx} \vert^{x_c}}{\frac{1}{2} \rho_i g H_c} = 1 - \frac{\int_{x_c}^{x_t} \left ( \tau_w + \tau_{bx} + \rho_i g H \partial_x b \right ) dx}{\frac{1}{2} \rho_i g H^2_c} = 1 - B.
\end{equation}To achieve consistency with the MTG definition of the buttressing number (D.10), we assume a flat bed. For the MTG result to converge to the LTG result, the LTG must have constant thickness; similarly, the MTG result converges to the IS case when the IS has a flat ice shelf bed elevation b.
D.4. 3D force balance
Our HFB analysis can be extended to 3 spatial dimensions, which is not included in the paper, but would involve one more integral along y of (D.3). One should account for variable bed slope and thickness, yielding
\begin{align}
&\int_{y_L}^{y_R} \int_b^s \sigma_{xx} \vert^{x_c} dz dy = \int_{y_L}^{y_R} \left [ H \overline{R}_{xx} \vert^{x_t} - \frac{\rho_i g }{2} H^2 \vert^{x_t} \right. \nonumber\\
&\quad \left.- \int_{x_c}^{x_t} \left ( \tau_{bx} + \rho_i g H \partial_x b \right ) dx \right ] dy + \int_{x_c}^{x_t} H \overline{R}_{xy} \vert_{y_L}^{y_R} dx.
\end{align}If we had an arbitrary shape, we would use (3) and arrive at the same conclusion: the horizontal forces acting on the control volume must balance.
In the application of HFB to more complicated geometries, one must determine the control volume, which may potentially align with a curved crevasse plane and may not be a priori known. This may be an interesting avenue for future research.
Open access
