The drainage of glacier and ice sheet surface lakes

We consider a surface melt water stream with cross-sectional area $S$, with the base of the stream channel at an elevation $b$, and the velocity in the stream being $u$. Let $x$ be distance along the stream and $t$ time, and let $S$, $u$ and $b$ depend on $x$ and $t$ (figure 1). Assuming a Darcy-Weisbach law governing shear stress at the walls of the channel, we express conservation of mass and momentum using a St Venant model as [e.g. 53, chapter 4]

	$\displaystyle\rho_{w}\left[S_{t}+(uS)_{x}\right]=$	$\displaystyle\rho m,$		(1a)
	$\displaystyle\rho_{w}S(u_{t}+uu_{x})=$	$\displaystyle-\frac{\rho_{w}fu^{2}P(S)}{8}-\rho_{w}gS\left[b_{x}+h(S)_{x}\right]$		(1b)
where subscripts $x$ and $t$ denote partial derivatives, $m$ is melt rate at the channel wall, expressed as an area of ice melted per unit time and unit length of channel, $P(S)$ is the wetted perimeter of the channel and $h(S)$ is the elevation of the water surface above the bottom of the channel. $\rho_{w}$ is the density of water and $f$ is a friction coefficient depending on wall roughness in the channel. Note that by equating the source term $m$ with melt rate, we ignore seepage into or out of a firn aquifer, or substantial water input from tributary streams.

To simplify matters, we assume that the cross-sectional area can grow or shrink but retains a shape determined by its size alone. Our main interest is in downward incision of the channel, which we assume to be a slow process compared with the adjustment of channel shape, since we will assume below that water depth is much less than the typical amplitude of channel elevation $b$. Consequently, we treat $P$ and $h$ as non-decreasing functions of $S$, whose form depends on the geometry of the cross-section. At a minimum, we know that water depth must vanish when cross-sectional area does, so $h(0)=0$.

The simplest way to parameterize the cross-sectional shape of the channel is to treat it as a semi-circle (figure 2). In that case, the radius $r$ of the cross-section is $r=(2\pi^{-1}S)^{1/2}$ and

$$P(S)=\pi r=(2\pi S)^{1/2},\qquad h(S)=r=(2\pi^{-1}S)^{1/2}$$

(1c)

Alternatives would be to assume the channel is triangular with a fixed angle $\theta$ between the channel sides and the vertical

$$P(S)=\frac{2S^{1/2}}{\sin(2\theta)^{1/2}},\qquad h(S)=\frac{S^{1/2}\cos(\theta)}{\sin(2\theta)^{1/2}}$$

(1d)

or to fix a width $W$ that is much greater than water depth, and to put [as is done in 53]

$$P(S)=W,\qquad h(S)=\frac{S}{W}.$$

(1e)

Generically, this suggests we consider

$$P(S)=c_{1}S^{\alpha},\qquad h(S)=c_{2}S^{\beta}$$

(1f)

with $c_{1},\,c_{2}>0$, $\alpha\geq 0$, $\beta>0$ (we admit that width and therefore wetted perimeter may not depend on $S$, but water depth must, so $\beta$ cannot vanish while $\alpha$ can): (1c)–(1d) have $\alpha=\beta=1/2$ while (1e) puts $\alpha=0$, $\beta=1$. In fact, the examples above suggest that the product of wetted perimeter and water depth scale as the cross-sectional area, in which case $\alpha+\beta=1$, and that the exponents are not only positive but also satisfy $0\leq\alpha<1$, $0<\beta\leq 1$.

We assume that energy dissipated by the flow is instantly transferred to the wall of the channel and turned into latent heat, and that this is the dominant mechanism of channel incision. A more sophisticated model could track the temperature of the water and use a heat transfer model [see also the discussion in 54]; here we assume that heat transfer is highly efficient at the length scales under consideration. Letting $\mathcal{L}$ be the latent heat of fusion per unit mass of water, we put

$$\rho\mathcal{L}m=-\frac{\rho_{w}fu^{3}P(S)}{8},$$

(1g)

the right-hand side being the rate at which work is done per unit length of channel by water moving at velocity $u$ against the friction force $\rho fu^{2}P(S)/8$ on the channel wall. In order to model how fast the channel cuts into the ice, we assume that downward incision can be estimated by distributing melt equally over the wetted perimeter, leading to an incision rate of $m/P$. Future work will need to address both the channel shape parameterizations and the distribution of melt over the channel wall: related work on englacial channels [55] may serve as a template.

In addition, we assume that the ice surface is moving horizontally at a velocity $U$ and is subject to localized uplift or drawdown at a prescribed rate $w(x)$ due to flow of the glacier or ice sheet over bed topography [e.g. 47], where we will later assume for simplicity that $U$ is constant in space as well as time, as is appropriate for instance for rapidly sliding ice. Then

$$b_{t}+Ub_{x}=w-\frac{m}{P(S)}.$$

(1h)

We assume that the base of the channel is incised into an ice surface at elevation $s$, with $b\leq s$. In assuming that $w$ is constant in time, we are assuming not only that we can ignore localized, enhanced ‘creep closure’ around a deeply incised channel [56] as well as snow accumulation at the base of the channel during winter, but also that $s$ is in steady state [47], and itself satisfies

$$Us_{x}=w.$$

(1i)

We will generally use $b(x,0)=s(x)$ as an initial condition, representing a channel that is only just beginning to incise into the ice surface. A modification of the present model to a dynamically evolving ice surface $s$ will be presented elsewhere.

We envisage the channel draining a reservoir at its upstream end. For simplicity, we assume that the seal point of the lake (the maximum in $b$) is some distance downstream of $x=0$, and that we can relate lake volume directly to water level $h$ at $x=0$, and put

$$\dot{V}=q_{0}(t)-\left.(uS)\right|_{x=0},\qquad V(t)=V_{L}(h(S(0,t))),$$

(1j)

with the dot on $\dot{V}$ denoting an ordinary time derivative, $q_{0}$ being a prescribed rate of inflow to the lake due to surface melting in some larger upstream catchment. $V_{L}$ is an increasing function of $h$, dictated by the bathymetry of the lake. The bathymetry in turn is presumably determined by $U$ and $w$, but we do not consider that in detail here, nor do we consider the possibility that surface loading due to the lake could affect the motion of the ice. As with an evolving ice surface $s$, the latter complication will be studied in a separate paper.

