A simple stress-based cliff-calving law

Over large coastal regions in Greenland and Antarctica the ice sheet calves directly into the ocean. In contrast to ice-shelf calving, an increase in calving from grounded glaciers contributes directly to sea-level rise. Ice cliffs with a glacier freeboard larger than ≈100 m are currently not observed, but it has been shown that such ice cliffs are increasingly unstable with increasing ice thickness. This cliff calving can constitute a self-amplifying ice loss mechanism that may significantly alter sea-level projections both of Greenland and Antarctica. Here we seek to derive a minimalist stress-based parametrization for cliff calving from grounded glaciers whose freeboards exceed the 100 m stability limit derived in previous studies. This will be an extension of existing calving laws for tidewater glaciers to higher ice cliffs.

To this end we compute the stress field for a glacier with a simplified two-dimensional geometry from the two-dimensional Stokes equation. First we assume a constant yield stress to derive the failure region at the glacier front from the stress field within the glacier. Secondly, we assume a constant response time of ice failure due to exceedance of the yield stress. With this strongly constraining but very simple set of assumptions we propose a cliff-calving law where the calving rate follows a power-law dependence on the freeboard of the ice with exponents between 2 and 3, depending on the relative water depth at the calving front. The critical freeboard below which the ice front is stable decreases with increasing relative water depth of the calving front. For a dry water front it is, for example, 75 m. The purpose of this study is not to provide a comprehensive calving law but to derive a particularly simple equation with a transparent and minimalist set of assumptions.

1 Introduction

Ice loss from Greenland and Antarctica is increasingly contributing to global sea-level rise (Rignot et al., 2014; Shepherd et al., 2018; WCRP Global Sea Level Budget Group, 2018). A possible additional future mass loss from these ice sheets is of crucial importance for future sea-level projections (Slangen et al., 2017; Church et al., 2013; DeConto and Pollard, 2016; Kopp et al., 2017; Mengel et al., 2016; Ritz et al., 2015; Levermann et al., 2014). Ice sheets gain mass by snowfall. The question whether they contribute to changes in sea level is determined by the question how strongly this mass addition is compensated for or overcompensated for by mass loss. Ice sheets in both Greenland and Antarctica currently show a net ice loss. Calving accounts for roughly half the ice loss of the Antarctic ice shelves, the rest is lost by basal melt (Depoorter et al., 2013). For the Greenland ice sheet, calving accounted for two-thirds of the ice loss between 2000 and 2005, the rest is due to enhanced surface melting and runoff (Rignot and Kanagaratnam, 2006). Because surface melt increased faster than glacier speed, calving accounted for one-third of the Greenland ice sheet mass loss between 2009 and 2012 (Enderlin et al., 2014). In the future the melt elevation feedback might further increase surface melt (Levermann and Winkelmann, 2016).

Tidewater glaciers calve vigorously when they are near floatation thickness, producing icebergs with a horizontal extent smaller than the ice thickness. This has been expressed in semiempirical height-above-floatation calving laws (Meier and Post, 1987; van Der Veen, 1996; Vieli et al., 2002). Calving at ice-shelf fronts or floating glacier tongues has long rest periods interrupted by the calving of large tabular ice bergs (Lazzara et al., 1999) and is preceded by the formation of deep crevasses upstream (Joughin and MacAyeal, 2005). The distinction between these two kinds of calving is not always easy because a tidewater glacier can form or lose a floating tongue; this has for example been observed at the Columbia glacier in Alaska (Walter et al., 2010).

In order to model calving not just for single glaciers but for whole ice sheets, a calving parametrization is needed. Models describing the nucleation and spreading of crevasses in ice (Pralong and Funk, 2005) are computationally very intense and difficult to apply in simulations over long timescales and large spatial dimensions. In order to parametrize calving processes, several approaches have been used.

First, calving can be described as a function of strain rate and crevasse depth. Nye (1957) first described the formation of crevasses as a result of velocity gradients: the depth of the crevasse is determined by the strain rate and overburdening pressure of the ice. Observations show that ice velocities are greater near the calving front than upstream (Meier and Post, 1987); hence, crevasses form mainly at the calving front. When crevasses are deep enough, icebergs are then separated from the glacier and calve off. Benn et al. (2007) proposed a calving law with the assumption that a glacier calves where crevasses reach the water level, Nick et al. (2010) proposed calving when surface and basal crevasses meet. These calving laws have been applied successfully in higher-order flow-line models (Nick et al., 2010) and in a 3-D Stokes model (Todd et al., 2018).

Second, a number of approaches have been taken to analyze calving processes via the stress balance. Bassis and Walker (2011) analyzed depth-averaged stresses at the calving front. Considering tensile and shear failure, they found that there is an upper limit for the thickness of stable ice cliffs: an ice cliff is only stable if the glacier's freeboard (ice thickness minus water depth) is lower than 200 m. The limit decreases to 100 m if weakening of the ice through crevasses is also considered. Krug et al. (2014) used damage and fracture mechanics to model calving. This approach, using linear elastic fracture mechanics, has recently been analyzed by Jiménez and Duddu (2018), who found that it can be applied to floating shelves but not to grounded glaciers. Morlighem et al. (2016) give a calving rate in terms of ice velocity and the von Mises stress. Recent works by Ma et al. (2017) and Benn et al. (2017) solved the 2-D full-Stokes equation at the calving front with finite element methods. Ma et al. (2017) found that while sliding glaciers calve through tensile failure, for glaciers frozen to the bed, shear failure dominates. Benn et al. (2017) used finite element models to solve the stress balance and a discrete element model to simulate fracture formation. They modeled a range of calving mechanisms including calving driven by buoyancy and melt-undercutting, but did not give parameterizations of calving rates.

Finally, Mercenier et al. (2018) analyzed tensile failure with a 2-D finite element model and derived a calving law for tidewater glaciers. They analyzed crevasse formation at the glacier terminus, determined the distance of the crevasse to the front and the time to failure until the crevasse penetrates the whole glacier and the iceberg in front of the crevasse calves off. Together this gives an equation for the calving rate as a function of water depth and ice thickness.

