Boundary layer models for calving marine outlet glaciers

Abstract

We consider the flow of marine-terminating outlet glaciers that are laterally confined in a channel of prescribed width. In that case, the drag exerted by the channel side walls on a floating ice shelf can reduce extensional stress at the grounding line. If ice flux through the grounding line increases with both ice thickness and extensional stress, then a longer shelf can reduce ice flux by decreasing extensional stress. Consequently, calving has an effect on flux through the grounding line by regulating the length of the shelf. In the absence of a shelf, it plays a similar role by controlling the above-flotation height of the calving cliff. Using two calving laws, one due to Nick et al. (2010) based on a model for crevasse propagation due to hydrofracture and the other simply asserting that calving occurs where the glacier ice becomes afloat, we pose and analyse a flowline model for a marine-terminating glacier by two methods: direct numerical solution and matched asymptotic expansions. The latter leads to a boundary layer formulation that predicts flux through the grounding line as a function of depth to bedrock, channel width, basal drag coefficient, and a calving parameter. By contrast with unbuttressed marine ice sheets, we find that flux can decrease with increasing depth to bedrock at the grounding line, reversing the usual stability criterion for steady grounding line location. Stable steady states can then have grounding lines located on retrograde slopes. We show how this anomalous behaviour relates to the strength of lateral versus basal drag on the grounded portion of the glacier and to the specifics of the calving law used.

Introduction

In the theory of laterally unconfined marine ice sheet flow, a standard result is that flux through the grounding line is an increasing function of bedrock depth . This leads to the conclusion that grounding lines can have stable steady states only when the ice sheet bed has sufficiently steep down-flow slope : a slight advance in grounding line position into deeper water leads to an increase in flux through the grounding line, causing the ice sheet to retreat back to its original position. Analogously, a slight retreat leads to a reduction in flux through the grounding line and a re-advance of the ice sheet.

There are a number of mechanisms that can alter the flux-to-bedrock-depth relationship. These include the appearance of “hoop stresses” in an ice shelf fringing the ice sheet (see , and , though these may require unrealistically large ice shelves), the fact that bedrock elevation can actually change due to loading and unloading of the lithosphere , thermomechanically mediated changes in basal friction and lateral drag due to geometrical confinement of the flow into a channel . The last of these is probably the most significant mechanism; when ice flows through a channel, drag can be generated by the side walls of the channel. Drag at the channel side walls of the floating ice shelf reduces the extensional stress acting at the grounding line and, therefore, reduces the grounding line flux.

and demonstrate that sidewall drag can lead to the formation of stable steady states with grounding lines on upward-sloping (or “retrograde”) beds. Both papers have channels of uniform width and fix the edge of the ice shelf, which suggests the following physics: for a steady-state grounding line on an upward-sloping bed, a slight retreat in grounding line position will cause an increase in ice thickness at the grounding line. In isolation, this would lead to increased discharge and continued ice sheet retreat. However, the retreat in grounding line position also leads to a longer shelf, which therefore experiences more lateral drag and reduces extensional stress at the grounding line. This process is known as “buttressing” and, by itself, would lead to reduced discharge and a re-advance of the ice sheet. Which of the two mechanisms dominates presumably depends on the specifics of the channel and the fixed shelf front position.

An open question is whether an evolving calving front can lead to a similar stabilization, as we are no longer guaranteed that a retreat in the grounding line position leads to the same increase in shelf length and, therefore, to the same increase in lateral drag. To answer that question conclusively, we would need a universal “calving law” that can robustly predict the location of the calving front. Such a calving law is currently not available.

We investigate how two particular calving laws that are relatively widely used in models for tidewater glaciers affect buttressing in a simplified flowline model. The ice flow model itself lacks the sophistication of models that resolve the cross-channel dimension. Instead, it relies on a parameterization of lateral drag in terms of the mean along-channel velocity . The chief advantages of the model are that it allows flux through the grounding line to be computed rapidly as a function of ice thickness through the use of a boundary layer theory and that the role of different physical mechanisms becomes comparatively easy to trace. Future work will be required to address whether our results are reproduced qualitatively by more sophisticated (and more computationally intensive) models, and we hope that this paper can motivate such work.

The rationale for the calving models used here is described in greater detail in Sect. . One calving law simply states that calving occurs at the local flotation thickness at the grounding line. The calving front is at flotation when calving occurs but no floating ice shelf ever forms. We use this calving law as a simple reference case that should give results analogous to previous work on unconfined marine ice sheet flow , where the formation of a floating shelf has no effect on flux through the grounding line. The other calving law that we investigate is the “CD” calving law due to . Simulations of outlet glaciers in Greenland with this calving law have predicted stabilization of grounding lines on areas of upward-sloping bed , suggesting that it may indeed predict a relationship between flux and bedrock depth that differs from theories for unconfined marine ice sheet flow. We do, however, stress that our aim is not an exhaustive survey of all calving models. We anticipate that the analysis presented below can be applied to other calving models, but doing so is beyond the scope of our paper.

Model

Ice flow