Before we continue, we note some important limitations of the model. First, by assuming a fixed flow path and not modelling tortuosity, we are not considering the effect of meanders on flow and channel incision, even though meandering is known to be a common feature of glacier surface streams [e.g. 57, 58]. Second, the one-dimensional nature of the model implies not only that there is a single outflow from the lake, which is ultimately likely as two competing outflow channels are presumably prone to instability, with the larger channel persisting while the smaller is abandoned. It also implies that, if flow in the channel were to cease temporarily due for instance to seasonal variations in water supply $q_{0}$, the same channel will be re-occupied when flow recommences. We return to this in section 5.

In addition, we also neglect surface lowering due to melt driven by insolation or a warm atmosphere, or freezing due to heat fluxes into the ice. This may be reasonable for the incision of the channel relative to the rest of the ice surface $s$, is more questionable for the lake itself. Here, enhanced absorption of incoming radiation in the lake water is likely to lead to preferential melting of the deeper portions of the lake [see also 59]. By the same token, we also neglect the possibility that the lake water could be warmed relative to the melting point by incoming solar radiation [see also 35]. That said, by omitting externally-driven melting, our model allows us to focus purely on the coupled effects of ice and water flow in the erosion of the channel and its effect on lake drainage.

We assume that a length scale $[x]$ can be determined from the uplift field $w(x)$, and that scales for vertical and horizontal ice velocities $[w]$ and $[U]$ are also known. In terms of these, we define scales $[t]$, $[S]$, $[u]$, $[b]$ and $[V]$ through

$$\frac{\rho_{w}g[b][S]}{[x]}=\frac{\rho_{w}f[u]^{2}P([S])}{8},\qquad\frac{\rho_{w}f[u]^{3}}{8\rho\mathcal{L}}=[W]=\frac{[U][b]}{[x]},\qquad[U][t]=[x].$$

Our choice of scales here reflects the following: we are interested in significant channel incision over a single advective time scale $[t]=[x]/[U]$ for the ice surface, so that uplift, advection and incision of the channel naturally compete with each other. Given a natural surface topography $[b]=[w][x]/[U]$, that sets a dissipation rate and therefore a water flux scale $[u][S]$. It is possible to construct the same problem as below for a much shorter channel incision time scale, corresponding to greater dissipation and therefore water fluxes; that however precludes the generation of a lake, which requires the lake seal to be generated through uplift of the ice surface.

From these scales we obtain the dimensionless groups

$$\nu=\frac{h([S])}{[b]},\qquad Fr^{2}=\frac{[u]^{2}}{gh([S])},\qquad\varepsilon=\frac{g[b]}{\mathcal{L}},\qquad\delta=\frac{[U]}{[u]}.$$

(2)

These have straightforward interpretations: $\nu$ is the ratio of water depth to ice surface topography, the Froude number $Fr$ is the usual square root of the ratio of kinetic to gravitational energy, $\varepsilon$ is the ratio of gravitational potential energy to latent heat, and $\delta$ is the ratio of ice to water velocity. With the possible exception of $Fr$, we expect all of these parameters to be small: if water moved at speeds comparable to the ice, then surface drainage would presumably be of no interest, while the surface topography scale would have to be around 30 km with a terrestrial gravitational field $g\approx$ 10 m s^-1 in order for gravitational potential energy and latent heat $\mathcal{L}\approx 3.35\times 10^{5}$ J kg^-1 to be comparable. We also expect the water depth in a glacially-dammed lake to be larger than the flow depth in the stream draining it, except possibly during a very rapid outburst flood or for shallow lakes.

In fact, for realistic values of $[U]=100$ m a^-1, $[b]=10$ m, $[x]=1$ km, $g=9.8$ m s^-2, $f=0.05$, we obtain, with $h(S)$ and $P(S)$ given by (1c)

$$[u]=1.2\,\mbox{m s}^{-1},\qquad[S]=0.47\,\mbox{m}^{2},$$

values that are realistic for surface streams with gentle $[b]/[x]\approx 0.01$ slopes. With the choice of scales defined through (2), we define dimensionless variables through

$$x=[x]x^{*},\qquad t=[t]t^{*},\qquad u=[u]u^{*},\qquad S=[S]S^{*},\qquad b=[b]b^{*},$$

(3)

and define

$$P^{*}(S^{*})=\frac{P(S)}{P([S])},\qquad h^{*}(S^{*})=\frac{h(S)}{h([S])},$$

(4)

and also put $U=[U]U^{*}$, $w=[w]w^{*}$. Then, in dimensionless form, dropping the asterisks on the dimensionless variables immediately, the model becomes

	$\displaystyle\delta S_{t}+(uS)_{x}=$	$\displaystyle\varepsilon u^{3}P(S),$		(5a)
	$\displaystyle\nu Fr^{2}S(\delta u_{t}+uu_{x})=$	$\displaystyle-u^{2}P(S)-Sb_{x}-\nu Sh(S)_{x},$		(5b)
	$\displaystyle b_{t}+Ub_{x}=$	$\displaystyle w-u^{3}.$		(5c)

Following the discussion above, we assume that $\delta\ll 1$, $\nu\ll 1$ and $\varepsilon\ll 1$. At leading order in these small parameters, we obtain from (5)

$$(uS)_{x}=0,\qquad u^{2}P(S)=-Sb_{x},\qquad b_{t}+Ub_{x}=w-u^{3}.$$

(6)

Water flux along the channel is constant in space, velocity is controlled by a balance of friction at the channel wall and the downslope component of gravity acting on the water in the channel, and the channel bottom evolves due to advection, uplift, and melting driven by local dissipation of heat in the flow of water.

The reduced model is subject to the caveat that the local Froude number $Fr_{loc}=Fr\,u/(\beta S^{\beta})^{1/2}$ remain less than unity. Where $Fr_{loc}>1$, the channel becomes unstable to bedform formation at short wavelength, while for $Fr_{loc}>2/(1-\alpha)$, roll waves form in the flow (see §3 of the supplementary material, also sections 4.4.4–4.5.2 and chapter 5 of Fowler [53]). A reduced model that does not explicitly resolve these phenomena but focuses on channel incision at the larger scale may still be possible, but would presumably require a multiple scales expansion [60]. We leave this to future work.

Persisting with (6), we find that $q=uS$ is independent of position. Hence $u$ depends on flux $q$ and slope $-b_{x}$ through from (6)₂ as