All these approaches agree on the basic physics of glacier calving: thicker ice at the terminus leads to higher stresses and larger calving rates. Glaciers terminating in water are stabilized by the water's back-pressure and have smaller calving rates.

The stability limit derived by Bassis and Walker (2011) lead to the formation of the marine ice cliff instability hypothesis. If cliff calving from ice cliffs whose freeboards exceed the stability limit is initiated in an overdeepened basin, e.g., in East Antarctica, it can lead to runaway cliff calving where higher ice cliffs are exposed the further the grounding line retreats, causing even larger cliff-calving rates.

Pollard et al. (2015) and DeConto and Pollard (2016) incorporated cliff calving in Antarctica projections by assuming a linear relation between freeboard exceeding the stability limit and calving rate and showed that the marine ice cliff instability can lead to much faster sea-level rise than found in previous approaches. Bassis et al. (2017) rewrote the condition that the glacier freeboard should not exceed the stability limit as a lower bound on the rate of terminus advance or equivalently an upper bound on the calving rate. More research, and especially a more physically based cliff-calving law, is needed. Studies by Ma et al. (2017), Benn et al. (2017) and Mercenier et al. (2018) were made for tidewater glaciers not exceeding the stability limit and might not be applicable to glaciers exceeding the stability limit.

In this study, we analyze stresses at the calving front by solving the 2-D Stokes equation with a finite element model in order to propose a simple cliff-calving law. The purpose of this study is not to provide a comprehensive analysis. By contrast, we seek a minimalistic set of assumptions that paths the way to a simple stress-based cliff-calving law.

Figure 1 Geometrical setup of the stress computation: two-dimensional plane, flat glacier frozen to the bedrock with a calving front at its terminus. The glacier length L is 6 times as large as the glacier height H in order to ensure that the boundary condition on the right does not significantly influence the stress field at the terminus on the left. The ice thickness is denoted H, ice thickness below the water level is D and the free-board is denoted F.

2 Stress balance near the calving front

2.1 Problem setup: 2-D Stokes equation and boundary condition

In this study we consider a plane, flat glacier of constant thickness H terminating in water of depth D in a one-dimensional (flow-line) model with horizontal coordinate x and vertical coordinate z (Fig. 1).

In order to compute the stress field near the calving front we set the glacier to be grounded (relative water depth $w\equiv D/H<\mathrm{0.9}$) and frozen to the bed. The numerical domain has a length of $L=\mathrm{6}\cdot H\gg H$. The factor 6 was chosen as a compromise to reduce computational effort while ensuring that the upstream boundary does not effect stresses at the glacier terminus. L could have been chosen to be truly “much larger” than H, but that would have required a lot of computation time without significantly benefiting the precision of the calculation. The flow-line assumption is justified, for example, in situations where the glacier is wide in comparison to its length and thickness. In these cases lateral stresses can often be neglected. The flow-line assumption is a strong constraint that neglects, for example, any buttressing effects within the ice sheet. However, the considered geometry with the width of the glacier much larger than the horizontal extent in the flow-line direction $L=\mathrm{6}\cdot H$ is internally consistent and applicable to a number of situations observed both in Greenland and Antarctica. The assumption of a flat ice thickness is justifiable on a horizontal scale of several hundred meters to a few kilometers.

The ice flow and the stresses within the ice are governed by the Stokes equations,

$$\begin{array}{}\text{(1)}& {\displaystyle }{\partial }_x{\mathit{\sigma }}_{xx}+{\partial }_z{\mathit{\sigma }}_{xz}& {\displaystyle }=\mathrm{0},\text{(2)}& {\displaystyle }{\partial }_x{\mathit{\sigma }}_{zx}+{\partial }_z{\mathit{\sigma }}_{zz}& {\displaystyle }=f,\end{array}$$

and the continuity equation,

$$\begin{array}{}\text{(3)}& {\displaystyle }\mathrm{\nabla }\cdot \mathit{u}={\partial }_x{u}_x+{\partial }_z{u}_z{\displaystyle }=\mathrm{0},\end{array}$$

with the Cauchy stress tensor σ and the gravitational force f. The Cauchy stress tensor can be split into an isotropic pressure P (also called cryostatic pressure) and the deviatoric stress tensor S, such that

$$\begin{array}{}\text{(4)}& {\mathit{\sigma }}_{ij}=-P\cdot {\mathit{\delta }}_{ij}+S_{ij},\end{array}$$

where δ_ij is the Kronecker delta. Ice rheology is assumed to be given by Glen's flow law (van der Veen, 1999),

$$\begin{array}{}\text{(5)}& {\dot{\mathit{ϵ}}}_{ij}=A{S}_{\text{e}}^{n-\mathrm{1}}{S}_{ij},\end{array}$$

with the strain rate tensor ${\dot{\mathit{ϵ}}}_{ij}=\frac{\mathrm{1}}{\mathrm{2}}\left({\partial }_i{u}_j+{\partial }_j{u}_i\right)$ and the effective stress $S_{\text{e}}=\sqrt{\frac{\mathrm{1}}{\mathrm{2}}{S}_{xx}^{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{2}}{S}_{zz}^{\mathrm{2}}+S_{xz}^{\mathrm{2}}}$.

The surface boundary is assumed to be traction-free. At the calving front boundary, we assume traction continuity to the water pressure and no traction above the water line. At the glacier bed, a no-slip boundary condition is assumed, which corresponds to a glacier frozen to its bed. No inflow is assumed at the upstream boundary.