We consider a flowline model for a rapidly sliding, channelized outlet glacier with a parameterized representation of lateral drag. The model has the same essential ingredients as those in , , , , and . Figure  shows the physical domain. Mass accumulates over the glacier and is advected seaward by ice flow. Mass is ultimately lost by flow across the grounding line and eventual calving of icebergs. Our notation is summarized in a table given in the Supplement: let $x$ be the along-flow coordinate and $t$ time, while $u(x,t)$ and $h(x,t)$ are width-averaged ice velocity and thickness, respectively. If $b\left(x\right)$ is bed elevation and $w\left(x\right)$ the width of the outlet channel, each assumed constant in time, then we model force balance and mass conservation as $$\begin{array}{cc}{\displaystyle }& {\displaystyle }\mathrm{2}\left(\stackrel{\mathrm{‾}}{B}h\right|u_x{|}^{\mathrm{1}/n-\mathrm{1}}{u}_x{)}_x-{\stackrel{\mathrm{‾}}{B}}^{\prime }{w}^{-\mathrm{1}/n-\mathrm{1}}h|u{|}^{\mathrm{1}/n-\mathrm{1}}u-\mathit{\theta }C|u{|}^{m-\mathrm{1}}u\\ & {\displaystyle }& {\displaystyle }-{\mathit{\rho }}_{\mathrm{i}}g\left(\mathrm{1},-,(,\mathrm{1},-,\mathit{\theta },),{\mathit{\rho }}_{\mathrm{i}},/,{\mathit{\rho }}_{\mathrm{w}}\right)h(h_x+\mathit{\theta }{b}_x)=\mathrm{0}& {\displaystyle }& {\displaystyle }w{h}_t+(wuh{)}_x=w\left(a,-,(,\mathrm{1},-,\mathit{\theta },),\mathit{\mu }\right)\end{array}$$ for $\mathrm{0}<x<x_{\mathrm{c}}\left(t\right)$, where subscripts $x$ and $t$ denote partial derivatives. Here, ${\mathit{\rho }}_{\mathrm{i}}$ and ${\mathit{\rho }}_{\mathrm{w}}$ are the densities of ice and water, respectively, and $g$ is acceleration due to gravity, while $a$ is surface mass balance and $\mathit{\mu }$ is the melt rate at the base of the floating ice. The indicator function $\mathit{\theta }$ is given by $${\displaystyle }\mathit{\theta }\left(x\right)=\left\{\begin{array}{ll}\mathrm{1}& \text{if}{\mathit{\rho }}_{\mathrm{i}}h\left(x\right)\ge -{\mathit{\rho }}_{\mathrm{w}}b\left(x\right),\\ \mathrm{0}& \text{otherwise};\end{array}\right.$$ in other words, $\mathit{\theta }=\mathrm{1}$ if and only if the ice thickness is above flotation and the glacier is grounded.

Note that we have included the melt rate $\mathit{\mu }$ in Eq. (1) for completeness only. While the numerical code included with the Supplement permits computations of steady-state solutions in Sect.  and of the boundary layer problem described in Sect.  with a non-zero basal melt rate under floating ice, a full exploration of the extended parameter space would make this paper unmanageable. In what follows, we set $\mathit{\mu }=\mathrm{0}$ throughout, and we will address the effect of incorporating non-zero $\mathit{\mu }$ in a separate paper.

The parameters $\stackrel{\mathrm{‾}}{B}$ and $n$ are the usual parameters in the Glen's law rheology for ice . We neglect any complications associated with the dependence of ice viscosity on temperature or moisture content and treat $\stackrel{\mathrm{‾}}{B}$ as well as $n$ as constant. $C$ is a drag coefficient in a power-law basal friction law, with $m$ being the corresponding exponent e.g.; on theoretical grounds, it is often assumed that $m=\mathrm{1}/n$. Note that other friction laws have also been considered in boundary layer models for unconfined ice sheets . The basal drag term only applies where ice is grounded, corresponding to $\mathit{\theta }=\mathrm{1}$. For simplicity, we neglect the possibility that $C$ may depend on additional degrees of freedom such as basal water pressure or temperature.

The second term ${\stackrel{\mathrm{‾}}{B}}^{\prime }{w}^{-\mathrm{1}/n-\mathrm{1}}h|u{|}^{\mathrm{1}/n-\mathrm{1}}u$ is a parameterization of lateral drag, with ${\stackrel{\mathrm{‾}}{B}}^{\prime }$ being another constant. We assume that lateral drag is exerted on both grounded and floating ice. A more sophisticated treatment of lateral drag would require a domain with two horizontal dimensions and an additional equation representing force balance in the direction perpendicular to the channel axis e.g..

It is worth noting however that Eq. () becomes accurate in one of two mutually incompatible limits: (i) a wide channel where lateral drag vanishes (this is the one-dimensional flow case previously studied in e.g. ) or (ii) a long and narrow glacier where extensional stress is insignificant and there is no significant flow transverse to the channel axis. By “extensional stress”, we mean the non-cryostatic part of normal stress on a vertical surface placed across the flow, that is, $\mathrm{2}\stackrel{\mathrm{‾}}{B}|u_x{|}^{\mathrm{1}/n-\mathrm{1}}{u}_x$. In the narrow-channel case, assuming no slip at the channel side walls, Eq. (1) predicts the correct width-averaged velocity if we put (see also, e.g., for details on flows dominated by lateral shear) $${\stackrel{\mathrm{‾}}{B}}^{\prime }=(n+\mathrm{2}{)}^{\mathrm{1}/n}{\mathrm{2}}^{(n+\mathrm{1})/n}\stackrel{\mathrm{‾}}{B};$$ smaller values of ${\stackrel{\mathrm{‾}}{B}}^{\prime }$ can be justified if there is actual sliding at the lateral margins of the ice. We use Eq. () even when neither of the two limits above apply. As discussed in , this is a simplification that works reasonably well and allows at least semi-analytical progress to be made. The simplicity of the model has also led to a large number of authors adopting versions of it. We proceed in that spirit, analysing the model at face value.

Schematic of the model domain and variables used. Not shown is the lateral dimension: the glacier occupies a channel of width $w$, potentially variable along the glacier flowline. The upper ice surface shown is at an elevation $s=(\mathrm{1}-(\mathrm{1}-\mathit{\theta }){\mathit{\rho }}_{\mathrm{i}}/{\mathit{\rho }}_{\mathrm{w}})h+\mathit{\theta }b$, and the lower surface at $l=s-h$. The lower two diagrams are enlargements of the vicinity of a grounded terminus (left) and a floating terminus (right), illustrating the meaning of calving front position $x_{\mathrm{c}}$, grounding line position $x_{\mathrm{g}}$, calving front height $h_{\mathrm{c}}$, grounding line thickness $h_{\mathrm{g}}$ and flotation thickness $h_{\mathrm{f}}$. The dashed line indicates where the upper ice surface would need to be at a given position $x$ in order for the ice to be just about to float.

[Figure omitted. See PDF]

We denote the glacier terminus by $x_{\mathrm{c}}\left(t\right)$; this is the location where ice cover ends. Since $x_{\mathrm{c}}$ is a free boundary, two boundary conditions are required. One is needed to close the elliptic problem Eq. () and another to determine the evolution of $x_{\mathrm{c}}$. The former is a condition on extensional stress at the ice front e.g.Appendix B: $$\begin{array}{cc}{\displaystyle }& {\displaystyle }\mathrm{2}\stackrel{\mathrm{‾}}{B}h|u_x{|}^{\mathrm{1}/n-\mathrm{1}}{u}_x={\mathit{\rho }}_{\mathrm{i}}\left(\mathrm{1},-,(,\mathrm{1},-,\mathit{\theta },),{\mathit{\rho }}_{\mathrm{i}},/,{\mathit{\rho }}_{\mathrm{w}}\right)h^{\mathrm{2}}/\mathrm{2}-\mathit{\theta }{\mathit{\rho }}_{\mathrm{w}}g{b}^{\mathrm{2}}/\mathrm{2}\\ & {\displaystyle }& {\displaystyle }\text{at}x=x_{\mathrm{c}}\left(t\right).\end{array}$$ We take the second condition to be a calving law. While a stress condition is sufficient to close the force balance model Eq. (), a calving model can be understood as fixing the free-boundary location. The next section describes the different choices of calving laws used here.

Calving model

The process of calving remains relatively poorly understood, but several calving laws have been developed on theoretical grounds. Our aim is to illustrate how different calving laws can lead to qualitatively, rather than quantitatively, different dynamics in the outlet glacier. We consider two possible calving laws. The first is the CD model due to , and the second is a calving law that states that ice breaks off when the glacier reaches its flotation thickness. To streamline our notation, we refer to the latter as the “FL” calving law.

Calving laws. In all three panels, a solid line refers to the CD calving law, and a dashed line to the FL law. The grey shaded region refers to parts of parameter space where the calving front is grounded for the CD law, and white background to a floating calving front. $r=\mathrm{0.9}$ throughout. (a) The calving-front-thickness-to-bedrock-depth ratio $\mathit{\varphi }$ as a function of normalized water depth $-d_{\mathrm{w}}/b$. For the FL law, the ratio is simply $r^{-\mathrm{1}}$. (b) Normalized calving front thickness $h_{\mathrm{c}}/d_{\mathrm{w}}$ as a function of normalized flotation thickness $h_{\mathrm{f}}/d_{\mathrm{w}}=-b/\left(r{d}_{\mathrm{w}}\right)$ at the grounding line. (c) Normalized grounding line thickness $h_{\mathrm{g}}/d_{\mathrm{w}}$ as a function of normalized flotation thickness $h_{\mathrm{f}}/d_{\mathrm{w}}$.

[Figure omitted. See PDF]

The CD model works based on the assumption that water in surface crevasses affects the depth to which those crevasses can penetrate. When they penetrate deeply enough to connect with basal crevasses, calving occurs. When they do not, there is no calving and the ice front simply moves at the velocity of the ice. Algebraic manipulation of the CD model shows that connections with basal crevasses occur instantly in the model when ice thickness is at (or below) a value $h_{\mathrm{c}}$. In other words, the evolution of the calving front $x_{\mathrm{c}}$ satisfies $$\begin{array}{}& {\displaystyle }& {\displaystyle }\text{either}h=h_{\mathrm{c}}\text{at}x=x_{\mathrm{c}}\text{if }{\dot{x}}_{\mathrm{c}}\le u\left(x_{\mathrm{c}}\right)& {\displaystyle }& {\displaystyle }\text{or }{\dot{x}}_{\mathrm{c}}=u\text{at}x=x_{\mathrm{c}}\text{if}h>h_{\mathrm{c}},\end{array}$$ where the dot indicates differentiation with respect to time. Note that the domain lies to the left of $x_{\mathrm{c}}$, so ${\dot{x}}_{\mathrm{c}}<u$ implies that ice is removed by calving. $h_{\mathrm{c}}$ itself can be written as a function of the crevasse water depth $d_{\mathrm{w}}$ and of local bedrock depth $-b$. This function can be expressed as $$\begin{array}{cc}& {\displaystyle }& {\displaystyle }{h}_{\mathrm{c}}=-b\mathit{\varphi }\left(-,d_{\mathrm{w}},/,b\right),{\displaystyle }& {\displaystyle }\text{where}\mathit{\varphi }(-d_w/b)=\\ & {\displaystyle }& {\displaystyle }\left\{\begin{array}{ll}\mathrm{2}({\mathit{\rho }}_w/{\mathit{\rho }}_i)(-d_w/b)& \text{when}-d_w/b<\mathrm{1}/\mathrm{2},\\ \mathit{\nu }+\sqrt{{\mathit{\nu }}^{\mathrm{2}}-({\mathit{\rho }}_w/{\mathit{\rho }}_i)}& \text{when}-d_w/b\ge \mathrm{1}/\mathrm{2},\end{array}\right.\end{array}$$ and $\mathit{\nu }$ is $${\displaystyle }\mathit{\nu }=\mathrm{1}+\left[(,{\mathit{\rho }}_{\mathrm{w}},/,{\mathit{\rho }}_{\mathrm{i}},),-,\mathrm{1}\right](-d_{\mathrm{w}}/b).$$ The function $\mathit{\varphi }$ is the ratio of calving front thickness to depth to bedrock; the form of $\mathit{\varphi }$ is illustrated in Fig. a.

As an alternative to the CD model, in which the function $\mathit{\varphi }$ is defined through Eq. (), we consider the FL law, in which the glacier calves “at flotation”. This means that calving occurs when $h=-({\mathit{\rho }}_{\mathrm{w}}/{\mathit{\rho }}_{\mathrm{i}})b$ at the calving front but not when $h$ is larger. This is easy to incorporate into the calving framework Eqs. ()–() above: we simply have to replace the definition of $\mathit{\varphi }$ in Eq. () with the simpler $$\mathit{\varphi }={\mathit{\rho }}_{\mathrm{w}}/{\mathit{\rho }}_{\mathrm{i}}.$$ This condition is effectively what applies in previous work on marine ice sheet flow without sidewall drag as considered in, e.g., .

Note that do not formulate their calving law directly in the form Eqs. ()–(); a derivation of the thickness condition Eq. () based on their formulation is given in the Supplement. There are two cases in Eq. (): $-d_{\mathrm{w}}/b<\mathrm{1}/\mathrm{2}$ corresponds to a floating terminus, and $-d_{\mathrm{w}}/b\ge \mathrm{1}/\mathrm{2}$ to a grounded terminus. Note that $\mathit{\varphi }$ and $h_{\mathrm{c}}$ are continuous (in fact, continuously differentiable, as shown in Fig. ) at $-d_{\mathrm{w}}/b=\mathrm{1}/\mathrm{2}$, where calving occurs at the critical floatation thickness $h_{\mathrm{c}}=({\mathit{\rho }}_{\mathrm{w}}/{\mathit{\rho }}_{\mathrm{i}})b$. Once the calving front is afloat, $h_{\mathrm{c}}$ no longer depends on bedrock depth: when $-d_{\mathrm{w}}/b<\mathrm{1}/\mathrm{2}$, we simply have $h_{\mathrm{c}}=\mathrm{2}({\mathit{\rho }}_{\mathrm{w}}/{\mathit{\rho }}_{\mathrm{i}})d_{\mathrm{w}}$. In other words, for a fixed water depth parameter $d_{\mathrm{w}}$ and sufficiently large bedrock depths the CD model is actually a calving law that simply states that ice breaks off a floating glacier shelf at a critical thickness that is determined purely by the parameter $d_{\mathrm{w}}$ (Fig. b).

To complete the notation for our model, we also define the grounding line position $x=x_{\mathrm{g}}$. For a glacier with a floating shelf, this is the point $x=x_{\mathrm{g}}\left(t\right)$ at which $\mathit{\theta }$ changes discontinuously: $$h\left(x_{\mathrm{g}}\right(t),t)=-({\mathit{\rho }}_{\mathrm{i}}/{\mathit{\rho }}_{\mathrm{w}})b\left(x_{\mathrm{g}}\right(t\left)\right).$$ For the FL calving law, there is no floating shelf; the grounding line and calving front coincide. The CD calving model goes further, permitting cases where there is no floating ice shelf and the flotation condition $h=-({\mathit{\rho }}_{\mathrm{i}}/{\mathit{\rho }}_{\mathrm{w}})b$ is attained nowhere along the glacier, including at the calving front. To keep our terminology as simple as possible, we identify the grounding line in that case with the terminus location: $$x_{\mathrm{g}}=x_{\mathrm{c}}.$$ For later convenience, we also define the ice thickness $h_{\mathrm{g}}$ at the grounding line and the flotation thickness $h_{\mathrm{f}}$ at the grounding line through $$h_{\mathrm{g}}=h\left(x_{\mathrm{g}}\right(t),t),h_{\mathrm{f}}=-({\mathit{\rho }}_{\mathrm{i}}/{\mathit{\rho }}_{\mathrm{w}})b\left(x_{\mathrm{g}}\right(t\left)\right).$$ For glaciers with a floating shelf, we always have $h_{\mathrm{g}}=h_{\mathrm{f}}$; for that reason, existing theories for marine ice sheets generally do not make a distinction between $h_{\mathrm{g}}$ and $h_{\mathrm{f}}$ e.g.. The distinction becomes relevant when there is no floating ice shelf, in which case we only have $h_{\mathrm{g}}\ge h_{\mathrm{f}}$. We will use $h_{\mathrm{f}}$ frequently below, as it is a simple function of bedrock depth at the grounding line and therefore determined purely by the geometry of the glacier channel. $h_{\mathrm{g}}$ additionally depends on the calving process as shown in Fig. c.

In our view, the CD model is a cartoon version of the linear elastic fracture mechanics explored in by and . These papers consider the “mode 1” propagation of vertical cracks into ice under tensile (extensional) stresses. This is done by computing stress intensity factors at the crack tips from known Green's functions for parallel-sided elastic slabs with cracks penetrating from the upper or lower surfaces, accounting for the pressure exerted by water in the cracks, and applying a fracture toughness criterion. The CD model by contrast assumes that extensional stress increases with depth in the ice in a linear, cryostatic fashion. The model then computes crevasse penetration as being the distance from the upper and lower surfaces at which that extensional stress becomes sufficiently negative (that is, sufficiently compressive) to overcome the pressure exerted by water at the same depth. The CD model therefore does not compute stress with the same level of sophistication as the papers by and but follows the same basic approach of computing crevasse propagation based on a known ice geometry, extensional stress and crevasse water pressure, and it has the advantage of tractability.

The basic method in in principle allows for a constraint to be computed that links ice thickness, applied extensional stress, crevasse water level and fracture toughness at the moment that surface and basal crevasses together first penetrate through the entire ice thickness. Given that extensional stress is a function of ice thickness through Eq. (), this constraint could be converted into a criterion for the thickness $h_{\mathrm{c}}$ at which calving occurs, giving a more sophisticated version of the CD model. However, the papers by van der Veen do not deal with the case in which both surface and basal crevasses are present and interact with each other (so the relevant Green's functions are not given), and he does not explicitly compute a condition for calving that could be put in the form Eq. (). As a result, we confine ourselves to the simpler CD model here.

Note that there is one inconsistency in the calving law at small $-b/d_{\mathrm{w}}$: here the CD law predicts values of $h_{\mathrm{c}}<d_{\mathrm{w}}$ (see Fig. b). Obviously, this implies a greater water level in surface crevasses than the ice thickness, which is physically impossible. We are led to conclude that, for small enough $-b/d_{\mathrm{w}}$ (which will later correspond to small enough values of the parameter $\mathrm{\Lambda }$, defined in Sect. ), the calving model continues to lead to computable results but breaks down physically. Note that having a flotation height (equal to $-{\mathit{\rho }}_{\mathrm{w}}b/{\mathit{\rho }}_{\mathrm{i}}$) that is smaller than $d_{\mathrm{w}}$ is not problematic: the calving front thickness can be greater than the flotation height and therefore allow for such large water depths.

A second practical pitfall of the CD model is that it predicts no calving at all if $d_{\mathrm{w}}=\mathrm{0}$ and surface crevasses are free of water. It is possible that this is an artefact of the simple representation of stress in the CD model, where the tensile stress driving crevasse propagation is assumed to have the same dependence on depth below the ice surface regardless of whether a crevasse is present or not. In reality, the formation of crevasses that penetrate through a significant fraction of the ice shelf leads to extensional stress becoming more concentrated around the crack tips than for shallow crevasses (see for instance Fig. 4 of ). This represents a positive feedback on crack propagation and could lead to calving even for the case of water-free surface crevasses see also. In addition, the stress fields considered by , and are relatively simple and apply only at distances significantly greater than a single ice thickness from the calving front. In the calving of shorter, taller icebergs, torques near the calving front may allow calving when purely extensional stresses experienced further upstream do not.

More recently, others have extended the linear elastic fracture mechanics approach of and to include effects such as the role of distributed damage due to the formation of microcracks in initiating crevasse formation, the blunting of cracked tips due to viscous deformation and the presence of significant torques near the calving front . The complexity of these processes however makes them difficult to parameterize in a model that does not resolve the scale of individual crevasses, and we do not consider them here.

The CD calving model, along with the work of and , is based on tensile failure. We can contrast this with the shear failure model of (see also , and ). The CD model requires $d_{\mathrm{w}}>\mathrm{0}$ and predicts calving for any $h$ below the value given by Eq. (), instantaneously removing all parts of the glacier shelf that are too thin. By contrast, the shear failure model of predicts that calving will start at a critical calving front thickness and not occur below that thickness, so the inequality in Eq. () would need to be reversed. It also predicts that, once initiated, the calving front will continue to fracture as it moves into thicker ice inland. This is the basis of the catastrophic calving cliff instability mechanism for marine ice sheet collapse advocated by but cannot be captured by an analogue of Eq. (). It is clear that ice sheets whose calving cliff is larger than the critical thickness for shear failure simply cannot persist: they are guaranteed to disintegrate completely or to stabilize in some shape where the calving front thickness is below the critical thickness for shear failure, and the shear failure model by itself does not provide a timescale for that disintegration. We exclude such shear failure from consideration here and focus purely on the CD calving model.

Even when taking the CD model at face value, as we do here, the sensitivity to the parameter $d_{\mathrm{w}}$ remains problematic. In fact, one of our results below will be that flux through the grounding line is more sensitive to $d_{\mathrm{w}}$ than to any other model parameter. At present, we do not have a surface hydrology model that can predict $d_{\mathrm{w}}$. It is plausible that a future hydrology model could compute a water table height near the calving cliff ($d_{\mathrm{s}}-d_{\mathrm{w}}$ if measured relative to the local ice surface, where $d_{\mathrm{s}}$ is the depth of surface crevasses as discussed in the Supplement) rather than using $d_{\mathrm{w}}$ itself. Such a model would likely be based on drainage being driven by gradients in hydraulic head, but this awaits future development. We persist with the basic model, treating $d_{\mathrm{w}}$ as given.

Solution of the model

Non-dimensionalization

In the remainder of this paper, we will consider the problem Eq. (1) in dimensionless form. The purpose of doing so is two-fold. Non-dimensionalization (i) reduces the number of free parameters and (ii) allows systematic approximations based on the relatively small size of some dimensionless parameters. We assume that we know scales $\left[a\right]$ for accumulation rate and $\left[x\right]$ and $\left[w\right]$ for glacier length and width, respectively. We choose scales $\left[u\right]$, $\left[h\right]$ and $\left[t\right]$ based on the balances $$\begin{array}{cc}{\displaystyle }& {\displaystyle }{\stackrel{\mathrm{‾}}{B}}^{\prime }\left[w{]}^{-\mathrm{1}/n}\right[h\left]\right[u{]}^{\mathrm{1}/n}={\mathit{\rho }}_{\mathrm{i}}g[h{]}^{\mathrm{2}}/[x],\\ {\displaystyle }& {\displaystyle }\left[u\right]\left[h\right]=\left[a\right]\left[x\right],\left[u\right]\left[t\right]=\left[x\right].\end{array}$$ We define dimensionless variables as $u=\left[u\right]u^{*}$, $h=\left[h\right]h^{*}$, $x=\left[x\right]x^{*}$, $t=\left[t\right]t^{*}$ and $x_{\mathrm{c}}=\left[x\right]x_{\mathrm{c}}^{*}$, and we also put $$\begin{array}{cc}{\displaystyle }& {\displaystyle }\mathit{\epsilon }={\displaystyle \frac{\stackrel{\mathrm{‾}}{B}[w{]}^{\mathrm{1}/n+\mathrm{1}}}{\mathrm{2}{\stackrel{\mathrm{‾}}{B}}^{\prime }[x{]}^{\mathrm{1}/n+\mathrm{1}}}},\mathit{\gamma }={\displaystyle \frac{C\left[w{]}^{\mathrm{1}/n+\mathrm{1}}\right[u{]}^{m-\mathrm{1}/n}}{B^{\prime }\left[h\right]}},\\ & {\displaystyle }& {\displaystyle }\mathit{\lambda }={\displaystyle \frac{d_{\mathrm{w}}}{\left[h\right]}},r={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{i}}}{{\mathit{\rho }}_{\mathrm{w}}}},& {\displaystyle }& {\displaystyle }{a}^{*}={\displaystyle \frac{a}{\left[a\right]}},w^{*}={\displaystyle \frac{w}{\left[w\right]}},b^{*}={\displaystyle \frac{b}{\left[h\right]}}.\end{array}$$ Dropping asterisks on the dimensionless variables immediately, we obtain $$\begin{array}{cc}{\displaystyle }& {\displaystyle }\mathrm{4}\mathit{\epsilon }\left(h\right|u_x{|}^{\mathrm{1}/n-\mathrm{1}}{u}_x{)}_x-w^{-\mathrm{1}/n-\mathrm{1}}h|u{|}^{\mathrm{1}/n-\mathrm{1}}u-\mathit{\gamma }\mathit{\theta }|u{|}^{m-\mathrm{1}}u\\ & {\displaystyle }& {\displaystyle }-(\mathrm{1}-r+\mathit{\theta }r)h(h_x+\mathit{\theta }{b}_x)=\mathrm{0},& {\displaystyle }& {\displaystyle }w{h}_t+(wuh{)}_x=wa\end{array}$$ for $\mathrm{0}<x<x_{\mathrm{c}}\left(t\right)$, with $\mathit{\theta }$ being the indicator function for flotation $$\mathit{\theta }=\mathrm{1}\text{if }rh\ge -b,\mathit{\theta }=\mathrm{0}\text{otherwise}$$ and the boundary conditions at the terminus being $$\begin{array}{cc}{\displaystyle }& {\displaystyle }\mathrm{4}\mathit{\epsilon }h|u_x{|}^{\mathrm{1}/n-\mathrm{1}}{u}_x=(\mathrm{1}-r)h^{\mathrm{2}}/\mathrm{2}+\mathit{\theta }\left(r^{\mathrm{2}},h^{\mathrm{2}},-,b^{\mathrm{2}}\right)/(\mathrm{2}r)\\ & {\displaystyle }& {\displaystyle }\text{at}x=x_{\mathrm{c}}\left(t\right),\end{array}$$ and either $$\begin{array}{}& {\displaystyle }& {\displaystyle }h=-b\mathit{\varphi }(-\mathit{\lambda }{b}^{-\mathrm{1}})\text{at}x=x_{\mathrm{c}}\text{ if }{\dot{x}}_{\mathrm{c}}\le u\left(x_{\mathrm{c}}\right)& {\displaystyle }& {\displaystyle }\text{or }{\dot{x}}_{\mathrm{c}}=u\text{at }x=x_{\mathrm{c}}\text{ if }h>-b\mathit{\varphi }\left(-,d_{\mathrm{w}},/,b\right),\end{array}$$ with $\mathit{\varphi }$ given by Eq. () for the CD calving model, or by Eq. () (which states that $\mathit{\varphi }\equiv r^{-\mathrm{1}}$) for the FL calving law.

Direct numerical solution

The system Eq. (4) can be solved numerically as posed. In this paper, we focus on solutions of the steady-state version of the problem by a shooting method, which provides a straightforward alternative to a solution by more established time-stepping methods. As our method has not been used previously in this context, we sketch it here for completeness; results are presented at the end of this section and in Figs.  and .

We can write the steady-state problem as a four-dimensional, first-order autonomous system of differential equations if, in addition to $h$, we define the phase space variables $q$, $\mathit{\sigma }$ and $\mathit{\chi }$ through $$q=uhw,\mathit{\sigma }=|u_x{|}^{\mathrm{1}/n-\mathrm{1}}{u}_x,\mathit{\chi }=x.$$ For technical reasons associated with singular behaviour of the steady-state problem near ice divides, we also define a new independent variable $\mathit{\eta }$ through $${\displaystyle }{x}_{\mathit{\eta }}=q$$ to obtain a first-order system of differential equations from Eq. (4): $$\begin{array}{cc}{\displaystyle }& {\displaystyle }{h}_{\mathit{\eta }}=-h^{\mathrm{2}}w\left(\mathit{\chi }\right)|\mathit{\sigma }{|}^{n-\mathrm{1}}\mathit{\sigma }+hw(\mathit{\chi }\left)a\right(\mathit{\chi })\\ & {\displaystyle }& {\displaystyle }-hq{w}^{\prime }\left(\mathit{\chi }\right)w(\mathit{\chi }{)}^{-\mathrm{1}},{\displaystyle }& {\displaystyle }{\mathit{\sigma }}_{\mathit{\eta }}=\left(\mathrm{4}\mathit{ϵ}{)}^{-\mathrm{1}}{h}^{-\mathrm{1}/n}\right|q{|}^{\mathrm{1}/n+\mathrm{1}}w(\mathit{\chi }{)}^{-\mathrm{2}/n-\mathrm{2}}\\ {\displaystyle }& {\displaystyle }+\mathit{\gamma }\mathit{\theta }{h}^{-m-\mathrm{1}}\left|q{|}^{m+\mathrm{1}}w\right(\mathit{\chi }{)}^{-m}\\ {\displaystyle }& {\displaystyle }+(\mathrm{4}\mathit{ϵ}{)}^{-\mathrm{1}}\left[\mathrm{1},-,(,\mathrm{1},-,\mathit{\theta },),r\right]\left[-,h^{\mathrm{2}},w,\left(,\mathit{\chi },\right),|,\mathit{\sigma },{|}^{n-\mathrm{1}},\mathit{\sigma }\right.\\ {\displaystyle }& {\displaystyle }\left..+,h,w,\left(,\mathit{\chi },\right),a,\left(,\mathit{\chi },\right),-,h,q,w^{\prime },\left(,\mathit{\chi },\right),w,(,\mathit{\chi },{)}^{-\mathrm{1}}\right]\\ {\displaystyle }& {\displaystyle }-\left(\mathrm{4}\mathit{ϵ}{)}^{-\mathrm{1}}\mathit{\theta }q{b}^{\prime }\right(\mathit{\chi })+h|\mathit{\sigma }{|}^{n+\mathrm{1}}w\left(\mathit{\chi }\right)-\mathit{\sigma }w\left(\mathit{\chi }\right)a\left(\mathit{\chi }\right)\\ & {\displaystyle }& {\displaystyle }+\mathit{\sigma }{w}^{\prime }\left(\mathit{\chi }\right)w(\mathit{\chi }{)}^{-\mathrm{1}},& {\displaystyle }& {\displaystyle }{q}_{\mathit{\eta }}=qaw\left(\mathit{\chi }\right),& {\displaystyle }& {\displaystyle }{\mathit{\chi }}_{\mathit{\eta }}=q,\end{array}$$ with $\mathit{\theta }(h,\mathit{\chi })=\mathrm{1}$ if $h>-b\left(\mathit{\chi }\right)/r$ and $\mathit{\theta }=\mathrm{0}$ otherwise; here $a$, $w$ and $b$ are treated as prescribed functions, and the prime simply indicates their first derivative.

We assume there is an ice divide at $x=\mathrm{0}$, where $u=q=\mathrm{0}$. Technically, the ice divide then becomes a fixed point of the system Eq. (6) approached as $\mathit{\eta }\to -\mathrm{\infty }$, at which $(h,\mathit{\sigma },q,\mathit{\chi })=(h_{\mathrm{0}},[a\left(\mathrm{0}\right)/h_{\mathrm{0}}{]}^{\mathrm{1}/n},\mathrm{0},\mathrm{0})$ with the ice divide thickness $h_{\mathrm{0}}>\mathrm{0}$ not known a priori. The trick is to determine the value of $h_{\mathrm{0}}$ for which the boundary conditions at the glacier terminus are satisfied by means of a shooting method. Given $h_{\mathrm{0}}$, the fixed point has a unique orbit that emerges from it. In other words, if $h_{\mathrm{0}}$ is known, then the solution to Eq. (6) can be computed uniquely. A constraint on $h_{\mathrm{0}}$ therefore arises from imposing the boundary conditions at the glacier terminus, which are of the form $$\begin{array}{cc}{\displaystyle }& {\displaystyle }\mathrm{4}\mathit{\epsilon }h\mathit{\sigma }=(\mathrm{1}-r)h^{\mathrm{2}}/\mathrm{2}+\mathit{\theta }\left(r^{\mathrm{2}},h^{\mathrm{2}},-,b,(,\mathit{\chi },{)}^{\mathrm{2}}\right)/\left(\mathrm{2}r\right),\\ & {\displaystyle }& {\displaystyle }h=-b\mathit{\varphi }(-\mathit{\lambda }b(\mathit{\chi }{)}^{-\mathrm{1}})\end{array}$$ at some finite $\mathit{\eta }={\mathit{\eta }}_{\mathrm{c}}$. Equation () is dealt with simply by integrating along the orbit until the first condition is satisfied. The second then acts as a single constraint on the degree of freedom $h_{\mathrm{0}}$ that uniquely determines the solution. The code used to compute solutions here is included in the Supplement.

Steady-state profiles, with $w\equiv \mathrm{1}$, $n=\mathrm{1}/m=\mathrm{3}$, $r=\mathrm{0.9}$, $\mathit{\gamma }=\mathrm{0}$, $\mathit{\epsilon }=\mathrm{0.01}$, $\mathit{\lambda }=\mathrm{7.0304}r{\mathit{\epsilon }}^{n^{\mathrm{2}}/(n+\mathrm{1}{)}^{\mathrm{2}}}$, $b\left(x\right)={\mathit{\epsilon }}^{n^{\mathrm{2}}/(n+\mathrm{1})}(\mathrm{1}-\mathrm{2}{x}^{\mathrm{2}})$. The apparently contrived form in which $\mathit{\lambda }$ and $b$ are written is designed to make comparison with Sect.  simpler, whereby the corresponding variables $\mathrm{\Lambda }$ and $B$ are then easy to extract. Panel (e) shows grounding line positions $x_{\mathrm{g}}$ against surface accumulation rate $a$ for the CD calving model (solid line, grounded where the background is shaded grey) and FL calving law (dashed line). (a–d) Steady-state profiles, same colour scheme as in Fig. , same horizontal axis as in (e). (a) The FL calving law, $a=\mathrm{2.08}$. (b–d) CD calving law with $a=$ 0.134 (b), 4.18 (c), 3.54 (d). Red diamonds in (e) refer to the steady state in the panel indicated by the letter label.

[Figure omitted. See PDF]

Steady-state profiles, with $w\equiv \mathrm{1}$, $n=\mathrm{1}/m=\mathrm{3}$, $r=\mathrm{0.9}$, $\mathit{\gamma }=\mathrm{0}$, $\mathit{\epsilon }=\mathrm{0.001}$, $b\left(x\right)={\mathit{\epsilon }}^{n^{\mathrm{2}}/(n+\mathrm{1})}(\mathrm{9.4665}-\mathrm{28.3710}{x}^{\mathrm{2}}+\mathrm{13.3975}{x}^{\mathrm{4}}-\mathrm{1.97021}{x}^{\mathrm{6}})$. Panel (b) shows the geometry of the bed, which is a scaled version of the sixth-order polynomial bed shape used in . Panel (a) shows grounding line positions $x_{\mathrm{g}}$ against surface accumulation rate $a$ (note the logarithmic scale) for the FL model (dashed line) and the CD calving model (solid line, shown in black where the calving front is grounded, and blue where it is afloat) at values of $\mathit{\lambda }=\mathrm{4.7486}r{\mathit{\epsilon }}^{n^{\mathrm{2}}/(n+\mathrm{1}{)}^{\mathrm{2}}}$, $\mathrm{5.5122}r{\mathit{\epsilon }}^{n^{\mathrm{2}}/(n+\mathrm{1}{)}^{\mathrm{2}}}$, $\mathrm{6.2758}r{\mathit{\epsilon }}^{n^{\mathrm{2}}/(n+\mathrm{1}{)}^{\mathrm{2}}}$, $\mathrm{7.0304}r{\mathit{\epsilon }}^{n^{\mathrm{2}}/(n+\mathrm{1}{)}^{\mathrm{2}}}$, $\mathrm{7.8030}r{\mathit{\epsilon }}^{n^{\mathrm{2}}/(n+\mathrm{1}{)}^{\mathrm{2}}}$ and $\mathrm{8.5666}r{\mathit{\epsilon }}^{n^{\mathrm{2}}/(n+\mathrm{1}{)}^{\mathrm{2}}}$. To identify the curves, note that, for a given grounding line position, larger values of $\mathit{\lambda }$ invariably correspond to larger values of $a$.

[Figure omitted. See PDF]

Figure  shows a synthetic example, not based on trying to emulate any specific glacier geometry. Parameter values are given in the figure caption: here, we hold constant the bed and channel geometries; material properties such as $m$, $n$ and $r$; and the calving parameter $\mathit{\lambda }$ fixed and vary accumulation rate $a$. Importantly, the bed slopes downward monotonically in $x$. For “unbuttressed” glaciers subject to basal but not lateral drag e.g., this would lead to steady-state grounding line position $x_{\mathrm{g}}$ increasing monotonically with accumulation rate $a$: increased accumulation inland must be balanced by increased discharge of ice across the grounding line, which happens for unbuttressed glaciers when the grounding line moves into deeper water. The most notable feature in panel (b) is that this behaviour persists if we use the FL model (dashed line). For the CD model, we see a partial reversal of this behaviour: for accumulation rates larger than a certain value, there appear to be no steady-state solutions at all. For smaller accumulation rates, there are two steady-state solutions: (i) a large ice sheet for which $x_{\mathrm{g}}$ shrinks as $a$ increases and (ii) a small ice sheet for which $x_{\mathrm{g}}$ increases with $a$. The larger solution branch also contradicts existing understanding of marine ice sheet dynamics, precisely because an increase in surface mass balance causes the grounding line to retreat into shallower water. Such steady states are likely to be unstable (see , and Sect.  below)

Figure  shows analogous calculations to Fig.  but for an overdeepened bed shape based on that used in . For the FL model, we invariably see that an increase in accumulation rate makes the grounding line advance on a downward slope and retreat on a retrograde slope. This is again analogous to the unbuttressed case studied in , where the grounding line is then unstable when located on an upward slope. For the CD model, the behaviour becomes more complicated. We see that the grounding line can either advance or retreat with increasing accumulation rate, on both the downward- and upward-sloping parts of the bed. Qualitatively, shallow water depths at the grounding line are more commonly associated with the standard, “unbuttressed” behaviour (that is, an increase in accumulation tends to cause the grounding line to advance on downward slopes and retreat on retrograde slopes). The reverse behaviour is typically associated with larger water depths at the grounding line. We also see that a decrease in $\mathit{\lambda }$ leads to the “reverse” behaviour being observed down to shallower water depths at the grounding line and, in particular, through more of the overdeepened section. Note also that the solid (CD model) solution curves in Fig.  end at finite values of $a$ at $x_{\mathrm{g}}=\mathrm{0.6386}$, the location where water depth goes to zero: as shown in Fig. b, an oddity of the CD model is that it predicts a non-zero calving front thickness even when the water depth is zero, and hence there is a non-zero calving flux even where the “grounding line” is on dry land; we have only computed solutions where the grounding line remains in the water.

Our aim in what follows is to explain the results in Figs. – using the same boundary layer approach as in . In particular, we will show that flux through the grounding line can be computed to leading order in the parameter $\mathit{\epsilon }$ as a function of depth to bedrock and channel width at the grounding line, as well as of the calving parameter $\mathit{\lambda }$, friction coefficient $\mathit{\gamma }$ and the remaining physical parameters ($r$, $m$, $n$). Given such a relationship, it is then possible to determine how the grounding line location in a steady state depends on accumulation rates, purely by balancing net accumulation over the domain with outflow of ice through the grounding line.

Approximation: a small lateral aspect ratio

A local-force-balance version of the model

If we take $\mathit{\epsilon }$ as defined in Eq. () with ${\stackrel{\mathrm{‾}}{B}}^{\prime }$ given by Eq. (), we have $${\displaystyle }\mathit{\epsilon }=(n+\mathrm{2}{)}^{\mathrm{1}/n}{\mathrm{2}}^{\mathrm{1}/n}{\left({\displaystyle \frac{\left[w\right]}{\left[x\right]}}\right)}^{(n+\mathrm{1})/n}.$$ In other words, $\mathit{\epsilon }$ is a measure of the lateral aspect ratio $\left[w\right]/\left[x\right]$. A narrow channel ensures that $\mathit{\epsilon }$ is small, which is the basis for our approximation scheme. With $\mathit{\epsilon }$ small, we can neglect gradients of the depth-integrated extensional stress in Eq. () (that is, gradients of $\mathrm{4}h|u_x{|}^{\mathrm{1}/n-\mathrm{1}}{u}_x$) everywhere except close to the terminus, and we find $$\begin{array}{}& {\displaystyle }& {\displaystyle }-w^{-\mathrm{1}/n-\mathrm{1}}h|u{|}^{\mathrm{1}/n-\mathrm{1}}u-\mathit{\gamma }|u{|}^{m-\mathrm{1}}u-h(h_x+b_x)=\mathrm{0},& {\displaystyle }& {\displaystyle }w{h}_t+(wuh{)}_x=wa.\end{array}$$ For the case of $m=\mathrm{1}/n$ (which arises naturally from theories of hard-bed sliding; see, e.g., ; ), we get a diffusive model for ice thickness evolution: $$w{h}_t-{\left[{\displaystyle \frac{w^{n+\mathrm{2}}{h}^{n+\mathrm{1}}}{(h+\mathit{\gamma }{w}^{(n+\mathrm{1})/n}{)}^n}},\left|,\right(,h,+,b,{)}_x,{|}^{n-\mathrm{1}},(,h,+,b,{)}_x\right]}_x=wa.$$ This is essentially analogous to “shallow-ice” models in ice sheet flow : we have a local balance of forces and an ice flux that depends on ice thickness and ice surface slope see also. Other choices of $m>\mathrm{0}$ also imply an ice flux $uh$ that is an increasing function of width $w$, thickness $h$ and surface slope $-(h_x+b_x)$, but that flux cannot be computed in closed form.

The grounding line boundary layer

The model Eq. (7) holds everywhere except near the grounding line and in the floating ice shelf. Following , we can use the method of matched asymptotic expansions to capture the behaviour of ice flow in that region. This requires us to rescale the dimensionless model Eq. (4) to bring back extensional stress at leading order while maintaining an $O\left(\mathrm{1}\right)$ ice flux $q=uh$. The appropriate rescaling turns out to be $$\begin{array}{cc}{\displaystyle }& {\displaystyle }X={\mathit{\epsilon }}^{-n/(n+\mathrm{1})}(x-x_{\mathrm{g}}),T=t,H={\mathit{\epsilon }}^{-n^{\mathrm{2}}/(n+\mathrm{1}{)}^{\mathrm{2}}}h,\\ {\displaystyle }& {\displaystyle }U={\mathit{\epsilon }}^{n^{\mathrm{2}}/(n+\mathrm{1}{)}^{\mathrm{2}}}u.\end{array}$$ By contrast with , we also have to include a potentially non-zero ice shelf length here, so we put $${\displaystyle }{X}_{\mathrm{c}}={\mathit{\epsilon }}^{-n/(n+\mathrm{1})}(x_{\mathrm{c}}-x_{\mathrm{g}}).$$ We treat $H$ and $U$ as functions of $(X,T)$ and $X_{\mathrm{c}}$ as a function of $T$. The rescaling in $H$ implies that ice thickness at the grounding line must be small compared with the interior of the glacier. If there is a floating portion, the glacier must however also reach its flotation thickness at the grounding line. We assume that the glacier is at least near flotation if it has a calving front that remains above flotation. This implies that bed elevation must be small compared with ice thickness in the interior. Specifically, we rescale $${\displaystyle }B={\mathit{\epsilon }}^{-n^{\mathrm{2}}/(n+\mathrm{1}{)}^{\mathrm{2}}}b$$ and assume that $B\sim O\left(\mathrm{1}\right)$; the analogous case of laterally unconfined flow discussed in also requires shallow bed topography. In addition, we assume that thickness $b$ and width $w$ change significantly only over length scales associated with the glacier length as a whole. Over the short length scale associated with the boundary layer coordinate $X$, we treat $B={\mathit{\epsilon }}^{-n^{\mathrm{2}}/(n+\mathrm{1}{)}^{\mathrm{2}}}b\left(x_{\mathrm{g}}\right)$ and width $W=w\left(x_{\mathrm{g}}\right)$ as constant. These additional constraints are again analogous to those made in and imply that we should treat $b$ and $b_x$ as small in the outer problem Eq. (7). With $B$ constant at leading order in the boundary layer, we can also define a scaled flotation thickness $${\displaystyle }{H}_{\mathrm{f}}=-r^{-\mathrm{1}}B,$$ which we will use throughout the rest of the paper as a proxy for water depth to bedrock.

In order for the rescaling in $H$ above to be consistent, we also require that the calving front thickness be similarly small. This turns out to require that $\mathit{\lambda }\sim O\left({\mathit{\epsilon }}^{n^{\mathrm{2}}/(n+\mathrm{1}{)}^{\mathrm{2}}}\right)$, and we define the calving parameter $${\displaystyle }\mathrm{\Lambda }=r^{-\mathrm{1}}{\mathit{\epsilon }}^{-n^{\mathrm{2}}/(n+\mathrm{1}{)}^{\mathrm{2}}}\mathit{\lambda }$$ assuming that $\mathrm{\Lambda }=O\left(\mathrm{1}\right)$; all this implies is that water depths in surface crevasses are not so large as to create calving cliff heights much larger than the expected depth to bedrock at the grounding line.

The result is a boundary layer model at leading order in $\mathit{\epsilon }$. We do not give the detailed derivation here but merely state its form: $$\begin{array}{cc}{\displaystyle }& {\displaystyle }\mathrm{4}\left(H\right|U_X{|}^{\mathrm{1}/n-\mathrm{1}}{U}_X{)}_X-W^{-\mathrm{1}/n-\mathrm{1}}H|U{|}^{\mathrm{1}/n-\mathrm{1}}U\\ {\displaystyle }& {\displaystyle }-\mathit{\gamma }{\mathit{\epsilon }}^{n[\mathrm{1}-n(m+\mathrm{1}\left)\right]/(n+\mathrm{1}{)}^{\mathrm{2}}}\mathit{\theta }|U{|}^{m-\mathrm{1}}U\\ & {\displaystyle }& {\displaystyle }-\left[\mathrm{1},-,(,\mathrm{1},-,\mathit{\theta },),r\right]H{H}_X=\mathrm{0},& {\displaystyle }& {\displaystyle }(HU{)}_X=\mathrm{0},\end{array}$$ for $X<X_{\mathrm{c}}$, where $$\mathit{\theta }=\mathrm{1}\text{for }H\ge H_{\mathrm{f}},\mathit{\theta }=\mathrm{0}\text{otherwise.}$$ The additional boundary condition at the calving front takes the form $$\begin{array}{cc}{\displaystyle }& {\displaystyle }\mathrm{4}H|U_X{|}^{\mathrm{1}/n-\mathrm{1}}{U}_X=(\mathrm{1}-(\mathrm{1}-\mathit{\theta })r)H^{\mathrm{2}}/\mathrm{2}-\mathit{\theta }r{H}_{\mathrm{f}}^{\mathrm{2}}/\mathrm{2}\\ & {\displaystyle }& {\displaystyle }\text{at }X=X_{\mathrm{c}},& {\displaystyle }& {\displaystyle }H=r{H}_{\mathrm{f}}\mathit{\varphi }\left(\mathrm{\Lambda },H_{\mathrm{f}}^{-\mathrm{1}}\right)\text{at }X=X_{\mathrm{c}}.\end{array}$$ Note that the boundary layer is in a pseudo-steady state. This is again analogous to ; the timescale for dynamic adjustment of ice thickness in the boundary layer and of calving front position relative to the grounding line is much shorter than the timescale relevant to the evolution problem Eq. (7) see also. We emphasize that there is no assumption here that the glacier as a whole is in steady state.

In order to make the balances in Eq. (9) work with $\mathit{\epsilon }\ll \mathrm{1}$, we have to deal with the remaining coefficient that contains a power of $\mathit{\epsilon }$ in Eq. (9). We now make further assumptions about the physics of the flow near the grounding line. Our fundamental assumption will be that lateral drag $-w^{-\mathrm{1}/n-\mathrm{1}}h|u{|}^{\mathrm{1}/n-\mathrm{1}}u$ plays a leading order role in force balance at the grounding line but that the floating ice shelf, if it exists, is not so long as to fully buttress the grounding line. By this, we mean that the depth-integrated extensional stress $\mathrm{4}h|u_x{|}^{\mathrm{1}/n-\mathrm{1}}{u}_x$ is comparable in magnitude to $h^{\mathrm{2}}$ all the way up to the grounding line, as is the case at the terminus $x_{\mathrm{c}}$ by dint of the boundary condition Eq. (), but that sidewall drag cannot be neglected.

In that physical regime (termed a “distinguished limit”, in which all physical processes are potentially active), we have to assume that the basal drag coefficient in Eq. () and the calving coefficient in Eq. () are both of $O\left(\mathrm{1}\right)$, meaning that the parameter $\mathrm{\Gamma }$ defined through $$\mathrm{\Gamma }={\mathit{\epsilon }}^{n[\mathrm{1}-n(m+\mathrm{1}\left)\right]/(n+\mathrm{1}{)}^{\mathrm{2}}}\mathit{\gamma }$$ is of $O\left(\mathrm{1}\right)$. We confine our analysis to parameter regimes where this is the case. Note that with $m>\mathrm{0}$ and $n\ge \mathrm{1}$ this implies, strictly speaking, that $\mathit{\gamma }\ll \mathrm{1}$, and basal friction upstream of the boundary layer is formally small in the parameter regime we are considering.

Asymptotic matching is the mathematical formalism by which the boundary layer problem and the “outer” problem Eq. (7) for the dynamics of the bulk of the glacier are connected . With the assumptions on $\mathrm{\Gamma }$ in place, this leads to so-called matching conditions between the boundary layer and the outer problem: $$\begin{array}{cc}{\displaystyle }& {\displaystyle }\underset{X\to -\mathrm{\infty }}{lim}UH=Q=\underset{x\to x_{\mathrm{g}}^{-}}{lim}(-w^{n+\mathrm{1}}h|h_x{|}^{n-\mathrm{1}}{h}_x),\\ {\displaystyle }& {\displaystyle }{W}^{-\mathrm{1}/n-\mathrm{1}}Q|U{|}^{\mathrm{1}/n-\mathrm{1}}\sim -(Q/U\left)\right(Q/U{)}_X,\\ & {\displaystyle }& {\displaystyle }U\to \mathrm{0}\text{as}X\to -\mathrm{\infty }.\end{array}$$ Here $Q$ is the flux at the boundary of the domain of the outer problem, to be determined through the solution of the boundary layer problem. Physically, the first condition states that the flux near the grounding line in the outer problem is the flux that enters the boundary layer at its upstream end. The second condition states that, near that upstream end of the boundary layer, extensional stress gradients have become insignificant and flux is given by a shallow-ice type formula (with $U=Q/H$, the condition can be re-written as $Q\sim -W^{n+\mathrm{1}}H|H_X{|}^{n-\mathrm{1}}{H}_X$, the appropriate local-force-balance formula in our case). Lastly, the third condition states that velocities in the interior of the boundary layer are large compared with those in the rest of the glacier. Structurally, the boundary layer problem above is very similar to that in , with additional physics due to lateral shearing and calving accommodated at the cost of a more complicated formulation.

From the perspective of the model Eq. () for the dynamics of the glacier as a whole, the purpose of the boundary layer is to provide the relevant boundary conditions at $x=x_{\mathrm{g}}$. As Eq. () is a diffusion model for $h$, it requires two boundary conditions at any moving boundary. Where a floating portion exists, a condition on $h$ at the grounding line arises straightforwardly from Eq. (); for a grounded calving front, an equivalent condition is provided by Eq. ().

As in previous work , we can show that the second boundary condition takes the form of a flux condition that can be found by solving the boundary layer problem: the problem Eq. (9) has a solution only if $Q$ satisfies a functional relationship with flotation thickness $H_{\mathrm{f}}$, width $W$, friction coefficient $\mathrm{\Gamma }$ and the calving parameter $\mathrm{\Lambda }$. It is important to emphasize again that $H_{\mathrm{f}}$ need not be the ice thickness at what we have termed the “grounding line”. Instead, $H_{\mathrm{f}}$ is the flotation thickness there, determined purely by bedrock depth, and ice thickness equals $H_{\mathrm{f}}$ at the grounding line only if the glacier has a floating shelf or is at the point of forming one. The flotation thickness $H_{\mathrm{f}}=-B/r$ is of course prescribed for any given grounding line position, as is the channel width $W$. We therefore end up with ice flux as a function of grounding line position, basal drag coefficient $\mathrm{\Gamma }$ and the calving parameter $\mathrm{\Lambda }$, itself a proxy for water depth in surface crevasses.

We give additional detail on how to compute that relationship between flux, geometry and model parameters in Appendix  and in the Supplement, and the code used to solve the problem is also included in the Supplement. Importantly, we are able to show the relationship takes the form $$Q=W{H}_{\mathrm{f}}^{n+\mathrm{1}}{G}_{\mathrm{\Lambda }}\left({\displaystyle \frac{H_{\mathrm{f}}}{\mathrm{\Lambda }}},,,{\displaystyle \frac{\mathrm{\Gamma }{W}^{(nm+n+m+\mathrm{1})/(n+\mathrm{1})}}{{\mathrm{\Lambda }}^{\mathrm{2}-nm}}},,,n,,,m,,,r\right).$$

The practical use of this form is that it reduces the complexity of the flux formula: for a given set of constants $n$, $m$ and $r$, what we primarily need to calculate is the dependence of flux on the first two arguments of the function $G_{\mathrm{\Lambda }}$. For the FL law, we can in fact go further and use the fact that flux cannot depend on the now redundant parameter $\mathrm{\Lambda }$ to simplify the expression to $$Q=W{H}_{\mathrm{f}}^{n+\mathrm{1}}{G}_{\mathrm{\Gamma }}\left({\displaystyle \frac{\mathrm{\Gamma }{W}^{(nm+n+m+\mathrm{1})/(n+\mathrm{1})}}{H_{\mathrm{f}}^{\mathrm{2}-nm}}},,,n,,,m,,,r\right).$$

Solutions of the boundary layer problem

Equation () allows us to collapse solutions for flux onto a one-parameter family of plots for each of the two calving laws considered (the FL and CD calving laws). Specifically, we can plot $Q/\left(W{\mathrm{\Lambda }}^{n+\mathrm{1}}\right)$ against $H_{\mathrm{f}}/\mathrm{\Lambda }$ for fixed values of $\mathrm{\Gamma }{W}^{(nm+n+m+\mathrm{1})/(n+\mathrm{1})}/{\mathrm{\Lambda }}^{\mathrm{2}-nm}$. Roughly speaking, we can think of this as plotting flux $Q$ against flotation thickness $H_{\mathrm{f}}$ for different values of the basal drag coefficient. Solutions are plotted in this way in Fig. a. The black curves signify solutions with vanishing basal friction ($\mathrm{\Gamma }=\mathrm{0}$), while coloured curves show solutions with non-zero values of $\mathrm{\Gamma }$ as specified in the figure caption. The dashed line in each case corresponds to the FL model, while the solid line corresponds to the CD calving model.

Solutions of the boundary layer problem. $r=\mathrm{0.9}$; $n=\mathrm{1}/m=\mathrm{3}$. Dashed lines represent the FL law, solid lines the CD calving law. Grey shading indicates a grounded calving front in the CD model, and the white background a floating shelf. (a) Normalized flux against normalized thickness at different basal friction parameter values, $\mathrm{\Gamma }{W}^{(n+\mathrm{1})/n}/\mathrm{\Lambda }=$ 0 (black line), 1 (blue), 5 (red), 25 (green) and 125 (magenta). (b) Logarithm of normalized flux against normalized calving parameter at different basal friction parameter values, $\mathrm{\Gamma }{W}^{(n+\mathrm{1})/n}/H_{\mathrm{f}}=$ 0 (black line), 1 (blue), 5 (red), 25 (green) and 125 (magenta). (c) Normalized flux against normalized channel width at different (grounded) ice thickness: $H_{\mathrm{f}}/\mathrm{\Lambda }=\mathrm{2}$ (black), $\mathrm{5}/\mathrm{3}$ (blue), $\mathrm{4}/\mathrm{3}$ (red), 1 (green) and $\mathrm{2}/\mathrm{3}$ (magenta). (d) Same as panel (c) but with floating ice thickness values $H_{\mathrm{f}}/\mathrm{\Lambda }=\mathrm{2}$ (black), $\mathrm{7}/\mathrm{3}$ (blue), $\mathrm{8}/\mathrm{3}$ (red) and $\mathrm{10}/\mathrm{3}$ (green).

[Figure omitted. See PDF]

As already suggested by the steady-state solutions to the full problem in Fig. , Fig.  confirms that flux is not a monotonically increasing function of flotation thickness $H_{\mathrm{f}}$ for the CD model. We have what we term an anomalous flux–flotation-thickness relationship for large enough values of $H_{\mathrm{f}}$: flux $Q$ actually decreases with increasing flotation thickness $H_{\mathrm{f}}$ for all but relatively small $H_{\mathrm{f}}$, at least for moderate or small basal drag coefficients. For large values of the basal drag coefficient (the $\mathrm{\Gamma }{W}^{(n+\mathrm{1})/n}/\mathrm{\Lambda }=\mathrm{25}$ and $\mathrm{125}$ cases shown), the relationship between $Q$ and $H_{\mathrm{f}}$ is even more complicated. We have the same anomalous flux–flotation-thickness relationship as for small basal drag while the calving front is grounded, but for a floating calving front we find that flux $Q$ increases again with $H_{\mathrm{f}}$ (this is even clearer in of Fig. a, which is a zoomed-in version of Fig. ). In all cases, the flux for the CD calving model approaches the same limit for large $H_{\mathrm{f}}$, independently of the calving law.

By contrast, the flux always increases with flotation thickness in the calving at floatation model, just as it does in laterally unconfined marine ice sheet flow . In fact, Eq. () already told us as much for the case of vanishing basal friction coefficient $\mathrm{\Gamma }$. Note that the flux curve for the CD model and for the FL model always has a point of intersection at $H_{\mathrm{f}}/\mathrm{\Lambda }=\mathrm{2}$. From the definition of $\mathit{\varphi }$ in Eq. (), it is easy to see that this is the point at which the calving front is just at flotation in the CD model. Therefore, the model produces the same result as the FL model. For smaller values of $H_{\mathrm{f}}/\mathrm{\Lambda }$, the CD model has a grounded calving front, while the calving front becomes the terminus of a floating ice shelf at larger values of $H_{\mathrm{f}}/\mathrm{\Lambda }$. Note that flux in the CD model is always a decreasing function of $H_{\mathrm{f}}/\mathrm{\Lambda }$ for $H_{\mathrm{f}}/\mathrm{\Lambda }$ slightly less than the critical value of 2 for changeover from a grounded to a floating terminus. This observation will be key to our interpretation of the physics involved in the anomalous flux–thickness relationship.

Other features of our solutions are also shown in panels (a)–(c) of Fig. . Each panel isolates one parameter ($H_{\mathrm{f}}$, $\mathrm{\Lambda }$, $W$) on the horizontal axis but normalizes it as dictated by Eq. (10) and plots it against an also normalized flux (again as dictated by Eq. 10) on the vertical axis. Apart from the dependence of $Q$ on $H_{\mathrm{f}}$, panel (a) also shows that flux always decreases with increasing friction coefficient $\mathrm{\Gamma }$, while panel (c) shows that flux increases with channel width $W$. This holds regardless of the calving model used and is what one expects intuitively: wider channels and lower basal drag ought to speed up ice flow and lead to larger ice discharge. Panel (b) shows that, for the CD calving law, flux also increases with the calving parameter $\mathrm{\Lambda }$: recall that the calving parameter $\mathrm{\Lambda }$ is a dimensionless version of water depth in surface crevasses, and larger values of $\mathrm{\Lambda }$ lead to taller calving cliffs and hence to larger extensional stresses near the grounding line. In fact, the flux is far more sensitive to changes in $\mathrm{\Lambda }$ than in any other parameter: notice that panel (b) plots the logarithm of the flux on the vertical axis.

We can also confirm our boundary layer results by direct comparison with numerical solutions of the full ice flow problem, computed by the method in Sect. . This is shown for the case of vanishing basal friction in Fig. . Here we use the same parameter values as in Fig.  but for two different values of $\mathit{\epsilon }$. Different solutions to the steady-state problem are again obtained by varying $a$. For each $a$, we plot ice flux across the grounding line in the steady-state solution, scaled as in Sect. , against flotation thickness at the grounding line, also scaled as in Sect. . As expected, the flux solutions obtained from the full steady-state problem Eq. (6) for the CD and FL models converge to those obtained from the boundary layer problem as $\mathit{\epsilon }$ is made smaller: for $\mathit{\epsilon }={\mathrm{10}}^{-\mathrm{4}}$, the flux curves are virtually indistinguishable, confirming the accuracy of the boundary layer solution.

There are two aspects of the CD model flux solution that we still need to explain in more detail: (i) why flux decreases with increasing flotation thickness $H_{\mathrm{f}}$ in some circumstances and (ii) why flux approaches a constant limit for large $H_{\mathrm{f}}$ and so becomes independent of depth to bedrock in the channel, depending instead only on the calving parameter $\mathrm{\Lambda }$. We turn to these problems next.

Comparison of the solution of the boundary layer problem with $\mathrm{\Gamma }=\mathrm{0}$ (black lines, solid for the CD model, dashed for FL) and solutions of the steady-state problem Eq. (6), solved with the parameter values given in the caption of Fig. . Plotted is scaled, normalized ice flux against scaled, normalized ice thickness for $\mathit{\epsilon }={\mathrm{10}}^{-\mathrm{2}}$ (blue symbols) and $\mathit{\epsilon }={\mathrm{10}}^{-\mathrm{4}}$ (red symbols). Diamonds show solutions for the CD model, and circles show solutions for the calving model.

[Figure omitted. See PDF]

The role of extensional stress at the grounding line

Key to the flux–flotation-thickness relationship is that flux depends on the extensional stress at the grounding line, and that extensional stress in turn depends on the geometry of the calving front and floating ice shelf. For relatively small extensional stresses $\mathrm{\Sigma }$ defined by $${\displaystyle }\mathrm{\Sigma }=\mathrm{4}|U_X{|}^{\mathrm{1}/n-\mathrm{1}}{U}_X$$ (“small” meaning that $\mathrm{\Sigma }$ is much smaller than $H$), it is possible to derive an approximate formula for flux in terms of ice thickness $H_{\mathrm{g}}=H\left(\mathrm{0}\right)$ and extensional stress ${\mathrm{\Sigma }}_{\mathrm{g}}=\mathrm{\Sigma }\left(\mathrm{0}\right)$ at the grounding line $X=\mathrm{0}$. For the remainder of this section, we will use the commonly assumed friction exponent $m=\mathrm{1}/n$. We obtain from Eq. () that $${\displaystyle }-H_X=H{U}_X/U,$$ and if we neglect gradients of $H\mathrm{\Sigma }$ in force balance, then Eq. () leads to $${\displaystyle }|U{|}^{(n+\mathrm{1})/n}\approx H^{\mathrm{2}}{U}_X/(\mathrm{\Gamma }+W^{-(n+\mathrm{1})/n}H).$$ Therefore, with the condition $\mathrm{4}|U_X{|}^{\mathrm{1}/n-\mathrm{1}}{U}_X={\mathrm{\Sigma }}_{\mathrm{g}}$, $H=H_{\mathrm{g}}$ at $X=\mathrm{0}$, $$Q\approx {\left({\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\right)}^{n^{\mathrm{2}}/(n+\mathrm{1})}{\displaystyle \frac{H_{\mathrm{g}}^{(\mathrm{3}n+\mathrm{1})/(n+\mathrm{1})}{\mathrm{\Sigma }}_{\mathrm{g}}^{n^{\mathrm{2}}/(n+\mathrm{1})}}{(\mathrm{\Gamma }+W^{-(n+\mathrm{1})/n}{H}_{\mathrm{g}}{)}^{n/(n+\mathrm{1})}}}.$$ This formula is essentially a modification of formula (29) in , and its derivation is a translation of Appendix A of to our modified boundary layer problem. The omission of the extensional stress gradient can also be formalized on the basis that the density difference $(\mathrm{1}-r)$ between ice and water is small, leading to gradients of $H\mathrm{\Sigma }$ being negligible in the balance of forces (see the Supplement).

For the FL model, it is easy to extract an analytical formula for flux as a function of channel width and depth to bedrock from Eq. (). Specifically, we have $H_{\mathrm{g}}=H_{\mathrm{f}}$ and ${\mathrm{\Sigma }}_{\mathrm{g}}=(\mathrm{1}-r)H_{\mathrm{f}}/\mathrm{2}$, so we get $$Q\approx {\left({\displaystyle \frac{\mathrm{1}-r}{\mathrm{8}}}\right)}^{n^{\mathrm{2}}/(n+\mathrm{1})}{\displaystyle \frac{H_{\mathrm{f}}^{(n^{\mathrm{2}}+\mathrm{3}n+\mathrm{1})/(n+\mathrm{1})}}{{\left(\mathrm{\Gamma },+,W^{-(n+\mathrm{1})/n},H_{\mathrm{f}}\right)}^{n/(n+\mathrm{1})}}},$$ which simply generalizes formula (3.51) in with $m=\mathrm{1}/n$ to the case of lateral as well as basal drag. With $\mathrm{\Gamma }=\mathrm{0}$, we can also immediately recognize a version of formula () with $$Q\approx {\left({\displaystyle \frac{\mathrm{1}-r}{\mathrm{8}}}\right)}^{n^{\mathrm{2}}/(n+\mathrm{1})}W{H}_{\mathrm{f}}^{n+\mathrm{1}}.$$

Figure b and c show that Eq. () performs well for $\mathrm{\Gamma }=\mathrm{0}$ and $\mathrm{\Gamma }=\mathrm{125}$. Clearly, Eq. () predicts flux increasing with flotation thickness; this is the result of both grounding line thickness $H_{\mathrm{g}}$ and extensional stress ${\mathrm{\Sigma }}_{\mathrm{g}}$ increasing with flotation thickness $H_{\mathrm{f}}$ at the grounding line. Next, we will use Eq. () to explain the anomalous behaviour with the CD calving law.

Limiting forms of ice flux near the critical ice thickness $H_{\mathrm{f}}/\mathrm{\Lambda }=\mathrm{2}$ and for large $H_{\mathrm{f}}/\mathrm{\Lambda }$. (a) A zoomed-in version of Fig. a, same colour scheme. The black, red and blue curves exhibit the anomalous relationship between flux and flotation thickness for a floating calving front, the green and magenta curves (with larger basal friction coefficients) do not. The “dominant stress” labels refer to terms in the force balance of the shelf that balance the depth-integrated extensional stress gradient $(H\mathrm{\Sigma }{)}_X$ for different basal drag coefficients $\mathrm{\Gamma }$; see Sect. . Red markers correspond to profiles shown in Fig. . (b) Solution to the boundary layer problem for $\mathrm{\Gamma }{W}^{(n+\mathrm{1})/n}/\mathrm{\Lambda }=\mathrm{0}$ plotted in black (same as in panel a). Formula () is shown as a blue dashed line, formulae () and () as blue solid lines, while the dot-dashed blue line indicates the long shelf limit Eq. (). Red markers correspond to the profiles shown in Fig. a and b. (c) Same as panel (b) but with the boundary layer solution for $\mathrm{\Gamma }{W}^{(n+\mathrm{1})/n}/\mathrm{\Lambda }=\mathrm{125}$ shown in magenta and with red markers corresponding to profiles in Fig. c.

[Figure omitted. See PDF]

Grounded calving fronts

Here, we are interested in the anomalous relationship between $Q$ and the flotation thickness $H_{\mathrm{f}}$. Recall that $H_{\mathrm{f}}=-B/r$ is given by depth to bedrock $B$ and is therefore prescribed for a given grounding line location. The actual ice thickness $H_{\mathrm{g}}$ at the grounding line is equal to $H_{\mathrm{f}}$ when the glacier has a floating ice shelf or calves at flotation; for a grounded calving cliff, $H_{\mathrm{f}}$ may exceed $H_{\mathrm{g}}$. However, as Fig.  shows, $H_{\mathrm{g}}$ always increases with $H_{\mathrm{f}}$. Equation () further shows that flux $Q$ increases with ice thickness $H_{\mathrm{g}}$ and extensional stress ${\mathrm{\Sigma }}_{\mathrm{g}}$ at the grounding line. The anomalous relationship must therefore hinge on ${\mathrm{\Sigma }}_{\mathrm{g}}$ decreasing sufficiently rapidly as flotation thickness $H_{\mathrm{f}}$ increases.

Note that the anomalous decrease in flux with increasing flotation thickness is most pronounced around the critical value $H_{\mathrm{f}}/\mathrm{\Lambda }=\mathrm{2}$, where the calving front goes from grounded to floating. We can understand the behaviour of ice flux near that value by considering the effect of small perturbations in $H_{\mathrm{f}}$ away from that critical value. Again, recall the actual thickness at the calving front is given by $H_{\mathrm{c}}=H\left(X_{\mathrm{c}}\right)=r{H}_{\mathrm{f}}\mathit{\varphi }\left(\mathrm{\Lambda }{H}_{\mathrm{f}}^{-\mathrm{1}}\right)$. Let the critical value of $H_{\mathrm{f}}$ be $H_{\mathrm{f}\mathrm{0}}=\mathrm{2}\mathrm{\Lambda }$, for which $H_{\mathrm{c}}=H_{\mathrm{f}\mathrm{0}}$. Now consider perturbing $H_{\mathrm{f}}$ slightly, say to $H_{\mathrm{f}\mathrm{0}}+H_{\mathrm{f}}^{\prime }$. We can use a first-order Taylor expansion of $\mathit{\varphi }$ to compute the perturbed calving front thickness as $$\begin{array}{cc}{\displaystyle }{H}_{\mathrm{c}}& {\displaystyle }=r(H_{\mathrm{f}\mathrm{0}}+H_{\mathrm{f}}^{\prime })\mathit{\varphi }\left({\displaystyle \frac{\mathrm{\Lambda }}{H_{\mathrm{f}\mathrm{0}}+H_{\mathrm{f}}^{\prime }}}\right)\\ {\displaystyle }& {\displaystyle }=r(H_{\mathrm{f}\mathrm{0}}+H_{\mathrm{f}}^{\prime })\left[\mathit{\varphi },\left({\displaystyle \frac{\mathrm{\Lambda }}{H_{\mathrm{f}\mathrm{0}}}}\right),-,{\mathit{\varphi }}^{\prime },\left({\displaystyle \frac{\mathrm{\Lambda }}{H_{\mathrm{f}\mathrm{0}}}}\right),{\displaystyle \frac{\mathrm{\Lambda }{H}_{\mathrm{f}}^{\prime }}{H_{\mathrm{f}\mathrm{0}}^{\mathrm{2}}}},+,O,\left(,{H_{\mathrm{f}}^{\prime }}^{\mathrm{2}},\right)\right],\end{array}$$ where the prime on $\mathit{\varphi }$ denotes an ordinary derivative. For the CD model, Eq. () shows that $\mathit{\varphi }(\mathrm{\Lambda }/H_{\mathrm{f}\mathrm{0}})=\mathit{\varphi }(\mathrm{1}/\mathrm{2})=r^{-\mathrm{1}}$ and ${\mathit{\varphi }}^{\prime }(\mathrm{\Lambda }/H_{\mathrm{f}\mathrm{0}})=\mathrm{2}{r}^{-\mathrm{1}}$, so $$H_{\mathrm{c}}=H_{\mathrm{f}\mathrm{0}}+O\left({H_{\mathrm{f}}^{\prime }}^{\mathrm{2}}\right).$$ In other words, at linear order, a small perturbation in bedrock depth has no effect on ice thickness at the calving front in the CD model. Note that because $\mathit{\varphi }$ in the CD model is continuously differentiable, this holds regardless of whether $H_{\mathrm{f}}^{\prime }$ is positive or negative, that is, regardless of whether the perturbation causes a calving cliff thicker than flotation or a floating ice shelf to form.

Maintaining constant ice thickness at the calving front while bedrock depth changes has a significant effect on the extensional stress at the grounding line. Consider the case of a grounded calving front when $H_{\mathrm{f}}^{\prime }<\mathrm{0}$; the grounding line thickness is the calving front thickness $H_{\mathrm{g}}=H_{\mathrm{c}}\approx H_{\mathrm{f}\mathrm{0}}$. The stress condition Eq. () can also be approximated to first order in $H_{\mathrm{f}}^{\prime }$ as $${\mathrm{\Sigma }}_{\mathrm{g}}=\mathrm{4}|U_X{|}^{\mathrm{1}/n-\mathrm{1}}{U}_X\approx (\mathrm{1}-r)H_{\mathrm{f}\mathrm{0}}/\mathrm{2}-H_{\mathrm{f}}^{\prime }\text{at }X=\mathrm{0}.$$ As we have assumed that the flotation thickness perturbation $H_{\mathrm{f}}^{\prime }$ is negative, we see an increase in extensional stress relative to the unperturbed flotation thickness $H_{\mathrm{f}\mathrm{0}}$. Even for small $H_{\mathrm{f}}^{\prime }$, the increase can be very significant: with $\mathrm{1}-r=\mathrm{0.1}$, $H_{\mathrm{f}}^{\prime }$ only needs to be $\mathrm{1}/\mathrm{20}$ of the unperturbed flotation thickness $H_{\mathrm{f}\mathrm{0}}$ in order for the stress perturbation to be the same size as the unperturbed stress $(\mathrm{1}-r)H_{\mathrm{f}\mathrm{0}}/\mathrm{2}$. An illustration of this effect is given in Fig. a, where a slight decrease in flotation thickness (panel a1) clearly leads to a substantial increase in extensional stress (panel a2).

Boundary layer solutions $r=\mathrm{0.9}$, $n=\mathrm{1}/m=\mathrm{3}$, $\mathrm{\Lambda }=W=\mathrm{1}$. (a1, b1, c1) Boundary layer ice geometry, same colour scheme as in Fig. . (a2, b2, c2) The corresponding extensional stress profiles. Dashed lines correspond to the FL model with $H_{\mathrm{f}}=\mathrm{2}$ (the “unperturbed” flotation thickness $H_{\mathrm{f}\mathrm{0}}$ in Sect.  and ), and solid profiles to perturbed flotation thicknesses. (a, b) $\mathrm{\Gamma }=\mathrm{0}$ (the dashed profiles are identical), with the solid line showing (a) $H_{\mathrm{f}}=\mathrm{1.9}$ (grounded cliff) and (b) $H_{\mathrm{f}}=\mathrm{2.1}$ (floating calving front). (c) $\mathrm{\Gamma }=\mathrm{125}$, with the solid line showing $H_{\mathrm{f}}=\mathrm{2.1}$.

[Figure omitted. See PDF]

The extensional stress perturbation occurs because the calving cliff ice thickness has not changed at first order, but bedrock depth is shallower. The calving cliff now protrudes further above the water line, and the depth-averaged normal stress exerted on it by the water is smaller. As a result, the extensional stress in the ice has to increase. This increase in stress is what leads to the increase in flux caused by the decrease in flotation thickness $H_{\mathrm{f}}$. In fact, for small $H_{\mathrm{f}}^{\prime }<\mathrm{0}$, Eq. () then becomes $$\begin{array}{cc}{\displaystyle }& {\displaystyle }Q\approx {\left({\displaystyle \frac{\mathrm{1}-r}{\mathrm{8}}}\right)}^{n^{\mathrm{2}}/(n+\mathrm{1})}\\ & {\displaystyle }& {\displaystyle }{\displaystyle \frac{H_{\mathrm{f}\mathrm{0}}^{(\mathrm{3}n+\mathrm{1})/(n+\mathrm{1})}{\left[H_{\mathrm{f}\mathrm{0}},-,\mathrm{2},(,\mathrm{1},-,r,{)}^{-\mathrm{1}},H_{\mathrm{f}}^{\prime }\right]}^{n^{\mathrm{2}}/(n+\mathrm{1})}}{{\left(\mathrm{\Gamma },+,W^{-(n+\mathrm{1})/n},H_{\mathrm{f}\mathrm{0}}\right)}^{n/(n+\mathrm{1})}}},\end{array}$$ and $Q$ increases as $H_{\mathrm{f}}^{\prime }$ becomes more negative.

This is consistent with the behaviour shown in Fig. . For grounded calving, we always find the anomalous relationship between $Q$ and $H_{\mathrm{f}}$, regardless of the basal friction parameter. Figure b and c also show that Eq. () is accurate only for very small $H_{\mathrm{f}}^{\prime }$; this is presumably a result of the fact that $\mathrm{1}-r=\mathrm{0.1}$ is not extremely small and of the fact that the quadratic term in Eq. () starts to become large enough to affect results (indeed, Fig. a indicates that a linearization of $\mathit{\varphi }$ is unlikely to be accurate for grounded calving fronts except very close to $H_{\mathrm{f}}/\mathrm{\Lambda }=\mathrm{2}$).

Floating calving fronts

We can conversely take the case of $H_{\mathrm{f}}^{\prime }>\mathrm{0}$, which leads to the formation of a floating ice shelf. As the calving front thickness does not change to first order, the extensional stress at the calving front remains equal to $$\mathrm{\Sigma }\left(X_{\mathrm{c}}\right)=\mathrm{4}|U_X{|}^{\mathrm{1}/n-\mathrm{1}}{U}_X=(\mathrm{1}-r)H_{\mathrm{f}\mathrm{0}}/\mathrm{2}.$$ Suppose that basal drag is not so large as to render lateral drag insignificant on the grounded portion of the boundary layer. In that case, even though the ice at the grounding line is slightly thicker than at the calving front (by $H_{\mathrm{f}}^{\prime }$), the driving stress in the floating ice shelf is small compared with the other forces acting on the shelf. In particular, the surface slope of the ice shelf is small because of the small density difference between ice and water. Most of the reduction in thickness between grounding line and calving front is accounted for by the bottom of the ice shelf sloping upwards (see Fig. b1). The surface slope of the ice shelf (which causes the driving stress) is only $(\mathrm{1}-r)/r$ times the bottom slope. Since the driving stress is weak, the dominant balance of forces on the ice shelf is then between the gradient of depth-integrated extensional stress $(H\mathrm{\Sigma }{)}_X=\mathrm{4}(H|U_X{|}^{\mathrm{1}/n-\mathrm{1}}{U}_X{)}_X$ and the lateral drag $W^{-(n+\mathrm{1})/n}H|U{|}^{\mathrm{1}/n-\mathrm{1}}U$: in other words, Eq. () becomes approximately $$(H\mathrm{\Sigma }{)}_X-W^{-(n+\mathrm{1})/n}H|U{|}^{\mathrm{1}/n-\mathrm{1}}U\approx \mathrm{0},$$ with the driving stress an $O(\mathrm{1}-r)$ correction. It follows that the floating ice shelf acts to reduce extensional stress at the grounding line relative to its value at the calving front; this is the buttressing effect of the ice shelf.

For small $H_{\mathrm{f}}^{\prime }$, we obtain a short ice shelf and can treat $H\approx H_{\mathrm{f}\mathrm{0}}$ and $U\approx Q/H_{\mathrm{f}\mathrm{0}}$ in the shelf as constant, so that $${\displaystyle }\mathrm{\Sigma }\left(X\right)\approx (\mathrm{1}-r)H_{\mathrm{f}\mathrm{0}}/\mathrm{2}-W^{-(n+\mathrm{1})/n}{Q}^{\mathrm{1}/n}{H}_{\mathrm{f}\mathrm{0}}^{-\mathrm{1}/n}(X_{\mathrm{c}}-X),$$ with $X_{\mathrm{c}}$ being the length of the floating ice shelf. This effect – the linear reduction in extensional stress in the floating shelf with distance from the calving front – is illustrated in Fig. b.

The shelf length is dictated by $H_{\mathrm{f}}^{\prime }$. A larger flotation thickness requires a longer ice shelf before the calving front thickness $H_{\mathrm{f}\mathrm{0}}$ is reached, potentially leading to more buttressing. The ice shelf thickness gradient is $-H_X=H{U}_X/U\approx H_{\mathrm{f}\mathrm{0}}^{\mathrm{2}}{\mathrm{\Sigma }}^n/\left({\mathrm{4}}^{n}Q\right)$, and the shelf length $X_{\mathrm{c}}$ is constrained by the fact that the decrease in ice thickness between grounding line and calving front is $H_{\mathrm{f}}^{\prime }\approx {\int }_{\mathrm{0}}^{X_{\mathrm{c}}}-H_X\mathrm{d}X$. This allows us to compute $X_{\mathrm{c}}$ and hence the stress at the grounding line: $$\begin{array}{cc}{\displaystyle }& {\displaystyle }{\mathrm{\Sigma }}_{\mathrm{g}}=\mathrm{\Sigma }\left(\mathrm{0}\right)\approx \left\{{\left[{\displaystyle \frac{(\mathrm{1}-r)H_{\mathrm{f}\mathrm{0}}}{\mathrm{2}}}\right]}^{n+\mathrm{1}}\right.\\ & {\displaystyle }& {\displaystyle }{\left..-,{\displaystyle \frac{(n+\mathrm{1}){\mathrm{4}}^n{Q}^{(n+\mathrm{1})/n}{H}_{\mathrm{f}}^{\prime }}{W^{(n+\mathrm{1})/n}{H}_{\mathrm{f}\mathrm{0}}^{(\mathrm{2}n+\mathrm{1})/n}}}\right\}}^{\mathrm{1}/(n+\mathrm{1})}.\end{array}$$ Substituting Eq. (19) in Eq. (), we find that flux satisfies, for small $H_{\mathrm{f}}^{\prime }>\mathrm{0}$, $$\begin{array}{cc}{\displaystyle }& {\displaystyle }Q\approx {\displaystyle \frac{H_{\mathrm{f}\mathrm{0}}^{(\mathrm{3}n+\mathrm{1})/(n+\mathrm{1})}}{{\left(\mathrm{\Gamma },+,W^{-(n+\mathrm{1})/n},H_{\mathrm{f}\mathrm{0}}\right)}^{n/(n+\mathrm{1})}}}\left\{{\left[{\displaystyle \frac{(\mathrm{1}-r)H_{\mathrm{f}\mathrm{0}}}{\mathrm{8}}}\right]}^{n+\mathrm{1}}\right.\\ & {\displaystyle }& {\displaystyle }{\left..-,{\displaystyle \frac{(n+\mathrm{1})Q^{(n+\mathrm{1})/n}{H}_{\mathrm{f}}^{\prime }}{W^{(n+\mathrm{1})/n}{H}_{\mathrm{f}\mathrm{0}}^{(\mathrm{2}n+\mathrm{1})/n}}}\right\}}^{n^{\mathrm{2}}/(n+\mathrm{1}{)}^{\mathrm{2}}}.\end{array}$$ At first glance, it does not seem that Eq. () is of much use – it defines $Q$ implicitly. However, from Eq. (), it is not difficult to show that an increase in $H_{\mathrm{f}}^{\prime }$ leads to a decrease in extensional stress ${\mathrm{\Sigma }}_{\mathrm{g}}$ at the grounding line and, therefore, to a decrease in flux $Q$. The stress decreases because the ice shelf lengthens as $H_{\mathrm{f}}^{\prime }$ increases and the total amount of lateral drag on the ice shelf increases.

Again, we have given an ad hoc derivation for Eq. (). We can formalize that derivation as shown in the Supplement, once more based on the small density difference $\mathrm{1}-r$. Figure b shows that Eq. () is more qualitatively than quantitatively accurate for the case of no basal friction $\mathrm{\Gamma }=\mathrm{0}$. Again, this is presumably the result of $\mathrm{1}-r$ not being extremely small and of higher-order terms in the approximation scheme above becoming important.

However, as Fig. a also shows, the anomalous behaviour disappears entirely for floating ice shelves when the basal friction coefficient $\mathrm{\Gamma }$ becomes large. In such cases, the argument above must become qualitatively incorrect. The changeover from the anomalous behaviour for flow dominated by lateral drag to the “normal” behaviour obtained with significant basal drag occurs because, when the basal friction coefficient is large, ice velocities near the grounding line become small. This has two effects: it (i) reduces the lateral drag term $W^{-(n+\mathrm{1})/n}H|U{|}^{\mathrm{1}/n-\mathrm{1}}U$ in Eq. () and (ii) increases the thickness gradient and, therefore, the driving stress. Specifically, conservation of mass in the floating shelf dictates that $$H_X=-H{U}_X/U=-H|\mathrm{\Sigma }{|}^{n-\mathrm{1}}\mathrm{\Sigma }/({\mathrm{4}}^{n}U),$$ so that $H_X$ becomes large when $U$ is small. The driving stress $-(\mathrm{1}-r)H{H}_X$ can then no longer be ignored in Eq. (); Eq. () is no longer applicable, and neither is Eq. (). An increase in flotation thickness can now potentially cause an increase in ice flux, at least when the calving front is afloat. This is shown in Fig. c and is also described in more formal detail in the Supplement.

It is relatively straightforward to estimate how large $\mathrm{\Gamma }$ needs to be in order for driving stress to appear at leading order in the shelf. Note that thickness, velocity and stress are continuous across the grounding line. As a result of Eq. (), so is the thickness gradient $H_X$, but not the surface slope. With large $\mathrm{\Gamma }$ on the grounded side, driving stress balances basal shear, $-H{H}_X\sim \mathrm{\Gamma }|U{|}^{\mathrm{1}/n-\mathrm{1}}U$. In order for driving stress to appear at leading order in the shelf, it should be comparable to lateral drag, so $-(\mathrm{1}-r)H{H}_X\sim W^{-(n+\mathrm{1})/n}H|U{|}^{\mathrm{1}/n-\mathrm{1}}U$. It follows that $$\mathrm{\Gamma }{W}^{(n+\mathrm{1})/n}/H\sim \mathrm{\Gamma }{W}^{(n+\mathrm{1})/n}/\left(\mathrm{2}\mathrm{\Lambda }\right)\sim (\mathrm{1}-r{)}^{-\mathrm{1}}.$$ With $r=\mathrm{0.9}$, this corresponds to $\mathrm{\Gamma }{W}^{(n+\mathrm{1})/n}/\mathrm{\Lambda }\sim \mathrm{20}$, consistent with panel (a) in Fig. , where the green line corresponds to $\mathrm{\Gamma }{W}^{(n+\mathrm{1})/n}/\mathrm{\Lambda }=\mathrm{25}$.

Finally, consider the limiting case of very large basal friction coefficient (i.e. $\mathrm{\Gamma }{W}^{(n+\mathrm{1})/n}/H\gg (\mathrm{1}-r{)}^{-\mathrm{1}}$) combined with an ice shelf that has limited extent. By an extension of the argument above, this corresponds to a large driving stress and to lateral drag playing an insignificant role in force balance. In this case, we can make our theory agree with previous work for laterally unconfined flow in by simply ignoring lateral drag in Eq. (): $${\displaystyle }(H\mathrm{\Sigma }{)}_X-(\mathrm{1}-r)H{H}_X=\mathrm{0}.$$ Integrating and applying the boundary condition Eq. () shows that extensional stress at the grounding line is simply given by $${\displaystyle }{\mathrm{\Sigma }}_{\mathrm{g}}=(\mathrm{1}-r)H_{\mathrm{f}}/\mathrm{2}.$$ This is the same extensional stress at the grounding line as we would expect in the case of the FL model. The flux increases monotonically with floatation thickness when this is substituted into Eq. () and $\mathrm{\Gamma }$ is assumed large: $$\begin{array}{cc}{\displaystyle }Q& {\displaystyle }\approx {\left({\displaystyle \frac{\mathrm{1}-r}{\mathrm{8}}}\right)}^{n^{\mathrm{2}}/(n+\mathrm{1})}{\displaystyle \frac{H_{\mathrm{f}}^{(n^{\mathrm{2}}+\mathrm{3}n+\mathrm{1})/(n+\mathrm{1})}}{{\left(\mathrm{\Gamma },+,W^{-(n+\mathrm{1})/n},H_{\mathrm{f}}\right)}^{n/(n+\mathrm{1})}}}\\ & {\displaystyle }& {\displaystyle }\approx {\left({\displaystyle \frac{\mathrm{1}-r}{\mathrm{8}}}\right)}^{n^{\mathrm{2}}/(n+\mathrm{1})}{\displaystyle \frac{H_{\mathrm{f}}^{(n^{\mathrm{2}}+\mathrm{3}n+\mathrm{1})/(n+\mathrm{1})}}{{\mathrm{\Gamma }}^{n/(n+\mathrm{1})}}}.\end{array}$$ Since we are assuming that $m=\mathrm{1}/n$, this is actually nothing more than a scaled version of Eq. (3.51) in . As Fig. a shows, we do get agreement between the CD calving law results and the FL model for large basal friction coefficients, at least while $H_{\mathrm{f}}/\mathrm{\Lambda }$ remains close enough to the critical value of 2: the flux curves then agree well with each (as indicated by the arrow labelled “driving stress dominant”).

For larger $H_{\mathrm{f}}/\mathrm{\Lambda }$, this agreement ceases. The shelf gets long enough that, even with large enough basal friction on the grounded portion, lateral drag on the floating shelf cannot be ignored. The next section describes in more detail the mechanics of a very long ice shelf.

The finite flux limit for large flotation thickness

For a fixed value of $\mathrm{\Lambda }$ and large flotation thickness $H_{\mathrm{f}}$, the flux $Q$ appears to approach a finite limit in Fig. a. That limit is of the form $$Q\sim (\mathrm{1}-r{)}^{n}C(n\left)W\right(\mathrm{2}\mathrm{\Lambda }{)}^{n+\mathrm{1}},$$ with $C\approx \mathrm{8.67}\times {\mathrm{10}}^{-\mathrm{3}}$ for $n=\mathrm{3}$, regardless of the choice of basal drag parameter $\mathrm{\Gamma }$. The physics behind this are relatively simple: a large value of $H_{\mathrm{f}}$ corresponds to a large difference between ice thickness at the grounding line and at the calving front, which in turn requires a long ice shelf. With a long ice shelf, most of the floating ice shelf becomes fully “buttressed”, in the sense that extensional stress gradients are weak and there is a balance between driving stress and lateral drag as well as basal drag in the grounded part of the glacier. In other words, a local force balance persists around the grounding line, and extensional stresses only become significant in the floating shelf close to the calving front. Moreover, since we are assuming that the floating shelf is still short compared with the length of the glacier and do not include basal melting in our model, ice flux also varies insignificantly along the floating shelf. Hence, the flux through the grounding line is the same as the flux through the calving front. Importantly, we now have a situation in which ice flux through the grounding line is determined entirely by the calving parameter $\mathrm{\Lambda }$ and independent of depth to bedrock at the grounding line.

This situation was previously explored by and . These authors find that ice flux through the calving front is determined in a boundary layer around the calving front in which extensional stress is significant. In our notation, the boundary layer takes exactly the same form as Eq. (9) for floating ice ($\mathit{\theta }=\mathrm{0}$) but with different matching conditions: $$\begin{array}{cc}{\displaystyle }& {\displaystyle }\mathrm{4}\left(H\right|U_X{|}^{\mathrm{1}/n-\mathrm{1}}{U}_X{)}_X-W^{-(n+\mathrm{1})/n}H|U{|}^{\mathrm{1}/n-\mathrm{1}}U\\ & {\displaystyle }& {\displaystyle }-(\mathrm{1}-r)H{H}_X=\mathrm{0},& {\displaystyle }& {\displaystyle }(UH{)}_X=\mathrm{0}\end{array}$$ for $X<X_{\mathrm{c}}$, with $$\begin{array}{cc}{\displaystyle }& {\displaystyle }\mathrm{4}H|U_X{|}^{\mathrm{1}/n-\mathrm{1}}{U}_X=(\mathrm{1}-r)H^{\mathrm{2}}/\mathrm{2},\\ & {\displaystyle }& {\displaystyle }H=H_{\mathrm{c}}\text{at }X=X_{\mathrm{c}},\end{array}$$ where $H_{\mathrm{c}}$ is the prescribed calving front thickness, and $$\begin{array}{cc}{\displaystyle }& {\displaystyle }UH\to Q,\\ {\displaystyle }& {\displaystyle }{W}^{-\mathrm{1}/n-\mathrm{1}}Q|U{|}^{\mathrm{1}/n-\mathrm{1}}\sim -(\mathrm{1}-r\left)\right(Q/U\left)\right(Q/U{)}_X,\\ & {\displaystyle }& {\displaystyle }U\to \mathrm{0}\text{as}X\to -\mathrm{\infty }.\end{array}$$ The analysis of this boundary layer (a formal derivation of which is included in the Supplement) is much the same as for Eq. (9), and $Q$ satisfies a power-law relationship with ice thickness $H_{\mathrm{c}}$ and channel width $W$ at the calving front, of the form $Q\propto (\mathrm{1}-r{)}^{n}C(n)W{H}_{\mathrm{c}}^{n+\mathrm{1}}$. In the CD model, the ice thickness at a floating calving front is $H_{\mathrm{c}}=\mathrm{2}\mathrm{\Lambda }$, which gives the flux relationship Eq. (). A more detailed derivation of the relationship Eq. () is given in the Supplement, and a numerical value computed from the reduced boundary layer model presented there agrees very well with that given above, obtained from the solutions to the full boundary layer problem Eq. (9) for large values of $H_{\mathrm{f}}/\mathrm{\Lambda }$. In fact, convergence to that value is very rapid as $H_{\mathrm{f}}/\mathrm{\Lambda }$ increases, as shown in Fig. a. Figure  shows that, for even moderately large values of $H_{\mathrm{f}}/\mathrm{\Lambda }$ (so when the ice shelf still has relatively limited extent), the thickness and stress profile near the calving front are well approximated by the solution to Eq. (25), regardless of the amount of basal friction in the grounded part of the glacier. The fluxes in all three examples shown in Fig.  are almost identical.

Consider the special case of no basal drag on the grounded part of the glacier. We can show how Eq. () confirms that we expect an anomalous flux-depth-to-bedrock relationship due to buttressing in the ice shelf. Take a grounded calving cliff just at flotation, with thickness $H_{\mathrm{f}}=H_{\mathrm{c}}$. The flux is given by Eq. () with $\mathrm{\Gamma }=\mathrm{0}$. Compare this with the flux through a long floating ice shelf that terminates in a calving cliff of the same thickness $H_{\mathrm{c}}$. The solution Eq. () then predicts that the flux through the floating shelf is smaller than through the grounded calving front, even though the two have the same thickness. This is true at least when the density difference $\mathrm{1}-r$ is small, because the exponent $n^{\mathrm{2}}/(n+\mathrm{1})$ on $(\mathrm{1}-r)$ in the formula for flux for the grounded cliff in Eq. () is smaller than the exponent $n$ for the floating cliff in Eq. (). This provides further evidence for the buttressing action of the ice shelf leading to an anomalous flux–flotation-thickness relationship.

Solutions for a long shelf, same plotting scheme and parameter values as in Fig. . The dashed line shows the solution with $\mathrm{\Gamma }=\mathrm{125}$ and $H_{\mathrm{f}}=\mathrm{3.9}$, the solid line shows $\mathrm{\Gamma }=\mathrm{0}$ and $H_{\mathrm{f}}=\mathrm{3}$. Both are marked with red diamonds in Fig. a. The red line shows the solution to the near-calving-front boundary layer Eq. (25). All three solutions agree closely with this boundary layer near the calving front.

[Figure omitted. See PDF]

Discussion and conclusions

In this paper, we have applied the boundary layer analysis of to a model for channelized outlet glacier flow, incorporating a parameterized description of lateral drag and a simple calving law due to . The purpose of this work has been to show how calving and lateral drag can potentially combine to produce a very different relationship between ice flux at the grounding line and glacier bed geometry from that for laterally unconfined marine ice sheet flow. For the latter, ice flux is an increasing function of depth to bedrock, while for a channelized outlet glacier we find that an “anomalous” relationship in which flux decreases with increasing depth to bedrock is possible (Fig. a above).

Such an anomalous relationship has significant consequences for stable glacier margin positions. Consider the model Eq. () for the flow of the glacier as a whole. Two boundary conditions apply at the free boundary $x=x_{\mathrm{g}}$. One of these is a thickness condition, while the second is the flux condition Eq. (), which can be written in the form $${\displaystyle }q\left(x_{\mathrm{g}}\right(t),t)=Q_{\mathrm{g}}\left(-,b,\left(,x_{\mathrm{g}},\right(,t,\left),\right),,,w,\left(,x_{\mathrm{g}},\right(,t,\left),\right),,,\mathit{\lambda },,,\mathit{\gamma },,,\mathit{\epsilon },,,m,,,n,,,r\right),$$ where $q=-w^{n+\mathrm{1}}{h}^{n+\mathrm{1}}(h+w^{(n+\mathrm{1})/n}\mathit{\gamma }{)}^{-n}|(h+b{)}_x{|}^{n-\mathrm{1}}(h+b{)}_x$ is ice flux, and $Q_{\mathrm{g}}$ is the flux $Q$ predicted by the boundary layer problem, written in terms of the original dimensionless parameters and the channel geometry at the grounding line.

Steady states can now be computed easily from Eq. (). To determine their stability, the theory of and can be extended straightforwardly to the present case, the only modification required being the generalization of the thickness variable $h$ and flux variable $Q$ in to our $wh$ and $wq$, respectively. It then follows that a steady state is linearly stable if and only if (see, e.g., condition (5.1) in ) $$\begin{array}{cc}{\displaystyle }& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}{x}_{\mathrm{g}}}}\left[w,\left(,x_{\mathrm{g}},\right),Q_{\mathrm{g}},\left(-,b,\left(,x_{\mathrm{g}},\right),,,w,\left(,x_{\mathrm{g}},\right),,,\mathit{\lambda },,,\mathit{\gamma },,,m,,,n,,,r\right)\right]\\ & {\displaystyle }& {\displaystyle }>w\left(x_{\mathrm{g}}\right)a\left(x_{\mathrm{g}}\right).\end{array}$$ If $Q_{\mathrm{g}}$ does not necessarily increase with $b$, then steady grounding lines located on downward-sloping beds can become unstable. This is illustrated in Fig. , where steady grounding line positions on a downward-sloping bed are plotted against the accumulation rate $a$ over the ice sheet, which is assumed to be spatially uniform, as is channel width $w$. The steady-state grounding line position of Eq. () is defined implicitly by $${\displaystyle }wQ\left(-,b,\left(,x_{\mathrm{g}},\right),,,w,,,\mathit{\lambda },,,\mathit{\gamma },,,m,,,n,,,r\right)=aw{x}_{\mathrm{g}}.$$ Treating $x_{\mathrm{g}}$ as a function of $a$ and differentiating both sides with respect to $a$, we have $${\displaystyle }\left[{\displaystyle \frac{\mathrm{d}}{\mathrm{d}{x}_{\mathrm{g}}}},\left[w,Q_{\mathrm{g}},\left(-,b,\left(,x_{\mathrm{g}},\right),,,w,\left(,x_{\mathrm{g}},\right),,,\mathit{\lambda },,,\mathit{\gamma },,,m,,,n,,,r\right)\right],-,a,w\right]{\displaystyle \frac{\mathrm{d}{x}_{\mathrm{g}}}{\mathrm{d}a}}=w{x}_{\mathrm{g}},$$ and when the stability condition Eq. () is satisfied, $\mathrm{d}{x}_{\mathrm{g}}/\mathrm{d}a>\mathrm{0}$, so that a stable grounding line must advance when $a$ is increased: a stable ice sheet gets larger when it receives more surface snowfall. Figure  shows a solution branch for which this is not the case, even though the grounding line is located on a downward-sloping bed.