$$u^{3}q^{-1}P(u^{-1}q)=-b_{x},$$

(7)

where we assume that $b_{x}<0$ and $q>0$. With channel geometry given by (1f), specifically $P(S)=S^{\alpha}$ in dimensionless form, we obtain a dimensionless melt rate

$$M(-b_{x},q):=u^{3}=\left(-q^{1-\alpha}b_{x}\right)^{3/(3-\alpha)}.$$

(8)

The function $M$ here is monotonically increasing and convex function of $-b_{x}$ for $\alpha\geq 0$, strictly so if $\alpha>0$, and a monotonically increasing function of $q$ for $\alpha<1$ (see also figure 3). Our assumptions about channel geometry can be relaxed significantly while allowing these properties to be preserved: as shown in §2 of the supplementary material, monotonicity is assured if hydraulic radius $S/P(S)$ is an increasing function that vanishes when $S=0$, while (strict) convexity follows if $P$ is (strictly) concave.

At face value, substituting in (6)₃ yields the single evolution equation for $b$,

$$b_{t}=w-Ub_{x}-M(-b_{x},q).$$

(9)

Note that this is effectively the stream power equation for landscape evolution in [49, 52, 61], with some alterations as explained in section 3.1.

Our assumption that $b_{x}<0$ is however not always satisfied. Where such reverse slopes occur, the reduced model breaks down: water depths become large, flow velocities become small and melt rates vanish at leading order, and we put $M(-b_{x},q)=0$ if $b_{x}>0$ to account for this. That in itself does not however suffice, since reverse slopes can induce ponding even on downward slopes further upstream. We deal with this next.

Formally, the dynamics of a ‘ponded section’ of channel (where, once more, water is deep while surface slope, velocity and melt rates are small, see figure 1) can be captured formally by a rescaling as described in appendix A. Without engaging in that formalism, assume that water storage in the ponded sections remains negligible, so that we can treat flux $q$ as constant throughout the domain; this turns out to require $\delta\ll\nu^{1/\beta}$. When there is flow with $q>0$, then a ponded section comprises those points that lie below the high point or ‘seal’ in the channel bed at its downstream end, since the nearly flat water surface in the ponded section cannot exceed that seal height significantly. Hence the union of all ponded sections is the set $\{x:b(x,t)<\sup_{x^{\prime}>x}b(x^{\prime},t)\}$, which necessarily encompasses reverse slopes $b_{x}>0$, but will generally also include stretches of channel with a downward-sloping bed $b_{x}<0$.