$$\begin{array}{}\text{(6)}& {\displaystyle }& {\displaystyle }\text{Ice top:}& {\displaystyle }& {\displaystyle }\mathit{\sigma }\cdot \mathit{n}=\left(\begin{array}{c}{\mathit{\sigma }}_{xz}\\ {\mathit{\sigma }}_{zz}\end{array}\right)=\mathbf{0}\text{(7)}& {\displaystyle }& {\displaystyle }\text{Ice base:}& {\displaystyle }& {\displaystyle }\mathit{u}=\mathbf{0}\text{(8)}& {\displaystyle }& {\displaystyle }\text{Ice front:}& {\displaystyle }& {\displaystyle }\mathit{\sigma }\cdot \mathit{n}=\left(\begin{array}{c}-{\mathit{\sigma }}_{xx}\\ -{\mathit{\sigma }}_{xz}\end{array}\right)=\left\{\begin{array}{l}(-{\mathit{\rho }}_{w}gz,\mathrm{0}),z<D\\ (\mathrm{0},\mathrm{0}),z>D\end{array}\right.\text{(9)}& {\displaystyle }& {\displaystyle }\text{Upstream:}& {\displaystyle }& {\displaystyle }{u}_x=\mathrm{0}\end{array}$$

2.2 Numerical solution of the stress field

The boundary value problem was solved with the Finite Element package FEniCS (Alnæs et al., 2015) and stabilized with the Pressure Penalty method (Zhang et al., 2011). The numerical domain was divided into a regular triangular mesh with 100 vertical and 600 horizontal divisions.

Since the Stokes equations are linear in the stresses and the terminus boundary condition is linear in the ice thickness, the equations can be solved on a dimensionless domain and the stresses scaled to arbitrary ice thickness. Velocities do not scale linearly but can be obtained from the scaled stresses through the ice rheology equation. The water depth at the calving front was incorporated via the relative (dimensionless) water depth $w=D/H$.

In order to determine a suitable stress criterion for cliff calving we consider a number of commonly used stresses that have a clear physical role (Fig. 2). Generally, stresses increase with ice thickness, while the presence of water at the glacier terminus decreases the stresses and stabilizes the calving front.

The deviatoric normal stress, S_xx, corresponds to an outwards force at the calving front that has two maxima, one at the waterline and one at the foot of the terminus. The deviatoric shear stress, or Cauchy shear stress, (S_xz=σ_xz), translates to a bending moment that bends the top of the calving front forward and downward.

The different components of the deviatoric stress tensor are not invariants of the stress tensor, i.e., they depend on the coordinate system in which they are computed, and therefore they are not suitable as failure criteria. The largest principal stress,

$$\begin{array}{}\text{(10)}& {\mathit{\sigma }}_{\mathrm{1}}={\displaystyle \frac{{\mathit{\sigma }}_{xx}+{\mathit{\sigma }}_{zz}}{\mathrm{2}}}+\sqrt{{\left({\displaystyle \frac{{\mathit{\sigma }}_{xx}-{\mathit{\sigma }}_{zz}}{\mathrm{2}}}\right)}^{\mathrm{2}}+{\mathit{\sigma }}_{xz}^{\mathrm{2}}},\end{array}$$

is calculated as the largest eigenvalue of the Cauchy stress tensor and corresponds to the largest normal stress in a given point. When σ₁ is positive, it is tensile and crevasses can form. The maximum shear stress,

$$\begin{array}{}\text{(11)}& {\mathit{\tau }}_{\text{max}}=\sqrt{{\left({\displaystyle \frac{{\mathit{\sigma }}_{xx}-{\mathit{\sigma }}_{zz}}{\mathrm{2}}}\right)}^{\mathrm{2}}+{\mathit{\sigma }}_{xz}^{\mathrm{2}}},\end{array}$$

acts on a plane at an angle 45^∘ to the plane where the largest principal stress acts. It has its maximum at the foot of the calving front. The maximum shear stress can be related to brittle compressive failure (Schulson, 2001) and is therefore of particular interest for cliff failure.

The von Mises stress is the second invariant J₂ of the deviatoric stress tensor,

$$\begin{array}{}\text{(12)}& {\mathit{\sigma }}_{\text{Mises}}=\sqrt{{\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}\left(S_{xx}^{\mathrm{2}}+S_{zz}^{\mathrm{2}}+\mathrm{2}{S}_{xz}^{\mathrm{2}}\right)},\end{array}$$

and is used as a measure of deviatoric strain energy. It can also be related to material failure (Ford and Alexander, 1963) and has been used as a calving criterion by Morlighem et al. (2016). Since $S_{xx}=-S_{zz}$ due to the incompressibility of ice, the von Mises stress and the maximum shear stress differ by only a single factor: ${\mathit{\sigma }}_{\text{Mises}}=\sqrt{\mathrm{3}}{\mathit{\tau }}_{\text{max}}$.

Figure 2 Stress configurations at the calving front for different relative water depths (w=0, 0.5, 0.85) for a fixed ice thickness of 1000 m. The first column shows the deviatoric normal stress in the x direction, S_xx, the second column shows the Cauchy shear stress, σ_xz=S_xz, the third column shows the largest principal stress, σ_i, and the last column shows the maximum shear stress, τ_max.

3 Cliff failure criterion

As a first step we select a failure criterion, which then yields a failure region based on the computed stress fields. As a second step we decide on a timescale for the failure in order to derive a simple calving law.

3.1 Partial thickness failure through crevasses

Crevasses are a natural candidate for ice front failure. In the case of glaciers that are frozen to the ground, crevasses, generally, do not form from the base upward (Ma et al., 2017). Instead, surface crevasses can form in the upper part of the glacier down to the depth where the principal stress becomes compressive, i.e., attains negative values (Nye, 1957). The presence of water at the calving front reduces the stresses in the ice and decreases the depth to which surface crevasses can penetrate. Surface crevasses, generally, do not penetrate through the whole glacier thickness and so crevasses cannot be the sole cause for calving. We thus do not follow this path to determine a failure region.

Surface meltwater filling surface crevasses can increase their depth (hydrofracturing) (Weertman, 1973; Das et al., 2008; Pollard et al., 2015), but this is also not considered here. The presence of crevasses weakens the ice and is expected to enable failure even when the critical shear stress is not yet exceeded but also this is not further considered here.

3.2 Full thickness shear failure

Instead, we assume shear faulting to be the dominant process in ice-cliff failure. We could use the von Mises stress as a failure criterion instead and reach qualitatively the same result because they differ only by a factor of $\sqrt{\mathrm{3}}$.

The failure region is defined as the region close to the calving front where the maximum shear stress exceeds a critical shear stress of τ_c=1 MPa (Schulson et al., 1999; Schulson, 2001). While the specific value of the critical shear stress may be subject to uncertainties (values might be between 0.5 and 5 MPa), it is mainly a constant that will not alter the calving rate dependence on the freeboard and the water depth. The specific choice of the value is motivated by laboratory experiments and can only provide an order of magnitude of the calving rate. However, the uncertainty resulting from this choice is smaller than the uncertainty arising from the estimate of the failure time (see below).

Figure 3 Assuming Coulomb failure, the required cohesion, ${\mathit{\tau }}_{\mathrm{0}}=\sqrt{{\mathit{\mu }}^{\mathrm{2}}+\mathrm{1}}{\mathit{\tau }}_{\text{max}}-\mathit{\mu }P$, is shown for different friction parameters (μ=0, 0.3, 0.8). The failure region for a maximum cohesion of τ_max=1 MPa is encased by the black line.

3.3 Comparison to Coulomb failure

In general, brittle compressive failure happens through shear faulting (Schulson et al., 1999) and can be described with the Coulomb law (Weiss and Schulson, 2009): the shear stress τ acting on the future fault plane is resisted by material cohesion S₀ and by friction μσ with the friction coefficient μ and the normal stress across the failure plane σ. Failure happens when

$$\begin{array}{}\text{(13)}& \mathit{\tau }\ge S_{\mathrm{0}}+\mathit{\mu }\mathit{\sigma }.\end{array}$$

This expression depends on the direction of the fault plane. The failure condition can be expressed more generally in terms of the maximum shear stress τ_max and the isotropic pressure P as

$$\begin{array}{}\text{(14)}& \sqrt{{\mathit{\mu }}^{\mathrm{2}}+\mathrm{1}}{\mathit{\tau }}_{\text{max}}={\mathit{\tau }}_{\mathrm{0}}+\mathit{\mu }P,\end{array}$$

where τ₀ is another measure of cohesive strength related to S₀ (Weiss and Schulson, 2009).

Weiss and Schulson (2009) provide values of μ=0.3…0.8, depending on the temperature of the ice. Since friction increases the strength of the ice, this could stabilize rather large ice cliffs. Bassis and Walker (2011) looked at upper bounds of glacier stability with a depth-averaged shear stress for different values of μ (0.65, 0.4, 0) and a cohesion of τ₀=1 MPa. With a large friction coefficient, ice cliffs would be stable for freeboards of up to 600 m (see Fig. 3) Since this is not observed in nature, they concluded that the best model is the one without friction, which only allows freeboards of up to 200 m. Thus, with vanishing friction, the Coulomb failure criterion is equal to the maximum shear stress criterion used here.

Figure 4 Outline of the failure region for different ice thicknesses on a dimensionless domain and without water stabilizing the front (ice thickness = glacier freeboard). The background color shows the maximum shear stress on a dimensionless scale with darker areas signifying larger stress. The failure region is defined as the region close to the calving front where the maximum shear stress exceeds the critical shear stress τ_c anywhere in the ice column. The outline for H=1000 m is also shown in Fig. 3 in the top-left panel.

4 Failure region

We define the failure region as the region close to the calving front where the maximum shear stress exceeds the critical shear stress τ_c anywhere in the ice column. The failure distance L is the maximum distance of the failure region to the front and was determined for a range of ice thicknesses H and relative water depths w by solving the 2-D Stokes equation numerically and tracing the contour line where the maximum shear stress τ_max equals the critical shear stress τ_c (see Fig. 4).

For a given water depth, the failure distance L increases with the ice thickness H or the glacier freeboard $F=H-D$ (Fig. 5). For glacier freeboards smaller than approximately 100 m, the failure region vanishes: the critical shear stress is not exceeded anywhere in the ice and no shear failure takes place. This confirms results by Bassis and Walker (2011), which were derived analytically with some simplifications (see Appendix A1 for more details). The relative water depths influences the slope of the freeboard–failure distance relation: for large relative water depths, the failure distance grows more quickly with increasing freeboard. This is because, for a large relative water depth, the overall ice thickness is much larger than for a similar freeboard with a smaller relative water depth and so the failure region is larger.

Above a critical freeboard of about 1000 m (see Fig. 4 for w=0 and F=H), the failure region encompasses the whole ice thickness. Below this critical value the failure region contains only the lower part of the ice thickness, but once the lower part of the ice column fails the upper part lacks support and fails as well. The freeboard–failure distance relation has a steeper slope for large freeboards when the whole ice thickness fails. This leads to a bend at the critical freeboard, and hence the two parts require separate analytical fits. Here, we only consider values below the critical freeboard because that is the range of values most likely to occur in nature.

In Fig. 5 we provide an analytical fit with a power-law function of the form

$$\begin{array}{}\text{(15)}& {\displaystyle }L& {\displaystyle }={\left({\displaystyle \frac{F-F_{\text{c}}}{F_{\text{s}}}}\right)}^s\text{m},\text{(16)}& {\displaystyle }{F}_{\text{s}}& {\displaystyle }=\left(\mathrm{115}\cdot {\left(w-\mathrm{0.356}\right)}^{\mathrm{4}}+\mathrm{21}\right)\text{m},\text{(17)}& {\displaystyle }{F}_{\text{c}}& {\displaystyle }=\left(\mathrm{75}-w\cdot \mathrm{49}\right)\text{m},\text{(18)}& {\displaystyle }s& {\displaystyle }=\mathrm{0.17}\cdot {\mathrm{9.1}}^w+\mathrm{1.76},\end{array}$$

with $w\equiv D/H<\mathrm{0.9}$ and $F\equiv H-D=H\cdot (\mathrm{1}-w)$. At first L was fitted as a function of F for each value of w. Then the parameter functions F_s, F_c and s were fitted as functions of w.

Figure 5 shows the numerical results and the fit. Note that the fit has been optimized for relative error, so for large freeboards the fit is a little off, but it was considered more important to fit the onset of cliff calving correctly.

Figure 5 Size of shear failure region L as a function of glacier freeboard $F=H-D$ and relative water depth $w=D/H$. Numerical results are shown for smaller freeboards where the failure region does not encompass the whole ice thickness (filled dots) and for large freeboards, where the failure region contains the whole ice thickness (empty circles). A power law has been fitted to the numerical results for small freeboards (continuous line), which is given by Eq. (15). The fit has been optimized for relative error in order to get the onset of cliff calving right.

5 Failure time

There is a theory for damage evolution in ice for tensile damage (Pralong et al., 2003), from which the time to failure is derived as follows (Mercenier et al., 2018):

$$\begin{array}{}\text{(19)}& T_f={\displaystyle \frac{(\mathrm{1}-D_{\mathrm{0}}{)}^{k+r+\mathrm{1}}-(\mathrm{1}-D_{\text{c}}{)}^{k+r+\mathrm{1}}}{(k+r+\mathrm{1})B({\mathit{\sigma }}_{\mathrm{0}}-{\mathit{\sigma }}_{\text{th}}{)}^r}},\end{array}$$

with the rate factor for damage evolution B, material constants r and k, initial damage D₀, critical damage D_c, stress threshold for damage creation σ_th, and the working stress σ₀, which we assume to be the maximum shear stress τ_max. With these assumptions Eq. (19) can be written as

$$\begin{array}{}\text{(20)}& T_f=({\mathit{\sigma }}_{\mathrm{0}}-{\mathit{\sigma }}_{\text{th}}{)}^{-r}/B^{*},\end{array}$$

with σ_th=0.17 MPa, r=0.43 and $B^{*}=\mathrm{65}{\text{MPa}}^{-r}{\text{a}}^{-\mathrm{1}}$, as given in Mercenier et al. (2018). These parameters have been determined by calibrating a tensile failure calving model with data on calving rate, water depth and ice thickness for a variety of tidewater glaciers in the Arctic.

However, Eq. (20) is valid only for damages created through tensile creep. The difference between tensile and compressive damage is that under tension a single crack grows in an unstable fashion to cause failure, while in compression a large number of small cracks grow in a stable fashion until their interaction causes failure (Ashby and Sammis, 1990).

There is plenty of literature about compressive creep and failure in rocks (Brantut et al., 2013). Fatigue failure happens when a material is loaded with stresses below the failure stress and fails with a time delay due to the development of micro-cracks. There is an exponential law as well as a power law for the time to failure:

$$\begin{array}{}\text{(21)}& {\displaystyle }{t}_{\text{f}}& {\displaystyle }=t_{\mathrm{0}}\exp \left(-b{\displaystyle \frac{\mathit{\sigma }}{{\mathit{\sigma }}_{\mathrm{0}}}}\right),\text{(22)}& {\displaystyle }{t}_{\text{f}}& {\displaystyle }=t_{\mathrm{0}}^{\prime }{\left({\displaystyle \frac{\mathit{\sigma }}{{\mathit{\sigma }}_{\mathrm{0}}}}\right)}^{-b^{\prime }}.\end{array}$$

The power-law exponent is usually large, $b^{\prime }\approx \mathrm{20}$, so the power law is very similar to the exponential law. Once the major stress σ exceeds the instantaneous strength σ₀, immediate failure is assumed (t_f=0). Both time to failure relations fit the experimental data for rock well (Amitrano and Helmstetter, 2006). However, the constants depend on material properties, and there are to our knowledge no studies for time dependence of compressive creep failure in ice.

This leaves us with a dilemma: there have been no studies that determined the material properties of ice under time-dependent brittle compressive failure. Also, we cannot determine those material properties ourselves by fitting the resulting calving law to observations because, so far, cliff calving has not been observed as the major calving process in any glacier. That makes it impossible to estimate the time to failure using Eq. (21) or (22). Equation (20) and the value of its constants have been determined for tensile failure, which is microscopically very different from brittle compressive failure. So there is little reason to expect it to describe the timescale of shear failure well.

Nevertheless, we will use it as a starting point for our further analysis. For the stresses above the shear failure threshold, σ₀>1 MPa, the time to failure for tensile failure (given by Eq. 20) changes by only a factor of 2 (see Fig. 6). Hence, the calving relation can be further simplified by assuming that there is a characteristic time to failure, T_c, that is the same for all stresses and sizes of failure regions, T_c≈ 4 d. This characteristic time has been derived from parameters determined for tensile failure, so its application to shear failure comes with an uncertainty that is difficult to quantify.

Figure 6 Time to failure given by Eq. (20). For stresses above the shear failure threshold, σ₀>1 MPa, the time to failure changes only little (box).

6 Calving law

With a constant failure time, the calving rate is proportional to the size of the failure region

$$\begin{array}{}\text{(23)}& {\displaystyle }C& {\displaystyle }=C_{\mathrm{0}}\cdot {\left({\displaystyle \frac{F-F_{\text{c}}}{F_{\text{s}}}}\right)}^s,\text{(24)}& {\displaystyle }{F}_{\text{s}}& {\displaystyle }=\left(\mathrm{115}\cdot {\left(w-\mathrm{0.356}\right)}^{\mathrm{4}}+\mathrm{21}\right)\text{m},\text{(25)}& {\displaystyle }{F}_{\text{c}}& {\displaystyle }=\left(\mathrm{75}-w\cdot \mathrm{49}\right)\text{m},\text{(26)}& {\displaystyle }s& {\displaystyle }=\mathrm{0.17}\cdot {\mathrm{9.1}}^w+\mathrm{1.76},\text{(27)}& {\displaystyle }{C}_{\mathrm{0}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}\text{m}}{\mathrm{4}\text{d}}}=\mathrm{91.25}\mathrm{m}{\mathrm{a}}^{-\mathrm{1}},\end{array}$$