Conversely, we may see grounding lines attain stable steady-state positions on upward-sloping beds if $Q_{\mathrm{g}}$ decreases with depth to bedrock $-b$: Fig.  shows several examples in which the steady-state grounding line advances up a reverse bed slope as accumulation rates are increased. A second mechanism by which such stabilization on upward-sloping beds can occur is the dependence of discharge $w{Q}_{\mathrm{g}}$ on width $w$: a sufficiently narrow bottleneck in the channel could stabilize a grounding line on an upward slope even if $Q_{\mathrm{g}}$ did increase with depth $-b$, because $w{Q}_{\mathrm{g}}$ is an increasing function of $w$ this argument is due to. This second mechanism is however not responsible for the behaviour shown in Fig. , where channel width is constant along the domain. It is worth noting that simulations of Greenland outlet glaciers using the CD calving law have similarly produced steady states located on upward-sloping beds. Our work suggests that this may be due not only to narrowing of the channel but also to the calving law.

Our aim has not been to be authoritative in establishing the existence of an anomalous flux–depth relationship: our model contains at least three components that can be improved upon. First, the parameterized description of lateral drag should eventually be dispensed with, replacing our model with one that resolves the cross-channel dimension. The scaling that underlies our boundary layer model should still be applicable in that case, but the actual boundary layer model will consist of a set of coupled partial differential equations (as opposed to ordinary differential equations) and is likely to be much more onerous to solve for a large number of parameter combinations, as we have been able to do here.