If no water is flowing and $q=0$, actual ponded sections with standing water may be smaller. The purpose of identifying ponded sections is however mostly to impose an absence of melting there, and with $q=0$, we obtain $M=0$ regardless of whether we have identified ponding correctly. The appropriate modification of (9) to account for ponding is therefore via a ‘ponding function’ $c$,

	$\displaystyle b_{t}=$	$\displaystyle w-Ub_{x}-c(x,t)M(-b_{x},q),$		(10a)
	$\displaystyle c(x,t)=$	$\displaystyle\left\{\begin{array}[]{l l}1&\mbox{if }b(x,t)\geq\sup_{x^{\prime}>x}b(x^{\prime},t),\\ 0&\mbox{otherwise}.\end{array}\right.$		(10d)

Note once more that the introduction of the ponding function is redundant where $b_{x}>0$ in ponded sections, since in that case $M(-b_{x},q)=0$ by definition.

We still need to deal with mass conservation equation (1j) for the lake to determine the flux $q(t)$, which is constant along the channel, but can change over time. Note that we have not rendered (1j) in a leading-order, dimensionless form yet. We do so next.

As we have to revisit the non-dimensionalization of the problem, we temporarily reintroduce asterisks on dimensionless variables. We assume that the lake at $x^{*}=0$ exists because it is upstream of a ponded section with a reverse slope. Let $\hat{h}^{*}=\nu^{-1}h^{*}(S^{*})=h(S)/[b]$ be dimensionless water depth in that ponded section, scaled with ice surface topography $[b]$ (as is appropriate for ponded sections, see appendix A). Also let $h_{0}^{*}=\hat{h}^{*}(0,t^{*})+b^{*}(0,t^{*})$ be the water level of the lake. We define a dimensionless lake volume function and a dimensionless water supply function through

$$\hat{V}(h_{0}(t^{*}))=\frac{V_{L}(h(S(0,t)))}{[u][S][t]},\qquad Q^{*}(t^{*})=\frac{q_{0}(t)}{[u][S]}.$$

(11)

where the variables on the right-hand sides of both equalities are dimensional. We assume formally that $\hat{V}$ and $Q$ are $O(1)$ functions. By this, we mean that lake volume is comparable to (or less than) the volume $[u][S][t]$ typically carried by the channel in a single channel evolution time scale $[t]$, and inflow into the lake is comparable with or less than the flux scale $[u][S]$ that causes significant channel incision over the advective time scale of the ice.

We immediately revert to dropping asterisks on dimensionless variables. As in the previous section, water surface elevation $b+\hat{h}$ must be constant up to the lake seal (the end of the ponded section that extends downsteam from $x=0$), Water surface elevation also cannot exceed seal height at leading order (appendices A and B.2). Consequently, we find that water depth at the upstream end of the domain is either at the height of the seal point if water is flowing, or below that seal height if no water is flowing. We denote the seal height by $b_{m}(t)$, so that

	$$h_{0}\leq b_{m}(t)=\sup_{x>0}b(x,t).$$		(12a)
Similarly, we will use $x_{m}$ to denote the seal location, $x_{m}(t)=\sup\{x:b(x,t)<b_{m}(t)\}$.

With $uS=q$ constant throughout the domain, water balance of the lake $\dot{\hat{V}}=Q-(uS)|_{x=0}$ can therefore be written as

	$\displaystyle\gamma\dot{h}_{0}=$	$\displaystyle Q(t)-q,$		(12b)
	$\displaystyle q=$	$\displaystyle\left\{\begin{array}[]{l l}0&\mbox{if }h_{0}<b_{m},\\ \max\left(Q-\gamma\dot{b}_{m},0\right)&\mbox{if }h_{0}=b_{m},\end{array}\right.$		(12e)

where $\gamma(h_{0})=\mathrm{d}\hat{V}/\mathrm{d}h_{0}$ is storage capacity in the lake and overdots again denote time derivatives. Flux $q$ in the channel is the difference between inflow into the lake and the rate at which water is retained in the lake, and the latter is controlled by how the high point in the channel itself evolves due to uplift, advection, and incision.

If we treat $c(x,t)$ and $q(t)$ momentarily as known, then (10a) can be recognized as being of Hamilton-Jacobi form [49],

$$b_{t}=-\mathcal{H}(x,t,b_{x},q),$$

(13)

where the Hamiltonian $\mathcal{H}$ is given by

$$\mathcal{H}(x,t,p,q)=Up+c(x,t)M(-p,q)-w(x),$$

(14)

replacing $b_{x}$ (itself a function of $x$ and $t$) with $p$ for clarity in the meaning of derivatives of $\mathcal{H}$.

The method of characteristics [62, §3] allows us to write the solution to (10a) as follows: we define characteristics as curves of constant $\sigma$ in the transformation $(\sigma,\tau)\mapsto(x,t)$ given by

$$x_{\tau}=\mathcal{H}_{p}=U-cM_{-p}(-p,q),\qquad t(\sigma,\tau)=\tau,$$

(15)

where $\mathcal{H}_{p}(x,t,p,q)$ is the partial derivative of $\mathcal{H}$ with regard to its third argument, with $x$, $t$ and $p$ all treated as functions of $\sigma$ and $\tau$, while $M_{-p}(-p,q)$ is the partial derivative of $M$ with respect to its first argument. Along a given characteristic, $b(\sigma,\tau)$ and $p(\sigma,\tau)=b_{x}$ evolve as

$$b_{\tau}=-\mathcal{H}+\mathcal{H}_{p}p=w-c[M_{-p}(-p,q)p+M(-p,q)],\qquad p_{\tau}=-\mathcal{H}_{x}=w_{x}$$

(16)

subject to the given initial and boundary conditions. We take these to be $b(x,0)=b_{in}(x)$ at $t=0$ and $b(0,t)=b_{in}(0)$ at $x=0$, so elevation at the upstream end of the domain remains constant throughout. Note that we do not attempt to differentiate the piecewise constant function $c$ when forming $\mathcal{H}_{x}$ in equation (16). Instead, we regard discontinuities in $c$ as potential shocks, and treat them separately.

As already pointed out, the problem (13)–(14) is structurally very similar to ‘stream power models’ for landscape evolution [49]. Aside from complications due to the ponding function $c$, the main differences between our model and canonical stream power models is the way that water flux is controlled by lake drainage, and the fact that the Hamiltonian $\mathcal{H}$ here remains convex in the slope variable $p$, but need not be monotonically decreasing due to the advection term. As a result, characteristics and shocks can travel downstream as well as upstream, with major implications for breaching the seal of the lake and controlling the flux $q$.

We give a comprehensive account of shocks and discontinuities in $c$ below Figure 4 provides an overview of the different possibilities. We treat the majority of cases in sections 3.2–3.3 and derive formulae for flux $q$ in terms of shock geometry at the lake seal in section 3.4. We relegate the analysis of the upstream end of ponded sections to appendix C, where we show that the discontinuity in $c$ at such a location cannot generate a shock but may give rise to an expansion fan rather. The material below is fairly dense and, at this stage, abstract. On a first reading, it may be preferable to skip to section 4 to understand the zoology of features of the solution before filling in the theoretical background in sections 3.2–3.4 and appendix C.

Equation (13) breaks down when characteristics intersect at shocks. Intersections require characteristics to travel faster upstream of the seal than downstream. The melt rate $M$ is convex in slope $p$, and for $c=1$ on either side of a shock in a flowing section, so is the Hamiltonian $\mathcal{H}$. Denoting by superscripts ⁺ and ^- limits taken from above and below the shock at $x=x_{c}(t)$, respectively, we must have $b_{x}^{+}<b_{x}^{-}<0$ for a shock in a flowing section (figure 4(a)). These shocks represent ‘knickpoints’ in standard geomorphological parlance [49, 51].

We require that $b$ remain continuous across any shock or discontinuity in $c$ (see appendix B for the boundary layer structure of the full model around the different types of shock). Differentiating both sides of $b^{-}(x_{c}(t),t)=b^{+}(x_{c}(t),t)$ with respect to $t$, we obtain $b_{t}^{-}+b_{x}^{-}\dot{x}_{c}=b_{t}^{+}+b_{x}^{+}\dot{x}_{c}$, where the overdot denotes differentiation with respect to time. Equivalently (see also Royden and Perron [51])

$$\dot{x}_{c}=\frac{\mathcal{H}(x_{c},t,b_{x}^{+},q)-\mathcal{H}(x_{c},t,b_{x}^{-},q)}{b_{x}^{+}-b_{x}^{-}}=U+\frac{c^{+}M(-b_{x}^{+},q)-c^{-}M(-b_{x}^{-},q)}{b_{x}^{+}-b_{x}^{-}},$$

(17)

where of course $c^{+}=c^{-}=1$ for a shock in a flowing section; we retain $c^{+}$ and $c^{-}$ for later convenience. By the mean value theorem, a strictly convex $\mathcal{H}$ corresponds to $\dot{x}_{c}$ somewhere between the characteristic velocities on either side, with characteristics terminating at the shock from both sides as expected. In fact, knickpoint shocks between two parts of a flowing section can occur only if $\alpha>0$ in (8), so $\mathcal{H}$ are strictly convex for $p<0$: for $\alpha=0$, the characteristic velocity $\mathcal{H}_{p}=U-q$ in a flowing section is independent of slope and characteristics do not cross.

Shocks between a ponded section upstream and a flowing section downstream ($c^{-}=0$, $c^{+}=1$, $b_{x}^{-}>0$, $b_{x}^{+}<0$, see figure 4(b)) are structurally equivalent to knickpoint-type shocks. Equation (17) still holds, where now $c^{-}M(-b_{x}^{-},q)=0$; the discontinuity in $c$ is in fact a red herring her since $M=0$ for $b_{x}>0$ anyway (this differs from the upstream end of a ponded section in appendix C, where the discontinuity in $c$ is crucial). The important distinction with the knickpoint shocks of section 3.2 is that shocks between ponded and flowing sections can form even if $M$ is not strictly convex (that is, for $\alpha=0$), since the characteristic velocity $x_{\tau}^{-}=U$ upstream of the shock is larger than its counterpart $x_{\tau}^{+}=U-M_{-p}(-b_{x}^{+},q)$ downstream, and characteristics terminate at the shock from both sides.

The seal of the lake may take the form of a shock between ponded and flowing, and its motion then controls the flux $q$ as described in section 3.4 below. We refer to ‘breaching’ of the seal as incision into a seal that was previously in steady state, leading flux $q$ to increase and the lake to drain, and this requires a shock to pass the steady seal location as we will show in §4.

The transition from ponded to flowing need not correspond to a shock, however. A continuous slope with $b_{x}^{-}=b_{x}^{+}=0$ is possible if the transition point $x_{s}(t)$ is a local maximum of $b$ such that characteristics enter from one side and exit on the other, with no jump in $b$ or $b_{x}=p$, and with a continuous melt rate $cM$ and Hamiltonian $\mathcal{H}$ (figure 4(c)). Differentiating $b_{x}^{-}(x_{s}(t),t)=b_{x}^{+}(x_{s}(t),t)=0$ with respect to time and differentiating (10a) with respect to $x$, so $b_{xt}^{-}+Ub_{xx}^{-}=w_{x}(x_{s})$ with $c^{-}=0$ and $b_{xt}^{+}+(U-M_{-p}(0^{-},q))b_{xx}^{+}=w_{x}(x_{s})$ with $c^{+}=1$, we obtain

$$\dot{x}_{s}=U-\frac{w_{x}}{b_{xx}^{-}}=U-M_{-p}(0^{-},q)-\frac{w_{x}}{b_{xx}^{+}}.$$

(18)

Since the ponded section is upstream and $x_{s}$ is a maximum of $b$, we have $b_{xx}^{-}<0$ and $b_{xx}^{+}<0$. There are two cases, with $w_{x}$ negative and positive at the seal, respectively. Positive $w_{x}$ corresponds to rapid downslope motion of the seal with $\dot{x}_{s}>U$; this does not occur except for contrived initial conditions.

Assume therefore that $w_{x}<0$, so $\dot{x}_{s}<U$. Characteristics upstream of $x_{s}$ travel at speed $U>\dot{x}_{s}$, so a smooth slope requires characeteristics to emerge on the downstream side, where the characteristic speed is $x_{\tau}^{+}=U-M_{-p}(-p^{+},q)=U-M_{-p}(0^{-},q)$, with $0^{-}$ indicating the limit taken as $p=0$ is approached from below. Requiring $x_{\tau}^{+}>\dot{x}_{s}$ so that characteristics exit downstream, a smooth slope is possible provided

$$M_{-p}(0^{-},q)<\frac{w_{x}}{b_{xx}^{-}}\qquad\mbox{and}\qquad w_{x}<0.$$

(19)

Importantly, this type of smooth seal can correspond to a steady state. The steady state seal location is defined by $w(x_{s})=0$ with $w_{x}(x_{s})<0$. In steady state, $Ub_{x}=w$ upstream of the seal, which ensures that $b_{x}^{-}=0$ at the seal as required, while differentiating both sides yields $w_{x}=Ub_{xx}^{-}$, so $\dot{x}_{s}=0$ as needed for a steady state.

It is instructive to ask when (19) can be violated. For $M$ given by equation (8), $M_{-p}(0^{-},q)=0$ if $\alpha>0$, and (19) will not be violated provided $w_{x}<0$ and $b_{xx}^{-}<0$. By contrast, for $\alpha=0$, $M_{-p}(0^{-},q)=q$ and (19) can be violated for sufficiently large fluxes $q>U$. In practice, this implies that shocks breaching the seal form downstream of the seal for $\alpha>0$ and migrate upstream, while for $\alpha=0$, such knickpoint shocks cannot form and shocks that breach the seal form at the seal itself when $q>U$.

Key to the model is to understand how flux $q$ evolves, which requires the evolution of the seal point height $\dot{b}_{m}$ in (12e). There are two scenarios: either the seal point $x_{m}$ is a shock, in which case from (17)

$$\dot{b}_{m}=b^{-}_{t}+b_{x}^{-}\dot{x}_{m}=w(x_{m})+\frac{b_{x}^{-}M(-b_{x}^{+},q)}{b_{x}^{+}-b_{x}^{-}},$$

(20)

or the seal is a smooth transition point, and with $b_{x}=0$ at such a smooth transition,

$$\dot{b}_{m}=b_{t}+b_{x}\dot{x}_{m},=b_{t}+Ub_{x}=w(x_{m}).$$

(21)

We can now compute flux $q$. Assuming lake level is equal to seal height with $h_{0}=b_{m}$, (12e) leads to

$$q=\max\left(Q-\gamma w(x_{m})-\gamma\frac{b_{x}^{-}M(-b_{x}^{+},q)}{b_{x}^{+}-b_{x}^{-}},0\right)$$

(22)

if there is a shock and (20) holds, or $q=\max(Q-\gamma w(x_{m}),0)$ if the seal is smooth. Here, $Q-\gamma w(x_{m})$ is simply the base outflow rate that results if there is no incision into the seal due to melting.

Only the case (22) of a shock at the seal is non-trivial. Any positive outflow rate $q$ satisfies

$$q+\gamma\frac{b_{x}^{-}M(-b_{x}^{+},q)}{b_{x}^{+}-b_{x}^{-}}=Q-\gamma w(x_{m}).$$

(23)

For $\alpha>0$ and $p<0$, the function $M(-p,q)$ defined in equation (8) is an increasing, concave function of $q$ with $M_{q}(-p,0)=\infty$. With $b_{x}^{-}>0$ and $b_{x}^{+}<0$ for a shock at the seal, it follows that the left-hand side of (23) vanishes at $q=0$, decreases for small $q$, reaches a global minimum and then increases thereafter (figure 5(a)). If the base outflow rate $Q-\gamma w(x_{m})$ is positive, there is then a single positive root for $q$, where $q$ increases with base outflow rate $Q-\gamma w(x_{m})$ and with storage capacity $\gamma$ (for fixed base outflow rate). For $Q-\gamma w(x_{m})\leq 0$, the solution for $q$ can become multivalued: $q=0$ is a valid solution of the original problem (22), but there are in addition two non-zero solutions if $Q-\gamma w(x_{m})<0$ is larger than the global minimum with respect to $q$ of the left-hand side of (23) (figure 5(a)). For the melt rate $M$ given by (8), that situation arises when

$$0>Q-\gamma w(x_{m})\geq-\frac{2\alpha}{3-\alpha}\left(\frac{3(1-\alpha)\gamma b_{x}^{-}}{(3-\alpha)(b_{x}^{-}-b_{x}^{+})}\right)^{(3-\alpha)/(2\alpha)}(-b_{x}^{+})^{3/(2\alpha)}.$$

(24)

For even more negative $Q-\gamma w(x_{m})$, $q=0$ is the only solution. The multivalued solution $q$ in that case results from the fact that the incision rate due to sufficiently large flux can overcome an uplift rate that on its own is sufficient to re-seal the lake, and thus maintain that flux.

To resolve the multivaluedness of $q$ properly, we have to go to higher order and compute the difference between water level $h_{0}$ in the ponded region and seal height $b_{m}$ at leading order. Doing so requires solving a boundary layer problem around the seal as described in appendix B.2 or, in more detail, in sections 1.5–1.6 of the supplementary material. We obtain a regularized version of (12)

$$\gamma h_{0,t}=Q(t)-q,\qquad q=\mathcal{Q}_{s}(\nu^{-1}(h_{0}-b_{m}),b_{x}^{-},b_{x}^{+}).$$

(25)

Note that this replaces the cruder but structurally similar ordinary differential equation models for lake surface lowering in Raymond and Nolan [35], Kingslake et al. [38] and Ancey et al. [39]. The function $\mathcal{Q}_{s}$, is zero when the first argument is non-positive, and has an $O(1)$ derivative with respect to its first argument when the latter is positive. We also assume that the derivative is positive, since we expect a higher water level relative to the seal to lead to a larger flux. This property is confirmed numerically in appendix B.2 and can presumably be proven mathematically, although we have not tried.

Putting $h_{0}=b_{m}+\nu h_{0}^{\prime}$, we recover (12) at leading order, with the flux $q$ implicitly defining the correction $h_{0}^{\prime}$. The latter however represents a steady state solution of a transient that occurs at a faster, $O(\nu)$ time scale: substituting for $\dot{b}_{m}$ from (20) and rescalng $T=\nu^{-1}t$ in (25) yields

$$\gamma\nu h^{\prime}_{0,T}=Q-\gamma w-\frac{\gamma b_{x}^{-}M(-b_{x}^{-},\mathcal{Q}(h_{0}^{\prime},b_{x}^{-},b_{x}^{+}))}{b_{x}^{+}-b_{x}^{-}}-\mathcal{Q}(h_{0}^{\prime},b_{x}^{-},b_{x}^{+})$$

(26)

with $b_{x}^{-}$ and $b_{x}^{+}$ constant on the fast time scale. $h_{0}^{\prime}$ will evolve to a pseudo-steady state that satisfies either (23) or $h_{0}^{\prime}=0$, and which is stable in time $T$ to small perturbations in $h_{0}^{\prime}$. The latter constraint implies that, when there are three solutions to (23), only the largest solution (for which $q$ again increases with base outflow rate $Q-\gamma w(x_{m})$) and the zero solution are viable as indicated by solid lines in figure 5(a). In addition, the relevant, stable solution branch of the original leading-order model (12) is chosen by continuity of $q$ in time whenever such continuity is possible: this is true at least provided there are no significant variations in $Q$ on the short $\sim O(\nu)$ time scale over which the lake level correction $h_{0}^{\prime}$ adjusts. In practice, this could be a real consideration with diurnal water input fluctuations. Presumably these are generally insufficient in practice to lead to $h_{0}^{\prime}$ changing significantly, and do not affect the outflow rate $q$, but a more sophisticated approach is necessary if they do.

For the commonly encountered situation of the upstream side of the lake seal being in steady state, $w(x_{m})=Ub_{x}^{-}$, we can use the stability result to visualize flux $q$ as a multi-valued function of $b_{x}^{+}$ and $b_{x}^{-}$ in the limit of small $\nu$ (figure 6(a)). Here a zero solution $q=0$ is possible everywhere above the red dash-dotted line $b_{x}^{-}=Q/(\gamma U)$, and becomes the only solution in the area demarcated by a solid black curve. Flux $q$ is not continuous across either of those boundaries when transitioning between solution branches. Note also the region bounded to the left by the dashed black contour: here the non-zero flux $q$ is less than the inflow $Q$, with the seal advancing and rising in height, but water still flowing out of the lake.

For $\alpha=0$, the situation is more complicated: we have $M(-p,q)=-pq$ and hence (23) reads

$$q(b_{x}^{+}-b_{x}^{-}-\gamma b_{x}^{-}b_{x}^{+})/(b_{x}^{+}-b_{x}^{-})=Q-\gamma w(x_{m}).$$

(27)

Stability again requires $q$ to increase with $Q-\gamma w(x_{m})$, so the coefficient of $q$ on the left must be positive if $q$ is non-zero, leading to a solveability condition $(b_{x}^{+}-b_{x}^{-}-\gamma b_{x}^{-}b_{x}^{+})<0$ when $Q-\gamma w(x_{m})>0$. By contrast, a unique, stable solution $q=0$ exists when $Q-\gamma w(x_{m})<0$ (figure 5(b)). That, however, potentially leaves a region of the $(b_{x}^{+},b_{x}^{-})$ plane with no viable solution, since we may have $Q-\gamma w(x_{m})>0$ while the solvability condition is violated.

We can again visualize this situation, plotting flux for $\alpha=0$ as a function of $b_{x}^{+}$ and $b_{x}^{-}$ if $w(x_{m})=Ub_{x}^{-}$ (figure 6(b)). Unlike the case $\alpha>0$, $q$ is now not multivalued, but instead undefined in a ‘forbidden’ part of the $(b_{x}^{+},b_{x}^{-})$ plane if $Q>U$ as shown, reflecting the solvability condition.

As a result, breakdown of the model is a very real possibility if $\alpha=0$. If the seal migrates backwards with a steepening upstream slope, $b_{x}^{-}$ can approach the critical value at which the coefficient in (27) vanishes, and $q$ undergoes runaway growth (as the red boundary of the forbidden region is approached from below in figure 6(b)). Physically, incision into the seal becomes fast enough in the forbidden region that the flux $q$ out of the reservoir cannot keep water level close to the seal height (appendix B.2): a finite (but presumably short-lived) jump in water height across the seal will evolve, with very large flux and incision rates that cannot be captured by our reduced model. The lake drains abruptly, and a rescaling of water depth in the channel is required to capture this phenomenon.

We solve (10)–(12) numerically using the method of characteristics with a backward Euler step as described in appendix D. We use the regularized flux prescription (25) for $\alpha>0$, and at times for $\alpha=0$ in order to explore what happens ‘beyond’ the model failure identified at the end of the last subsection. When we do use (25), we treat $\mathcal{Q}_{s}$ simply as a regularization rather than trying to emulate the function shown in figure 17(b). Consequently we drop the slopes $b_{x}^{-}$ and $b_{x}^{+}$ as arguments from $\mathcal{Q}_{s}$. In practice, we use $\mathcal{Q}(h_{0}^{\prime})=\max(h_{0}^{\prime},0)^{2}$, and put $\nu=10^{-3}$.

Figures 7–11 illustrate the behaviour of lakes that are initially empty with $b(x,0)=s(x)$, where $s$ is the unincised ice surface, satisfying $Us_{x}=w$. This initial profile is a steady state solution in the absence of flowing water, and the profile $b$ therefore remains unchanged until the lake is full: only then does water begin to flow and the channel becomes incised on the downstream side of the lake seal. We compute results for an uplift velocity of the form

$$w(x)=U^{-1}\left\{-\overline{b_{x}}-2b_{1}\lambda(x-x_{0})\exp\left[-\lambda\left(x-x_{0}\right)^{2}\right]\right\}$$

with $x_{0}=1.5960$, $b_{1}=\lambda=1$, $\overline{b_{x}}=-0.25$ and $U=1$; this results in a steady state surface $s(x)$ in the form of a Gaussian $b_{1}\exp(-\lambda(x-x_{0})^{2})$ superimposed on a uniform downward slope $-\overline{b_{x}}x$ (see e.g. the top profile in figure 7(a1–a2)). The choice of $x_{0}$ ensures that $Us_{x}=w=0$ at $x=0$, so that the upstream end of the domain is the bottom of the lake. The steady state surface is shown as a dashed black line in figure 8(a2), or as one of the black curves in figure 7, panels a1 and a2. An alternative form of $w$ that asymptotes to a negative limiting value from above for large negative $x$ (meaning, the lake bed does not flatten out above the seal) is explored in section 4.2 the supplementary material.

We use two different choices of shape exponent, $\alpha=0$ (the fixed-width slot of equation (1e)) and $\alpha=1/2$ (the semicircles and triangles of equations (1c) and (1d)). For each of these, we compute solutions for different constant values of water input $Q$ and of storage capacity $\gamma$, treating the latter as independent of water level [see also 42, for a discussion of lake hypsometries]. Both of these assumptions are simplistic, but help understand the dynamics of surface lakes more clearly. Importantly, real ice sheet surface lakes are subject to time-dependent forcing due to seasonal and shorter melt cycles. In many cases, that forcing is quite rapidly varying, since the time scale $[t]$ for ice to traverse the length scale of the seal may be quite long: with $U=100$ m yr^-` and $[x]=$ 1 km, one unit of dimensionless time here corresponds to ten annual melt cycles. We explore the effect of rapidly varying water input in §5.1 of the supplementary material, where we find that it leaves the qualitative behaviour of the system largely unchanged. Here we persist with a constant water input in order to illustrate that qualitative behaviour more cleanly.

Depending on the values of $\alpha$, $Q$ and $\gamma$, different outcomes are possible, differentiated at the coarsest level by whether the lake drains or not. Figure 7 illustrates two possible end members for $\alpha=1/2$. In both cases, outflow from the lake commences once water level (blue curve in panel b) reaches the smooth seal (black line in panel b). For the low-flux example in column 1, with $Q=0.196$, the smooth seal (the maximum in the unincised ice surface at $x=\bar{x}_{m}=1.468$) remains in place, and hence (with $w=0$ at the steady seal), $q=Q$ from the paragraph following equation (22). The channel steepens downstream, but no backward-migrating shock forms. The steepened flowing section terminates in a ponded section that migrates downstream, eventually leaving the entire domain in a new steady state.

Column 2 shows a high-flux counterexample to the steady state of column 1, with $Q=1.570$ and $\gamma=4$. Here a shock forms quickly: the inset in panel b2 shows that the shock (magenta) forms downstream of the smooth seal (black) as predicted at the end of section 3.3, and subsequently migrates upstream to breach the lake. This causes flux $q$ to increase as stored lake water is released. The dowstream side of the shock steepens and a ponded section again forms further downstream. Although the steepening on the downstream side of the seal is eventually reversed (panel a2), the backward migration of the shock continues until the lake is fully drained.

Note that we may naively attribute the steepening of slopes near the smooth seal location $\bar{x}_{m}$ in both columns of figure 7 to melting after outflow from the lake commences. Characteristics offer a different perspective that will be important later: slope evolves as $p_{\tau}=w_{x}$ along characteristics, that is, as the result of differential uplift. Melt enters into the evolution of slope by determining how fast a given characteristic propagates, with $x_{\tau}=U-cM_{-p}$ by (15). Larger fluxes $q$ lead to increased $M_{-p}$ and hence to reduced characteristic velocities. The steepest downward slopes result when $x_{\tau}$ is near zero where the uplift rate derivative $w_{x}$ is most negative (which is indeed near the smooth seal location), causing characteristics to linger. That does not occur at the largest fluxes $q$, however, since $x_{\tau}$ then becomes progressively more negative. The latter effect, of large discharge $q$ flattening slopes, is evident in column 2 of figure 7, where slopes downstream of the shock become less steep as the shock approaches the upstream end of the domain. This flattening of downstream slopes will play a key role in section 4.3 below.

As a further aside, note that we choose not to introduce new characteristics at the bottom end of the domain when the characteristic velocity there is negative, as in figure 7(b2), where the contour plot is blank in the top right-hand corner. Adding characteristics would require an additional boundary condition, the choice of which however should not affect lake drainage unless the characteristics introduced there reach the lake seal.

Depending on storage capacity $\gamma$, more complex behaviour can occur at intermediate values of water supply $Q$ as illustrated in figure 8. The left-hand column shows a higher storage ($\gamma=4$, $Q=0.7850$) example, the right a lower storage ($\gamma=2$, same $Q$) case.

For the former, we see periodic oscillations being generated. As in column 2 of figure 7, a shock forms downstream of the smooth seal, steepening initially and breaching the lake. Lake drainage does not continue to completion, however. The seal stops migrating upstream before reaching the low point of the lake, with slopes downstream of the seal again flattening some time after the lake seal has been breached (panel a1). Outflow $q$ stops and the shock is advected downstream, allowing the initial $q=0$ steady state surface profile to re-establish itself and re-forming a smooth lake seal. The shock that caused the original drainage event migrates some distance downstream before the lake refills and outflow starts afresh (panel b1). Panel d1 shows that a new shock (break in the magenta curve) forms and repeats the drainage of the lake, with the channel profile in the entire domain undergoing periodic oscillations after several cycles of lake drainage and refilling (panels c1, e1).

For the lower storage case of column 2 in figure 8, we again see a shock breaching the seal, partially draining the lake and then stopping, with slopes downslope of the shock initially steepening and then flattening. The shock is again advected downstream and a smooth seal re-forms, but the refilling of the lake occurs more rapidly. The same shock that originally breached the lake is reactivated and breaches the seal again, this time migrating further upstream. The lake fully drains during the third such drainage cycle. We return in section 4.3 to a more detailed analysis of the mechanism by which lake drainage becomes oscillatory, and of the differences between the two cases in figure 8.

The range of drainage styles observed for $\alpha=0$ is more limited. At low water input $Q$, the channel develops into a steady state in much the same way as shown in figure 7, column 1. At larger $Q$, a shock once more forms, although this time at the seal in agreement with sections 3.2–3.3. The outcome of that shock migrating backwards into the lake leads to flux $q$ increasing and one of two outcomes, shown in figure 9.

Column 1 shows a case with more moderate water input $Q=1.1$ and $\gamma=2$. Panel c1 shows results for discharge $q(t)$ and seal position $x_{m}(t)$ from two computations: one without the regularization advocated in equation (25) (magenta and red), and one that is regularized (black and blue). The unregularized model has a singularity in finite time, as expected from the results in section 3.4 (see in particular figure 6(b)): this manifests itself in a very rapid rate of increase $\mathrm{d}q/\mathrm{d}t$ followed by the Newton solver used to compute backward Euler steps failing to find a solution.

The regularized solution instead undergoes very rapid drainage at a slightly later time ($t\approx 4.97$), the timing being different because the regularization in question involves water level in the lake having to rise further to reach the same flux. The singularity in flux in the unregularized model is averted because he regularized model allows water level $h_{0}$ to differ significantly from seal height $b_{m}$: consequently lake drainage can lag behind the rapid lowering of the seal that occurs for $\alpha=0$. That being said, seal incision continues after the very rapid drainage, and lake drainage continues to completion as in column 2 of figure 7.

Column 2 of figure 9, at a larger inflow rate $Q=2$ than column 1, shows a much more straightforward analogue to column 2 of figure 7, with seal incision leading to a peak in flux and continued seal incision until lake drainage is complete, without a (near-) singular peak flux. Importantly, we did not find instances of oscillatory lake drainage for $\alpha=0$.

A more systematic exploration of the effect of changing storage capacity $\gamma$ and water input $Q$ is shown in figure 10 for $\alpha=1/2$ and figure 11 for $\alpha=0$, where we use the unregularized version of the model for the latter. In both cases, inflow rate $Q$ alone determines whether the seal is breached, with the results suggesting that a critical value of $Q$ separates solutions that experience at least partial lake drainage from those that leave the seal intact. The fact that the initial seal breach does not depend on storage capacity $\gamma$ is trivial: until a backward-migrating shock has formed and breached the seal, the intact, steady-state smooth seal leads to outflow balancing inflow once the lake has filled, with $q=Q$, and the shock that breaches the seal must form at that value of flux $q$.

Once the seal is breached, the outcome of lake drainage depends on both $Q$ and $\gamma$. As already indicated above, for $\alpha=0.5$, moderate $Q$ and larger $\gamma$ favour oscillatory drainage of the lake, with smaller $Q$ and larger $\gamma$ ultimately also leading to periodic oscillations rather than divergent oscillations eventually leading to lake drainage. For $\alpha=0$, smaller $Q$ and larger $\gamma$ by contrast favour blow-up of the solution with singular outflow $q$; oscillatory lake drainage does not occur in any of the computations we have performed.

For constant water input to the lake $Q$ with $b$ in steady state upstream of the seal, there appears to be a critical value of $Q$ above which a shock either forms at the seal (if $\alpha=0$) or below the seal (if $1>\alpha>0$). The shock migrates backwards, leading to at least partial lake drainage.

Below, we identify breaches of the seal with parameter combinations for which there is no steady state solution to (10)–(12), since steady states are naturally the solution for large times $t$ if the lake seal is not breached. To see that steady states are a natural consequence of the seal remaining intact, note that the upstream end of the domain has constant $b(0,t)=b_{in}(0)$, and hence $b_{t}=-\mathcal{H}=0$ there. If the seal is not breached and $q$ is therefore constant, then the value of the Hamiltonian remains constant along any characteristic that does not cross the upstream end of a ponded section. Specifically, denote by $\tilde{\mathcal{H}}(\sigma,\tau)$ the value of the Hamiltonian as a function of the characteristic variables, so

$$\tilde{\mathcal{H}}_{\tau}=\mathcal{H}_{x}x_{\tau}+\mathcal{H}_{p}p_{\tau}+\mathcal{H}_{q}q_{\tau}=\mathcal{H}_{x}\mathcal{H}_{p}-\mathcal{H}_{p}\mathcal{H}_{x}+\mathcal{H}_{q}q_{\tau}=\mathcal{H}_{q}q_{\tau}$$

(28)

and $\tilde{\mathcal{H}}$ is constant if $q$ is. Note that equation (28) holds except possibly where a $cM$ changes discontinuously along a characteristic, which is possible only at the upstream end of a ponded section (see appendix C).

From (28) and $\mathcal{H}=0$ at $x=0$, it follows that any point on a characteristic that originates at the upstream end of the domain is also in a local steady state with $b_{t}=-\tilde{\mathcal{H}}=0$, provided the seal remains steady and $q$ is therefore constant. This behaviour is evident in column 1 of figure 7, where steady state conditions are established progressively down-flow as characteristics that cross the seal after lake outflow has commenced propagate downstream across the domain, with the ponded region progressively moving down-flow as well.