with $w\equiv D/H<\mathrm{0.9}$ and $F\equiv H-D=H\cdot (\mathrm{1}-w)$.
A dry cliff (w=0) reaches calving rates of C=50 km a⁻¹ at an ice thickness of $F=H\approx \mathrm{800}$ m, while an ice cliff that is close to floatation (w=0.8) reaches the same calving rate at a freeboard of F≈300 m, which corresponds to an ice thickness of H≈1500 m (see Fig. 7).

Table 1 Table of parameters in the cliff-calving relation (Eq. 23), giving the exponent s, critical freeboard F_c and scaling factor F_s for a range of relative water depth values w.

Figure 7 Cliff-calving rates C as a function of glacier freeboard $F=H-D$ and relative water depth $w=D/H$, as given by Eq. (23).

How do cliff-calving rates given by Eq. (23) compare to currently observed calving rates? A glacier enters the cliff-calving regime when its freeboard is larger than the critical freeboard F_c and the cliff-calving rate given by Eq. (23) becomes nonzero. Obviously, glaciers calve through tensile failure before and after they reach the cliff-calving regime, so we expect the overall calving rate to be larger than the cliff-calving rate, especially for glaciers that just entered the cliff-calving regime and are heavily crevassed.

Jakobshavn glacier in Greenland is one of the few glaciers that are currently in a cliff-calving mode. Jakobshavn glacier terminates in water with a depth of 800 m (Morlighem et al., 2014) and has a glacier freeboard of 100 m (Xie et al., 2018). Therefore, it can be considered to be at the beginning of the cliff-calving regime. Since the terminus is also heavily crevassed, we expect tensile calving to be the main contribution to the overall calving rate. Hence, this example can only give an upper bound on the possible cliff-calving rate.