Second, we have neglected the effect of basal melting on the shelf here. This is tractable in the framework we have developed here with a simple, prescribed basal melt rate, but doing so still introduces sufficient complications to lie beyond the scope of a single paper; a second manuscript that incorporates melting into our theory is in preparation.

Third, the calving law we have employed is relatively poorly constrained by observation and is based on a number of simple assumptions about how cracks form near a calving front. Furthermore, it relies entirely on water depth in surface crevasses as a control parameter that should itself be determined by additional physics governing the drainage of surface meltwater.

We have chosen to take the calving model at face value, simply prescribing the crevasse water depth as a control parameter. This is worth emphasizing as the dependence of calving cliff height on flotation thickness predicted by the calving model turns out to be key to the anomalous flux–depth relationship. It is likely that other, more sophisticated calving models (for instance one based on the formulation in ) can also be written in the form of a calving cliff height as a function of crevasse water depth, though presumably with a different specific from the CD model: as in the latter, surface hydrology becomes a key component in understanding calving.

For a floating ice shelf, calving cliff height in the CD model is simply proportional to crevasse water depth and independent of depth to bedrock. In other words, the CD model can then be thought of as a generic calving model that imposes a fixed thickness at the floating glacier terminus. Moving the grounding line to a location with greater flotation thickness (or, equivalently, depth to bedrock) therefore leads to a longer ice shelf forming before it can reach the prescribed calving cliff height. If the mechanical effect of the ice shelf is primarily to provide lateral drag, then a longer shelf leads to a greater reduction in extensional stress between calving front and grounding line, and therefore to lower ice flux despite a greater depth to bedrock at the grounding line. Whether this occurs or not is a function of basal drag on the grounded part of the glacier: if basal drag upstream of the grounding line is moderate compared with lateral drag, then the surface slope and driving stress of the floating shelf will be small, so the effect of the shelf is mostly to generate lateral drag. By contrast, if basal drag is large upstream of the grounding line, then the floating shelf will be relatively steeply sloped and lateral drag will play a lesser role in force balance there, leading to the possibility that the floating shelf does not cause a reduction in extensional stress and hence flux through the grounding line. The changeover between the two regimes happens when, in the notation of Sect. , the basal drag coefficient is approximately (see Eq.  above) $${\displaystyle }\mathit{\gamma }\sim (\mathrm{1}-r{)}^{-\mathrm{1}}{w}^{-(n+\mathrm{1})/n}b,$$ and $\mathit{\lambda }\sim b$.