It is difficult to determine calving rates directly. The ice flow velocity to the front of Jakobshavn is up to 12 km a⁻¹ (Joughin et al., 2012). The grounding line of Jakobshavn glacier retreats and advances seasonally about 6 km each year, but the maximum grounding line position has not changed much between 2012 and 2015 (Xie et al., 2018). Assuming a fixed grounding line, the calving rate plus frontal melt rate would equal the flow velocity. Hence, the averaged yearly calving rate is smaller than 12 km a⁻¹.

Inserting values of glacier freeboard and water depth given above into Eq. (23) gives a cliff-calving rate of C=750 m a⁻¹, which is well below the overall calving rate.

7 Discussion and conclusion

We solved the 2-D Stokes equation numerically for a flat glacier frozen to its bed in a flow-line model and investigated the stresses at the calving front.

The following four simplifications were made.

1. The model was solved in one horizontal direction, neglecting lateral shear effects. Without lateral shear effects, the result is independent of the topography of individual glaciers.

2. We assumed a basal boundary condition corresponding to a glacier frozen to its bed. Sliding was not considered.

3. The main failure mechanism was assumed to be shear faulting. We assumed brittle compressive failure according to the Coulomb law without friction stabilizing the ice cliff. Friction would allow glaciers with larger freeboards than observed to be stable.

4. A constant time to failure has been assumed.

Under these assumptions, crevasses cannot penetrate the whole glacier depth and shear failure was chosen as the main failure mechanism. The region where shear stresses exceed a critical shear stress of 1 MPa is called the failure region. The extent of this failure region, the failure distance, was determined for a range of glacier freeboards and relative water depths. For freeboards small enough for the failure region not to encompass the whole ice thickness, an analytical fit was made. Assuming a constant time to failure, a cliff-calving rate was derived. Resulting cliff-calving rates seem large compared to currently observed calving rates. Comparison with Jakobshavn glacier in Greenland shows that the cliff-calving rate is smaller than the overall calving rate; hence, we conclude that Eq. (23) probably does not overestimate cliff-calving rates.

7.1 Idealized setup vs. realistic conditions

The cliff-calving rate was derived using an idealized setup, given by the first two of the four assumptions described above. Realistic glaciers that might experience cliff calving sit in valleys where they experience lateral drag and may be sliding. The calving front may have a slope rather than a vertical cliff and there might be an undercut caused by frontal melt.

7.1.1 Sliding glaciers

First consider sliding with a constant velocity v (i.e., vanishing strain rate) for which the upstream boundary condition is an influx with velocity v, so u=v. The basal boundary conditions become u=v, w=0. Solving the Stokes' equations with these boundary conditions numerically with FeniCS gives the exact same stress fields as in the frozen case and the velocity field is simply shifted by the sliding velocity v. This is not surprising: a simple Galilean transformation takes this sliding glacier back to the frozen glacier previously considered without changing any of the physics.

In general, sliding velocities increase towards the glacier terminus. The steepest possible velocity gradient can be obtained with a free-slip basal boundary condition: we assume no influx at the upstream boundary, u=0, and at the bed we assume free slip in the horizontal direction, which only leaves a boundary condition for the vertical velocity, w=0. The basal velocity is zero at the upstream boundary and takes its maximum at the calving front. Due to this velocity gradient, the maximum shear stress is large throughout the whole numerical domain (see Fig. 8). For increasing ice thickness it becomes difficult to define a meaningful failure region because the critical shear stress is exceeded in the whole numerical domain – one must assume that the whole numerical domain will fail. Thus, in the case of a sliding glacier, the failure region is larger than in the case of a glacier frozen to its bed. Hence, the derived cliff-calving rate can serve as a lower bound for this kind of calving front.

To summarize, the derived cliff-calving law is valid for glaciers that are frozen to the bed or sliding with a constant velocity and vanishing strain rate. It serves as a lower bound on the calving rate for glaciers in which velocities increase towards the calving front.