As we have indicated, the thickness of floating calving fronts in the CD model is uniquely controlled by the crevasse water depth parameter and does not depend on depth to bedrock. The same generic relationship between ice flux and depth to bedrock at the grounding line will therefore be obtained for any other calving law that fixes the height of a floating calving front independently of depth to bedrock. By contrast, the CD model results are unlikely to be robust in the same way for grounded calving fronts. Specifically, for a grounded calving front, the calving law predicts that calving cliff height decreases relatively slowly when the calving front is moved to a location with shallower depth to bedrock. In turn, this leads to more of the calving cliff being exposed above the water line and consequently to larger extensional stresses acting on the calving front, and these larger extensional stresses cause ice flux to increase as depth to bedrock is decreased.

This contrasts with an alternative “calving-at-flotation” (FL) calving law, in which calving front height is always proportional to depth to bedrock and no floating shelf forms. In that case, extensional stress at the grounding line increases with depth to bedrock, and so does ice flux.

We close by noting that our approach can potentially be used to study the effect of other calving laws relatively simply in future, by replacing the function that specifies ice thickness at the calving front. Since an anomalous flux-depth-to-bedrock relationship may be possible and would have significant consequences for stable outlet glacier configurations, it may be worth testing this before embarking on simulations of actual glaciers using different calving laws.