Figure 8 Stress configurations at the calving front for different relative water depths (w=0, 0.5, 0.85) for a fixed ice thickness of 1000 m with a free-slip basal boundary condition, instead of the no-slip boundary condition used in the previous analysis (compare Fig. 2). The first column shows the deviatoric normal stress in the x direction, S_xx, the second column shows the Cauchy shear stress, σ_xz=S_xz, the third column shows the largest principal stress, σ_i, and the last column shows the maximum shear stress, τ_max. In contrast with the no-slip case, there is no definite failure region, as the maximum shear stress is large throughout the whole numerical domain.

7.1.2 Lateral drag

In order to investigate how lateral drag influences cliff calving, we will assume ice flow in a channel with a flow-line in the x direction. Ice is assumed to flow only in the x direction with a flow maximum in the middle of the channel. Since deviatoric stresses are connected to the strain rate, ${\mathit{\tau }}_{ij}=B{\dot{\mathit{ϵ}}}_e{\dot{\mathit{ϵ}}}_{ij}$, and the strain rate is given by the velocity gradients, ${\dot{\mathit{ϵ}}}_{ij}=\frac{\mathrm{1}}{\mathrm{2}}\left({\partial }_i{u}_j+{\partial }_i{u}_j\right)$, we get an additional deviatoric shear stress in the x–y plane, τ_xy. The other stress components in y vanish, ${\mathit{\tau }}_{yz}={\mathit{\tau }}_{yy}=\mathrm{0}$, because the respective velocity gradients vanish. The Cauchy stress tensor becomes

$$\begin{array}{}\text{(28)}& \mathit{\sigma }=\left(\begin{array}{ccc}P+{\mathit{\tau }}_{xx}& {\mathit{\tau }}_{xy}& {\mathit{\tau }}_{xz}\\ {\mathit{\tau }}_{xy}& P& \mathrm{0}\\ {\mathit{\tau }}_{xz}& \mathrm{0}& P-{\mathit{\tau }}_{xx}\end{array}\right).\end{array}$$

The principal stresses σ_i are defined as eigenvalues of σ, and the maximum shear stress τ_max is defined as the difference between the maximum and minimum principal stress. In 3-D, there is no simple analytical formula for the eigenvalues of a matrix, and therefore it is not feasible to get an analytical estimate on whether the introduction of nonzero τ_xy makes τ_max smaller or larger.

Assuming P(x,z), τ_xx(x,z) and τ_xz(x,z), as given by the FeniCS simulation with a constant τ_xy=1 MPa, we calculate the principal stresses and the maximum shear stress numerically. This shows that τ_max increases with increasing absolute value of τ_xy (see Fig. 9).

Hence, lateral shear increases the maximum shear, therefore increasing the size of the failure region and the cliff-calving rate. The derived cliff-calving rate can serve as a lower bound if lateral drag is present.

Figure 9 Maximum shear stress τ_max in the vicinity of the calving front in the case without lateral drag (a) and with a constant lateral drag of τ_xy=1 MPa (b).

7.1.3 Calving front slope

Other studies have shown that a calving front with a slope has significantly reduced stresses compared to a calving front with a vertical cliff (Benn et al., 2017; Mercenier et al., 2018). It is clear that a calving front slope also reduces the cliff-calving rate.

We have not analyzed this effect here because once cliff calving has been initiated, the full thickness calving probably prevents calving front slopes from forming. We aim to find a parametrization that can be implemented in ice sheet models capable of simulating the Antarctic ice sheet. These simulations are done on resolutions of several kilometers and cannot resolve calving front slopes on length scales of several tens or hundreds of meters.

7.1.4 Melt undercut

Undercut from melt would increase the stresses near the calving front (Benn et al., 2017) and hence increase the calving rate.

7.2 Uncertainties

Cliff calving is still a rather hypothetical process with a very limited scope of observations. Since there are currently no glaciers that are clearly in a cliff-calving regime, the calving rate cannot be fitted to observed calving rates. There is uncertainty in the maximum shear stress used to determine the failure distance as well as the time to failure.

Laboratory studies give a range of values between 0.5 and 5 MPa for the critical shear stress (Schulson et al., 1999; Schulson, 2001). A much larger uncertainty arises from the time to failure. There are studies that give time to failure relations and parameters for brittle compressive failure of rocks but none for ice. Time to failure of ice has only been studied for tensile failure. We use the time to failure relation used by Mercenier et al. (2018) as a first guess. Applying this time to failure for tensile failure to a process of shear failure is very uncertain. We guess that the time to failure could be up to an order of magnitude smaller or larger.

The scaling parameter C₀ in Eq. (23) should therefore be considered a free parameter. In any implementation of this cliff-calving relation, a range of values for C₀ should be tested for plausibility.

7.3 Comparison with other calving parametrizations

7.3.1 Other cliff-calving approaches

Bassis and Walker (2011) derived a stability limit for ice cliffs considering shear and tensile failure (their assumptions are analyzed further in the appendix). According to Eq. (23), cliff calving starts when the freeboard exceeds F≈75 m, this is close to the stability limit of F≈100 m given by Bassis and Walker (2011).

Pollard et al. (2015) and DeConto and Pollard (2016) implemented cliff calving in their ice sheet model by assuming a cliff-calving rate that is zero until the freeboard has reached ≈ 100 m, increases linearly up to 3 km a⁻¹ for a freeboard of about 150 m and stays constant after that. The calving relation is modified by factors representing back stress and additional wet crevasse deepening. Edwards et al. (2019) did an ensemble study with a range of values for the maximum cliff-calving rate from 0 km a⁻¹ (no cliff calving) up to 5 km a⁻¹. Depending on the scaling constant C₀, cliff-calving rates given by Eq. (23) have an equal range of magnitude but increase with a power-law dependence and have no upper bound.

Bassis et al. (2017) implemented cliff calving by requiring that ice cliffs cannot exceed the stability limit. This becomes a condition for the speed of grounding line retreat and advance. Equation (23) is easier to implement in ice sheet models because it can be implemented just like other calving parameterizations and does not need to be rewritten as a condition for the grounding line.