The MATLAB code used in the computations reported is included in the Supplement.

A note on direct solutions of the steady-state problem

At issue is the uniqueness of the orbit that emerges from the fixed point of the dynamical system Eq. (6), at which $(h,\mathit{\sigma },q,\mathit{\chi })=(h_{\mathrm{0}},[a\left(\mathrm{0}\right)/h_{\mathrm{0}}{]}^{\mathrm{1}/n},\mathrm{0},\mathrm{0})$: only when uniqueness is guaranteed does the shooting method of Sect.  make sense. Linearizing the dynamical system around the fixed point leads to a problem with eigenvalues $\mathrm{0}$ (repeated), $a\left(\mathrm{0}\right)w\left(\mathrm{0}\right)>\mathrm{0}$ and $${\displaystyle }-\left[\mathrm{4}{\mathit{ϵ}}^{-\mathrm{1}}n{h}_{\mathrm{0}}\right(a\left(\mathrm{0}\right)/h_{\mathrm{0}}{)}^{\mathrm{1}/n}-(n-\mathrm{1})\left]a\right(\mathrm{0}\left)w\right(\mathrm{0}).$$ The sign of this last eigenvalue is negative when $\mathit{ϵ}$ is small enough. In that case, the fixed point has a two-dimensional centre manifold, a stable manifold and an unstable manifold . The centre manifold has no dynamics (it consists of other fixed points, corresponding to different values of $h_{\mathrm{0}}$, or orbits that do not satisfy ${lim}_{\mathit{\eta }\to -\mathrm{\infty }}\mathit{\chi }=\mathrm{0}$), so that the fixed point must be approached in the limit $\mathit{\eta }\to -\mathrm{\infty }$ along the unique unstable manifold.