7.3.2 Other stress-based calving laws

Mercenier et al. (2018) derived a cliff-calving law for tidewater glaciers below the stability limit by solving the stresses in the vicinity of the front and assuming tensile failure through the formation of a large crevasse. In contrast, we assume shear failure (also called brittle compressive failure). The calving rate given by Mercenier et al. (2018) increases approximately linearly with the freeboard and has no lower bound, while the calving rate given by Eq. (23) grows with a power s(w)>1 for freeboards larger than the critical freeboard F_c(w) (see Fig. 10). Hence, we expect tensile failure to dominate for small freeboards and shear failure to dominate for large freeboards.

It is difficult to say at which glacier freeboard the tensile failure regime ends and the shear failure regime begins, not only due to uncertainty in the scaling parameter C₀. In practice, both failure modes will interact, with tensile stress damaging the ice through a few large crevasses originating from the surface of the ice and shear stress damaging the ice through a large number of small crevasses in the lower part of the cliff. This likely interaction of failure modes cannot be analyzed by assuming ice to be a continuous medium (like the approach used here and by Mercenier et al., 2018) but should be done with damage theory or a discrete element approach.

Figure 10 Comparison of the cliff-calving law given by Eq. (23) (continuous line) with the calving law for tidewater glaciers given by Mercenier et al. (2018), Eq. (22) (dotted line). Note that the cliff-calving rate could be scaled differently due to the uncertainty in C₀.

7.4 Conclusion

The calving law proposed here was derived under a number of constraining assumptions. First, it was assumed that friction plays no role in shear failure. Second, it was assumed that once the critical shear stress is exceeded, ice fails after a constant time to failure. An improved cliff-calving model might include friction and allow a stress-dependent time to failure.

If the Coulomb law with a friction component is used, the immediate failure region is smaller than in the no-friction case. Time to failure relations for compressive failure, as given by Eqs. (21) and (22), are valid for stresses below the critical shear stress. Failure is assumed to be instantaneous as soon as the critical shear stress is reached. Regions where the stress is below the failure stress would be assigned a stress-dependent failure time leading to a spatially distributed time to failure. Since friction is smaller at the top of the ice cliff, the top would fail earlier than the base, leaving a foot that would subsequently fail due to buoyant forces. There is no simple way to find a parametrization of the cliff-calving rate for these processes.

Another problem is that there are no laboratory studies on the parameters in the time to failure relations for ice. It is also not possible to calibrate the calving relation using observed calving rates because there are no glaciers currently available where cliff calving is the primary failure mechanism. Paleorecords might provide some means to calibrate cliff-calving rates as attempted in Pollard et al. (2015) and DeConto and Pollard (2016).

Paleorecords might not be constraining enough to provide a useful limit for the Antarctic sea-level contribution of the next 85 years. But even if it is difficult to constrain the rate of cliff-calving there are important qualitative consequences of a monotonously increasing cliff-calving dependence on ice thickness. The most important is the potential of a self-amplifying ice loss mechanism, which is not constrained by the reduction in calving but must be constrained by other processes. Without some kind of cliff-calving mechanism it is likely that ice sheet models are lacking an important ice loss mechanism.

Appendix A: Simplified stress balance

It is possible to solve the stress balance at the calving front analytically in a depth-averaged model with a simplifying assumption for the isotropic pressure. This has been used by Bassis and Walker (2011) and Pollard et al. (2015). It is interesting to compare this with the numerical stress field solution obtained above.

Bassis and Walker (2011) and Pollard et al. (2015) assumed the isotropic pressure is given by the gravitational pressure

$$\begin{array}{}\text{(A1)}& P(x,z)={\mathit{\rho }}_{i}g(H-z),\end{array}$$

where ρ_i is the density of ice. This assumption is actually only true over length scales that are large compared with the ice thickness and far from the ice margins (MacAyeal, 1989), which is not the case when stresses close to the calving front are analyzed. But making this assumption allows for an analytical solution of the depth-averaged stresses and does not require any ice rheology.

Together with incompressibility, which means that the trace of the strain rate disappears (${\dot{\mathit{ϵ}}}_{kk}=\mathrm{0}$) and implies $S_{xx}+S_{zz}=\mathrm{0}$, the 2-D Stokes equations become

$$\begin{array}{}\text{(A2)}& {\displaystyle }\mathrm{0}& {\displaystyle }={\displaystyle \frac{\partial S_{xx}}{\partial x}}+{\displaystyle \frac{\partial S_{xz}}{\partial z}},\text{(A3)}& {\displaystyle }\mathrm{0}& {\displaystyle }={\displaystyle \frac{\partial S_{xz}}{\partial x}}-{\displaystyle \frac{\partial S_{xx}}{\partial z}}.\end{array}$$

Assuming a traction-free surface boundary, traction-continuity at the terminus boundary and vanishing deviatoric stresses at the upstream boundary, as well as the bed boundary, a boundary value problem arises that can be solved numerically.

The resulting stresses are smaller than the stresses obtained in Sect. 2 for the 2-D Stokes equation with nonlinear ice rheology (Fig. A1). A failure region can be defined as in Sect. 3 and its size shows a very similar dependence on glacier freeboard and water depth, though it is smaller by about a factor of 3.

The biggest difference between the two approaches lies in the largest principal stress: in this simplified problem, the largest principal stress is negative in the whole ice volume; there is no region of tensile stresses, so no crevasses form. This is due to the assumption that the isotropic pressure is equal to the gravitational pressure, which is not actually the case in the vicinity of the glacier terminus.

Figure A1 Stress configurations at the calving front for different relative water depths (w=0, 0.5, 0.85) for a fixed ice thickness of 1000 m. The first column shows the deviatoric normal stress in the x direction, S_xx, the second column shows the Cauchy shear stress, σ_xz=S_xz, the third column shows the largest principal stress, σ_i, and the last column shows the maximum shear stress, τ_max.