The boundary layer problem

The simplified forms of the flux law in Eq. (10) can be derived by a transformation of the boundary layer problem Eq. (9), using $$\begin{array}{cc}{\displaystyle }& {\displaystyle }U=Q^{n/(n+\mathrm{1})}{W}^{\mathrm{1}/(n+\mathrm{1})}\mathcal{U},H=Q^{\mathrm{1}/(n+\mathrm{1})}{W}^{-\mathrm{1}/(n+\mathrm{1})}\mathcal{H},\\ & {\displaystyle }& {\displaystyle }X=W\mathcal{X}\end{array}$$ and $$\begin{array}{cc}{\displaystyle }& {\displaystyle }\mathcal{C}=\mathrm{\Gamma }{W}^{(m+n+\mathrm{3})/(n+\mathrm{1})}/Q^{(\mathrm{2}-nm)/(n+\mathrm{1})},\\ {\displaystyle }& {\displaystyle }{\mathcal{H}}_{\mathrm{f}}=H_{\mathrm{f}}(W/Q{)}^{\mathrm{1}/(n+\mathrm{1})},\mathcal{L}=\mathrm{\Lambda }(W/Q{)}^{\mathrm{1}/(n+\mathrm{1})},\\ & {\displaystyle }& {\displaystyle }{X}_{\mathrm{c}}=W{\mathcal{X}}_{\mathrm{c}}.\end{array}$$ With these definitions, it is easy to show that Eq. (9) is invariant under the transformation $$\begin{array}{cc}{\displaystyle }& {\displaystyle }(U,H,X,X_{\mathrm{c}},H_{\mathrm{g}},W,Q,\mathrm{\Gamma },\mathrm{\Lambda })\\ {\displaystyle }& {\displaystyle }\mapsto (\mathcal{U},\mathcal{H},\mathcal{X},{\mathcal{X}}_{\mathrm{c}},{\mathcal{H}}_g,\mathrm{1},\mathrm{1},\mathcal{C},\mathcal{L}).\end{array}$$ The parameters in this rescaled version of the model are $r$, $n$, $m$, ${\mathcal{H}}_{\mathrm{f}}$, $\mathcal{C}$ and $\mathcal{L}$, while $\mathcal{U}$ and $\mathcal{H}$ are the dependent variables, and $\mathcal{X}$ is the independent variable. It then can be shown that the transformed boundary layer problem has a solution if and only if the parameters $r$, $n$, $m$, ${\mathcal{H}}_{\mathrm{f}}$, $\mathcal{C}$ and $\mathcal{L}$ satisfy some functional relationship with each other. Using this fact, it is easy to show that the simplified flux laws Eq. (10) must hold.

Deriving that functional relationship between $r$, $n$, $m$, ${\mathcal{H}}_{\mathrm{f}}$, $\mathcal{C}$ and $\mathcal{L}$ can be done by a further coordinate transform of the dependent variables (see also Appendix A in ) $$\begin{array}{cc}{\displaystyle }& {\displaystyle }\mathcal{Q}=\mathcal{U}\mathcal{H},\mathrm{\Psi }={\mathcal{Q}}^{-\mathrm{1}}{\mathcal{U}}^{-(\mathrm{2}n+\mathrm{1})/n^{\mathrm{2}}}|{\mathcal{U}}_{\mathcal{X}}{|}^{\mathrm{1}/n-\mathrm{1}}{\mathcal{U}}_{\mathcal{X}},\\ & {\displaystyle }& {\displaystyle }\mathit{\xi }={\mathcal{U}}^{(n+\mathrm{1}{)}^{\mathrm{2}}/n^{\mathrm{2}}}\end{array}$$ and switching to an independent variable $\mathit{\zeta }$ defined through $${\displaystyle }\mathit{\zeta }=\underset{\mathrm{0}}{\overset{\mathcal{X}}{\int }}\mathcal{Q}\left({\mathcal{X}}^{\prime }{)}^{-\mathrm{1}}\mathrm{\Psi }\right({\mathcal{X}}^{\prime }{)}^n\mathit{\xi }({\mathcal{X}}^{\prime }{)}^{-\mathrm{1}/(n+\mathrm{1})}\mathrm{d}{\mathcal{X}}^{\prime }.$$ This transforms the boundary layer problem into a non-singular dynamical system in which the matching conditions Eq. () correspond to a fixed point $(\mathrm{\Psi },\mathit{\xi },\mathcal{Q})=(\mathrm{1},\mathrm{0},\mathrm{1})$ being attained as $\mathit{\zeta }\to -\mathrm{\infty }$, and there is a unique orbit along which this can happen; to prove the uniqueness of that orbit, an additional transformation to $\mathit{\nu }={\mathit{\xi }}^{\left(n^{\mathrm{2}}\right(m+\mathrm{1})-n)/\left(\mathrm{2}\right(n+\mathrm{1}{)}^{\mathrm{2}})}$ may be required, but the basic argument remains the same as in . The boundary conditions Eqs. ()–() then provide two further constraints: one on the location of the point along the orbit that corresponds to the calving front, and another to relate the parameters of the model, leading to a functional relationship between $r$, $n$, $m$, ${\mathcal{H}}_{\mathrm{f}}$, $\mathcal{C}$ and $\mathcal{L}$. In practice, this must be solved for numerically by integrating the transformed dynamical system. More complete details on the solution of the boundary layer problem can be found in the Supplement and the numerical code provided therein.

The Supplement related to this article is available online at https://doi.org/10.5194/tc-11-2283-2017-supplement.

The authors declare that they have no conflict of interest.

Acknowledgements

Christian Schoof acknowledges the support of NSERC grants 357193-13 and 446042-13. Andrew D. Davis's work is supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, under award number DE-SC0010518. Tiberiu V. Popa was supported by the NSERC CREATE AAP programme. We are grateful to the three anonymous referees and the editor for their thorough scrutiny. Edited by: Olivier Gagliardini Reviewed by: three anonymous referees