The drainage of glacier and ice sheet surface lakes

Abstract

Supraglacial lakes play a central role in storing melt water, enhancing surface melt and ultimately in driving ice flow and ice shelf melt through injecting water into the subglacial environment and through facilitating fracturing. Here, we develop a model for the drainage of supraglacial lakes through the dissipation-driven incision of a surface channel. The model consists of the St Venant equations for flow in the channel, fed by an upstream lake reservoir, coupled with an equation for the evolution of channel elevation due to advection, uplift and downward melting. After reduction to a ‘stream power’-type hyperbolic model, we show that lake drainage occurs above a critical rate of water supply to the lake due to the backward migration of a shock that incises the lake seal. The critical water supply rate depends on advection velocity and uplift (or more precisely, drawdown downstream of the lake) as well as model parameters such as channel wall roughness and the parameters defining the relationship between channel cross-section and wetted perimeter. Once lake drainage does occur, it can either continue until the lake is empty, or terminate early, leading to oscillatory cycles of lake filling and draining, with the latter favoured by large lake volumes and relatively small water supply rates.

1. Introduction

Large areas of the Greenland ice sheet experience surface melt water drainage (Poinar & Andrews 2021), while surface drainage is confined to lower elevations in Antarctica (Lenaerts et al. 2016; Kingslake et al. 2017; Stokes et al. 2019). Surface melt drives the evolution of subglacial drainage systems in Greenland (Das et al. 2008; Chandler et al. 2013; Cowton et al. 2013), which in turn control sliding speed (Shepherd et al. 2009; Schoof 2010; Palmer et al. 2011; Tedesco et al. 2013). Accumulation of surface water can also enhance surface melting by reducing albedo (Lüthje et al. 2006; Tedesco et al. 2012) and cause the break-up of floating ice shelves (Scambos et al. 2004,  2009; Banwell, MacAyeal & Sergienko 2013; Lai et al. 2021), while the injection of surface melt under a floating ice tongue can drive convection in fjords (Straneo et al. 2011; Mortensen et al. 2020), and enhance as well as localize melting at the base of the ice tongue (Dallaston, Hewitt & Wells 2015; Washam et al. 2019).

Lakes situated in local depressions on the ice surface are common features of drainage systems in both Greenland and Antarctica. These lakes store water, enhance surface melt and, in the case of Greenland, can cause short-lived acceleration of ice flow through abrupt drainage to the bed by hydrofracturing (van der Veen 2007; Shepherd et al. 2009; Stevens et al. 2015; Christoffersen et al. 2018). Much research has focused on the latter effect, even though a significant fraction of surface lakes in Greenland either drain slowly or not at all (Koenig et al. 2015; Lampkin et al. 2020; Law et al. 2020; Benedek & Willis 2021; Dunmire et al. 2021; Poinar & Andrews 2021), while there are no known surface lakes on the grounded part of the Antarctic ice sheet that drain to the bed (Bell et al. 2018).

Motivated by field observations made in Antarctica, we consider the case of lakes draining purely through channels incised into the surface of a grounded ice sheet (as opposed to a floating ice shelf). On a grounded ice sheet, the surface depression occupied by a lake is usually generated by ice flow over suitably uneven bed topography under the ice sheet, and lakes are often observed to occupy the same position for long periods of time (that is, over many summer melt seasons). However, a combination of remote sensing imagery and ground-based radar (Schaap et al. 2020) suggests that such surface lakes can also drain relatively quickly through near-surface channels, and can do so after a lengthy periods of apparently steady lake levels. The same observations also suggest that drainage can occur in winter, when there is presumably little to no water input over winter. The question that arises is: What controls lake drainage through such a channel? Similar overland drainage may also be relevant to higher elevations in Greenland (Benedek & Willis 2021), where few lakes drain through hydrofracture (Poinar & Andrews 2021). However, we neglect seepage into a firn aquifer in our work (Forster et al. 2013; Meyer & Hewitt 2017), which may be relevant for some of these Greenlandic lakes.

There have only been a handful of attempts to model lake drainage along glacier and ice sheet surfaces through thermal erosion of a channel through an ice dam (Walder & Costa 1996; Raymond & Nolan 2000; Mayer & Schuler 2005; Vincent, Auclair & Le Meur 2010; Kingslake, Ng & Sole 2015; Ancey et al. 2019). Most of these consider drainage along more steeply angled glacier surfaces, where flow is likely to be Froude supercritical. In all of these previous studies except Kingslake et al. (2015), surface lakes are considered as natural hazards, with the ultimate aim of computing hydrographs for rapid surface drainage. In addition, and in contrast with models for drainage along the glacier bed (Nye 1976; Spring & Hutter 1981; Clarke 1982; Ng 2000; Kingslake & Ng 2013; Stubblefield et al. 2019; Schoof 2020) none of the surface drainage models listed resolve channel incision (and therefore channel slope) as a function of position along the flow path, but instead take the form of ‘lumped’ models intended to describe conditions near the channel intake only.

Although the model we develop is in principle applicable to the outburst floods studied previously, our main goal differs substantially from these prior studies. We are interested primarily in whether water input to a lake, causing the lake to overflow, necessarily leads to lake drainage by channel incision, and whether drainage can be partial or must continue until the lake basin is completely empty. In a lake with a water supply, it is natural to assume that channel incision on its own should drain the lake completely, since the water input should ensure there is discharge in the channel and therefore continued erosion of the ice dam. On a moving ice sheet, however, advection will also carry the channel downstream and can potentially re-build the seal of the lake.

As a result, we focus on systems in which incision of the channel is quite slow, and competes with horizontal advection and vertical uplift of the ice surface due to the flow of the ice sheet over bed topography. As pointed out, these latter two processes are responsible for shaping the surface depression occupied by the lake to begin with (Schoof 2002), but they have rarely been considered in the context of surface lake dynamics (Darnell et al. 2013) and are not incorporated in the existing models for rapid lake drainage.

Among other consequences, the incorporation of advection forces us to employ a partial differential equation-based model, resolving position along the channel as well as time. The model we derive bears close resemblance to so-called stream power models for fluvial landscape evolution in non-glacial contexts (Luke 1972). The latter typically incorporate uplift (e.g. Whipple & Tucker 1999; Royden & Perron 2007; Kwang & Parker 2017), but the additional effect of horizontal advection is not commonly considered as part of fluvial landscape evolution.

The paper is organized as follows: in § 2.1 we define a basic model consisting of the St Venant equations for a surface stream coupled with an evolution equation for channel depth, based on local dissipation driving channel incision. In §§ 2.2–2.4, the model is reduced based on a small ratio of water depth in the channel to lake depth, and water velocities being much larger than ice velocities, while the local Froude number is assumed to remain subcritical. This results in a nonlinear hyperbolic evolution equation for channel evolution, coupled to an evolution equation for lake volume (§ 3.1). The formation of shocks in the model and how they control discharge from the lake is studied in §§ 3.2–3.5, with boundary layer solutions of the full model around the shocks relegated to Appendices A–B. Numerical solutions by the method of characteristics (Appendix E) are given in § 4, where we show that lake drainage occurs above a critical value of water supply to the lake (§ 5), which can result either in complete lake drainage, or in oscillatory cycles of lake drainage and refilling (§ 6).

2. Model

2.1. Model formulation

We consider a surface melt water stream with cross-sectional area $S$, with the base of the stream channel at an elevation $b$, and denote the mean velocity in a given cross-section of the stream by $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. Fowler 2011, Chapter 4)

(2.1a)\begin{gather} \rho_{w} \left[ S_t + (uS)_x \right] = \rho_i m, \end{gather}

(2.1b)\begin{gather}\rho_{w} S (u_t + u u_x) = - \frac{\rho_{w} f u^2P(S)}{8} - \rho_{w} g S \left[ b_x + h(S)_x \right], \end{gather}

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; g is acceleration due to gravity; $\rho _{i}$ and $\rho _{w}$ are the densities of ice and water, respectively, 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.

Figure 1. Geometry of the problem: water surface in blue, channel/lake bottom in black, ice surface as dashed black line. Some of the symbols used here ($b_{m}$, $x_{m}$ $x_{s}$ and $x_{p}$) are defined in the context of a leading-order model in §§ 2.2–3.

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{\rm \pi} ^{-1}S)^{1/2}$ and

(2.1c)\begin{equation} P(S) ={\rm \pi} r = (2 {\rm \pi}S)^{1/2}, \quad h(S) = r = (2 {\rm \pi}^{ - 1} S)^{1/2}. \end{equation}

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

(2.1d)\begin{equation} P(S) = \frac{2S^{1/2}}{\sin(2\theta)^{1/2}}, \quad h(S) =\frac{S^{1/2}\cos(\theta)}{\sin(2\theta)^{1/2}} \end{equation}

or to fix a width $W$ that is much greater than water depth, and to put (as is done in Fowler 2011)

(2.1e)\begin{equation} P(S) = W, \quad h(S) = \frac{S}{W}. \end{equation}

Generically, this suggests we consider

(2.1f)\begin{equation} P(S) = c_1 S^\alpha, \quad h(S) = c_2 S^\beta, \end{equation}

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): (2.1c)–(2.1d) have $\alpha = \beta = 1/2$ while (2.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$.

Figure 2. Cross-section shapes: (a) semi-circle ($\alpha = \beta = 1/2$), (b) triangular ($\alpha = \beta = 1/2$) and (c) fixed-width slot ($\alpha = 0$, $\beta = 1$). Water with cross-sectional area $S$ is shown in blue, wetted perimeter in heavy black. The qualitative nature of solution computed in § 4 depends on whether $\alpha = 0$.

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 Evatt et al. 2006; Ogier et al. 2021); 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

(2.1g)\begin{equation} \rho_{i} \mathcal{L} m = \frac{\rho_{w} f u^3 P(S)}{8}, \end{equation}

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 f u^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 (Dallaston & Hewitt 2014) 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. Schoof 2002), 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

(2.1h)\begin{equation} b_t + Ub_x = w -\frac{m}{P(S)}. \end{equation}

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 (Jarosch & Gudmundsson 2012) as well as snow accumulation at the base of the channel during winter, but also that $s$ is in steady state (Schoof 2002), and itself satisfies

(2.1i)\begin{equation} Us_x = w. \end{equation}

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

(2.1j)\begin{equation} \dot{V} = q_0(t) - \left.(uS)\right|_{x = 0}, \quad V(t) = V_{L}(h(S(0,t))), \end{equation}

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. Also, $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.

Note that (2.1j) places the upstream boundary of the model at a fixed position $x = 0$ rather than the moving seal location $x = x_m$ indicated by figure 1. The latter would certainly be appropriate, but a fixed upstream domain boundary $x = 0$ causes no inconsistencies here: we will find shortly that melt rate $m = 0$ upstream of the lake seal at leading order, so that the channel bed elevation $b$ predicted by (2.1h) simply follows the unincised ice surface $s$ up to the seal, and in addition, water flux $uS$ is independent of position at leading order, so $(uS)|_{x = x_m} = (uS)|_{x=0}$. Extending the domain upstream of $x_m$ allows for a simpler presentation of the physics of seal migration in the leading-order version of the model that we will derive next, since those physics are simply those of a shock (or slope discontinuity) that can equally form further downstream in the channel.

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. Karlstrom, Gajjar & Manga 2013; Fernández & Parker 2021). 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 § 7.

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$, but 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 Buzzard, Feltham & Flocco 2018). 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 Raymond & Nolan 2000). 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.

2.2. Non-dimensionalization and a reduced model

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]$ and $[b]$ through

(2.2a–c)\begin{equation} \frac{\rho_{w} g [b] [S]}{[x]} = \frac{\rho_{w} f [u]^2 P([S])}{8} , \quad \frac{ \rho_{w} f [u]^3}{8\rho_{i} \mathcal{L}} = [w] = \frac{[U][b]}{[x]}, \quad [U][t] = [x]. \end{equation}

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 scale $[b] = [w][x]/[U]$, the advective time scale sets a dissipation rate and therefore a water flux scale $[u][S]$ (effectively, as a distinguished limit). Setting the dissipation rate based on the advective time scale may seem contrived, since it should really be set by surface slope and water flow rates. The latter are controlled by water supple, and are not physically controlled by surface advection. We assume that we are in a parameter regime in which water supply produces a melt-driven incision rate that is comparable to uplift and advection. The alternative would be a much faster melt rate (which the model we construct can still capture if we prescribe a large dimensionless water supply). In that case, however, lakes generally cannot persist as the uplift that is necessary to create a lake seal cannot compete with the incision rate of the channel.

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

(2.3a–d)\begin{equation} \nu = \frac{h([S])}{[b]}, \quad Fr^2 = \frac{[u]^2}{g h([S])}, \quad \varepsilon = \frac{g[b]}{\mathcal{L}}, \quad \delta = \frac{[U]}{[u]}. \end{equation}

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$^{-2}$ 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 (2.1c)

(2.4a,b)\begin{equation} [u] = 1.2 \, \text{m}\,\text{s}^{ - 1}, \quad [S] = 0.47 \, \text{m}^2, \end{equation}

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

(2.5a–e)\begin{equation} x = [x]x^*, \quad t = [t]t^*, \quad u = [u]u^*, \quad S = [S]S^*, \quad b = [b]b^*, \end{equation}

and define

(2.6a,b)\begin{equation} P^*(S^*) = \frac{P(S)}{P([S])}, \quad h^*(S^*) = \frac{h(S)}{h([S])}, \end{equation}

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

(2.7a)\begin{gather} \delta S_t + (uS)_x = \varepsilon u^3 P(S), \end{gather}

(2.7b)\begin{gather}\nu Fr^2 S (\delta u_t + uu_x) = - u^2 P(S) - S b_x - \nu S h(S)_x, \end{gather}

(2.7c)\begin{gather}b_t + U b_x = w - u^3, \end{gather}

with $P(S) = S^\alpha$ and $h(S) = S^\beta$.

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 (2.7)

(2.8a–c)\begin{equation} (uS)_x = 0, \quad u^2 P(S) = - Sb_x, \quad b_t + Ub_x = w - u^3. \end{equation}

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 available at https://doi.org/10.1017/jfm.2023.130, also §§ 4.4.4–4.5.2 and chapter 5 of Fowler 2011). 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 (Holmes 1995). We leave this to future work.

Persisting with (2.8a–c), we find that $q = uS$ is independent of position, and we will assume below that $q > 0$, so the lake at the upstream end of the domain drains through the channel, but is not filled through a reverse flow. With fixed $q$, $u$ depends on flux $q$ and slope $-b_x$ through (2.8a–c)$_2$ as

(2.9)\begin{equation} u^3 q^{ - 1} P(u^{ - 1}q) = - b_x, \end{equation}

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

(2.10)\begin{equation} M( - b_x,q) := u^3 = \left( - q^{1-\alpha}b_x\right)^{3/(3-\alpha)}. \end{equation}

The function $M$ here is a 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.

Figure 3. Melt rate $M(-b_x,q)$ against $b_x$ for $q = 0$, $0.25$, $0.5$, $0.75$, $1$ for (a) $\alpha = 0.5$, (b) $\alpha = 0$; $M = 0$ when $b_x > 0$. Note that, although the two panels look similar, $M$ in panel (a) is strictly convex for $b_x < 0$ and continuously differentiable at $b_x = 0$, which has significant implications for shock formation in the model (2.11).

At face value, substituting in (2.8a–c)$_3$ yields the single evolution equation for $b$,

(2.11)\begin{equation} b_t = w - Ub_x - M( - b_x,q). \end{equation}

Note that this is effectively the stream power equation for landscape evolution in Luke (1972), Fowler, Kopteva & Oakley (2007) and Kwang & Parker (2017).

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 therefore more generally put $M = 0$ when $b_x > 0$, replacing (2.10) by

(2.12)\begin{equation} M( - b_x,q) = q^{3(1-\alpha)/(3-\alpha)}\left[\max\left( - b_x,0\right)\right]^{3/(3-\alpha)}, \end{equation}

to account for this. That in itself does not, however, suffice, since local maxima in the stream bed can induce ponding even on downward slopes further upstream. We deal with this next.

2.3. Ponding

Ponding occurs at a point $x$ when there is a downstream point $x' > x$ at which the base of the channel is higher than at $x$, $b(x') > b(x)$. The appropriate modification of (2.11) to account for ponding is therefore via a ‘ponding function’ $c$

(2.13a)\begin{gather} b_t = w - Ub_x - c(x,t) M( - b_x,q), \end{gather}

(2.13b)\begin{gather}c(x,t) = \left\{\begin{array}{ll} 1, & \text{if } b(x,t) \geq \sup_{x' >x} b(x',t), \\ 0, & \text{otherwise}. \end{array} \right. \end{gather}

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. It is also worth pointing out that, if there is no flow ($q = 0$), the ‘ponded’ sections of the bed (given by the set $\{ x: b(x,t) < \sup _{x'>x} b(x',t) \}$) may not be fully submerged by stagnant water, but this does not alter the evolution problem (2.13a) further since the absence of flow already ensures that $M = 0$.

We assume additionally that ponded sections of channel store negligible quantities of water, so that we can continue to treat $q$ as independent of position $x$. Formally, we can use a rescaling as described in Appendix A to show that negligible storage corresponds to the parameter regime $\delta \ll \nu ^{1/\beta }$. In that case, the lake generally stores much more water than the ponded sections, since drainage of the lake affects flux $q$, while drainage of a ponded section does not. Physically, this occurs because the lake is much wider than the channel, as it occupies a depression in the unincised ice surface $s$ (figure 1). Since we assume $s$ to be in steady state, the lake basin shape is unaffected by the evolution of $b$, although the water level within that basin does depend on channel evolution as we describe immediately below in § 2.4. To make the model self-consistent, we also avoid the possibility of multiple such lakes by insisting that the uplift function $w$ have a single root at some location $\bar {x}_m$, with $w < 0$ downstream of that. The only depression in the unincised ice surface given by $Us_x = w$ is then upstream of $\bar {x}_m$.

We still need to deal with mass conservation equation (2.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 (2.1j) in a leading-order, dimensionless form yet. We do so next.

2.4. Outflow at the lake

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$ is contained in a depression in the unincised ice surface, with the water level in the lake controlled by ponding at the upstream end of the channel (that is, by the highest point in the channel bed). Water level in the lake therefore scales with ice surface topography $[b]$. To account for this, define a dimensionless water depth scaled with $[b]$ as

(2.14)\begin{equation} \hat{h}^* = \nu^{ - 1}h^*(S^*) = h(S)/[b], \end{equation}

as is appropriate for ponded sections, see Appendix A; note that this differs from the scaling for water depth in channel further downstream. Then

(2.15)\begin{equation} h_0^* = \hat{h}^*(0,t^*) + b^*(0,t^*) \end{equation}

is the dimensionless water level of the lake, relative to the same datum as channel bottom elevation $b^*$. We define a dimensionless lake volume function and a dimensionless water supply function through

(2.16a,b)\begin{equation} \hat{V}(h_0(t^*)) = \frac{V_L(h(S(0,t)))}{[u][S][t]}, \quad Q^*(t^*) = \frac{q_0(t)}{[u][S]}, \end{equation}

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 to 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 downstream from the domain boundary at $x = 0$, the latter being upstream of the lake seal unless the lake drains completely). Water surface elevation also cannot exceed seal height at leading order (Appendices A and B.2). Consequently, we find that water level 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, in which case no water is flowing. We denote the seal height by $b_{m}(t)$, so that

(2.17a)\begin{equation} h_0 \leq b_{m}(t) := \sup_{x > 0} b(x,t). \end{equation}

Similarly, we will use $x_{m}$ to denote the seal location, defined such that $b_{m}(t) = b(x_{m}(t),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

(2.17b)\begin{gather} \gamma \dot{h}_0 = Q(t) - q, \end{gather}

(2.17c)\begin{gather}q = \left\{ \begin{array}{ll} 0, & \text{if } h_0 < b_{m}, \\ \max\left( Q - \gamma \dot{b}_{m},0 \right), & \text{if } h_0 = b_{m}, \end{array} \right. \end{gather}

where $\gamma (h_0) = \mathrm {d} \hat {V} / \mathrm {d} h_0$ is storage capacity in the lake (given by its surface area), 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.

3. Characteristics, shocks and the dynamics of the lake seal

3.1. Characteristics

If we treat $c(x,t)$ and $q(t)$ momentarily as known, then (2.13a) can be recognized as being of Hamilton–Jacobi form (Luke 1972),

(3.1)\begin{equation} b_t = - \mathcal{H}(x,t,b_x,q), \end{equation}

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

(3.2)\begin{equation} \mathcal{H}(x,t,p,q) = Up + c(x,t)M( - p,q) - w(x), \end{equation}

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 (Courant & Hilbert 1989, § 3) allows us to write the problem (3.1) in the form of Charpit's equations as follows: we define characteristics as curves of constant $\sigma$ in the transformation $(\sigma,\tau ) \mapsto (x,t)$ given by

(3.3)\begin{equation} x_\tau = \mathcal{H}_p = U - c M_{ - p}( - p,q), \quad t(\sigma,\tau) = \tau, \end{equation}

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. Equation (3.3) underlines a key difference from classical stream power models: here, characteristics can travel downstream as well as upstream, with major implications for breaching the lake seal and controlling flux $q$.

Along a given characteristic, $b(\sigma,\tau )$ and $p(\sigma,\tau ) = b_x$ evolve as

(3.4)\begin{equation} b_\tau = - \mathcal{H} + \mathcal{H}_p p = w - c [M_{ - p}( - p,q)p + M( - p,q)] , \quad p_\tau = - \mathcal{H}_x = w_x, \end{equation}

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. Prescribed $b$ at the upstream end of the domain is appropriate for characteristics entering the domain there (as is always the case when there is a ponded section at that upstream boundary, in which case the characteristic velocity $x_\tau = U$ there).

There are also situations in which the characteristic velocity $x_\tau$ can become negative at the downstream end of a fixed domain, requiring additional boundary conditions there. In practice, surface channels either terminate abruptly at near-vertical cracks (or moulins) in the ice, or at the downstream margin of the ice sheet. Neither situation is adequately described by our model, and we use the following, somewhat unsatisfactory device instead: we fix a downstream domain boundary at some $x = L$ with suitably large $L$, and truncate any characteristic that reaches that location from upstream. Conversely, if the characteristic velocity at $x = L$ becomes negative, we do not introduce new characteristics at $x = L$ but allow the domain to shrink at the characteristic velocity. Implicitly, we are assuming that none of the ‘missing’ characteristics are able to reach the lake seal, and that the omitted physics at the downstream end of the domain does not change the ponding function $c$ by creating a local maximum in $b$ somewhere downstream.

As already pointed out, the problem (3.1)–(3.2) is a modification of ‘stream power models’ for landscape evolution (Luke 1972). The main differences are the control of water flux $q$ through lake drainage rather than simply a prescribed precipitation rate, the inclusion of the ponding function $c$, and the presence of the advection term $Up$. The latter leads to a convex but non-monotone Hamiltonian, which allows characteristics to propagate downstream as well as upstream. As we will see in § 3.3, this feature of the model is key to understanding whether the seal of a lake is breached by incision of the channel or not.

Before we proceed further, there are a couple of technical points to make. First, note that we do not attempt to differentiate the piecewise constant function $c$ when forming the derivative $\mathcal {H}_x$ in (3.4). Instead, we treat discontinuities in $c$ separately as described in detail in § 3.3 and Appendix D. Of these, only the discontinuity at the upstream end of a ponded section described in Appendix D poses any real difficulties: near the downstream end of a ponded section (§ 3.3), $c$ is obsolete because $M$ and $c$ both vanish when bed slope $b_x$ is positive.

Formally, the inclusion of $c$ actually turns the model into a hyperbolic system not only in terms of $b$, but also of a second variable $\breve {b}(x,t) = \sup _{x'>x} b(x',t)$. Assuming that $b$ has an integrable derivative, $\breve {b}$ satisfies (we are grateful to one of the referees for pointing this out)

(3.5)\begin{equation} \breve{b}_x = \min(b_x,0)H(\breve{b}-b), \end{equation}

where $H$ is the usual Heaviside function (with $H(0) = 1$), with $\breve {b} = b$ at the downstream end of the domain. The characteristics of $\breve {b}$ are then lines of constant $t$, and in terms of $b$ and $\tilde {b}$,

(3.6)\begin{equation} c = H(b-\breve{b}). \end{equation}

The auxiliary variable $\breve {b}$, however, only affects the solution through discontinuities in $c$, and on either side of such a discontinuity, (3.1) and hence (3.3)–(3.4) hold with constant $c=0$ or $c = 1$.

Second, it is worth pointing out that the Hamiltonian structure of the problem furnishes a simple evolution equation of $\mathcal {H}$ along characteristics

(3.7)\begin{equation} \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; \end{equation}

the only situation where this fails is when $cM$ and hence $\mathcal {H}$ change discontinuously.

With our choice of constant upstream boundary conditions (so $b_t = -\mathcal {H} = 0$ at the upstream end of the domain), an important corollary of (3.7) is that $b$ always remains in steady state upstream of the seal of the lake, since $\mathcal {H}_q = c M_q = 0$ in a ponded section with $c = 0$, and hence $b_t = -\mathcal {H}$ remains zero along characteristics entering the domain from upstream, at least before they reach the seal of the lake. As we will see later, the ability of such characteristics to fill the entire domain ultimately determines much of the dynamics of the system.

Next, we give a comprehensive account of shocks (that is, discontinuities in slope $b_x$ along which characteristics intersect), and of discontinuities in $c$ (which need not correspond to shocks). Figure 4 provides an overview of the different possibilities. We treat cases (a–c) in the figure in §§ 3.2–3.3 and derive formulae for flux $q$ in terms of shock geometry at the lake seal in § 3.4. We relegate the analysis of the upstream end of ponded sections as shown in figure 4(d) to Appendix D, where we show that the discontinuity in $c$ at such a location cannot generate a shock but may give rise to an expansion fan. The material below is fairly dense and, at this stage, abstract. On a first reading, it may be preferable to skip to § 4 to understand the zoology of features of the solution before filling in the theoretical background in §§ 3.2–3.4 and Appendix D.

Figure 4. Different flavours of shocks and discontinuities in $c$: (a) ‘knickpoint’ shocks in flowing sections (§ 3.2), (b) seal shocks (§ 3.3), (c) smooth seals (§ 3.3) and (d) upstream ends of ponded sections, which can correspond to expansion fans but not to shocks (Appendix D).

3.2. Knickpoint shocks

Equation (3.1) breaks down when characteristics intersect at shocks. Intersections require characteristics to travel faster upstream of the shock 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 4a). These shocks represent ‘knickpoints’ in standard geomorphological parlance (Luke 1972; Royden & Perron 2007).

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. Solving for shock velocity $\dot {x}_c$, using (3.1) to substitute for $b_t^-$ and $b_t^+$, we obtain (see also Royden & Perron 2007)

(3.8)\begin{equation} \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^-}, \end{equation}

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 (2.12), so $\mathcal {H}$ is 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.

3.3. The downstream end of a ponded section

Shocks between a ponded section upstream and a flowing section downstream ($c^- = 0$, $c^+ = 1$, $b_x^- > 0$, $b_x^+ < 0$, see figure 4b) are equivalent to knickpoint-type shocks. Equation (3.8) still holds, where now $c^- M(-b_x^-,q) = 0$ (as pointed out before, the discontinuity in $c$ is a red herring here since $M(-b_x^-,q) = 0$ for $b_x^- > 0$ anyway). The important distinction with the knickpoint shocks of § 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 knickpoint between ponded and flowing, and its motion then controls the flux $q$ as described in § 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 $c M$ and Hamiltonian $\mathcal {H}$ (figure 4c). We refer to this as a smooth seal. The simplest situation in which to understand this is that of a steady state upstream of the seal (which we generally expect to be the case for the lake seal as discussed after (3.7)), in which case a smooth seal will also be stationary. Characteristics will then pass through the seal provided the characteristic velocity downstream remains positive

(3.9)\begin{equation} x_\tau^+ = U - M_{ - p}(0^-,q) > 0. \end{equation}

For $\alpha > 0$, $M_{-p}(0^-,q) = 0$, and (3.9) is always satisfied: in order to breach the lake seal in that case, a knickpoint must form downstream and migrate to the seal location. Conversely, for $\alpha = 0$, such knickpoints cannot form, but (3.9) is violated if $q > U$. In that case, a shock forms at the seal itself, causing the lake seal to be breached.

Note that the argument above can be generalized to the case of non-steady smooth seals by differentiating both sides of $b_x^+(x_s(t),t) = b_x^-(x_s(t),t)$ with respect to time, and using the fact that $b_{xx} < 0$ on either side of $x_s(t)$. Details are provided in Appendix C.

3.4. The lake drainage flux $q$

Key to the model is to understand how the flux $q$ evolves, which requires the evolution of the seal point height $\dot {b}_{m}$ in (2.17c). There are two scenarios. First, the seal point $x_{m}$ can be a shock (note that this differs from figure 1, which shows a smooth seal). Differentiate $b_m = b(x_m(t),t)$ and use the fact that the shock velocity $\dot {x}_m$ is given by (3.8) with $c^-M(-b_x^-,q) = 0$, while $b_t^- = w(x_m) - Ub_x^-$. We obtain

(3.10)\begin{equation} \dot{b}_{m} = b^-_t + b_x^- \dot{x}_{m} = w(x_{m}) + \frac{b_x^-M( - b_x^+,q)}{b_x^+ - b_x^-}. \end{equation}

Assuming lake level is equal to seal height with $h_0 = b_{m}$, (2.17c) then leads to the equation

(3.11)\begin{equation} q = \max\left(Q - \gamma w(x_{m}) - \gamma \frac{b_x^-M( - b_x^+,q)}{b_x^+ - b_x^-},0\right). \end{equation}

The flux $q$ is the sum of water supply to the lake $Q$, and of water discharged due to lowering of the lake seal by uplift $-\gamma w(x_m)$ and melt-driven incision into the lake $-\gamma M/(b_x^+-b_x^-)$. If that sum is negative, there is no outflow $q$ from the lake, and the seal of the lake will in fact temporarily rise above the level of the lake. Importantly, (3.11) determines the flux implicitly, since $q$ appears on the right-hand side as a result of the melt rate being dependent on flux, and that melt rate in turn dictates the rate at which the lake seal is lowered.

Alternatively, the seal can be a smooth transition point, with $b_x^- = b_x^+ = 0$. Then

(3.12)\begin{equation} \dot{b}_{m} = b_t + b_x \dot{x}_{m} = b_t + U b_x = w(x_{m}). \end{equation}

In that case, (2.17c) leads to the explicit formula

(3.13)\begin{equation} q = \max(Q-\gamma w(x_{m}),0). \end{equation}

In what follows, we refer to $Q - \gamma w(x_{m})$ as the ‘base outflow rate’ that results if there is no incision into the seal due to melting (that is, if we simply put $M = 0$ at the seal in (3.11)). A negative base outflow rate signifies that uplift at the seal occurs faster than the refilling of the lake, and must be compensated by a positive incision rate of the seal in order for any outflow to occur at all.

Only the case (3.11) of a shock at the seal is non-trivial, precisely because $q$ appears on both sides of the equation. Since $M \geq 0$ for all $q$ and $M = 0$ if $q = 0$, we can equivalently write

(3.14)\begin{equation} \left. \begin{array}{ll@{}} Q - \gamma w(x_{m}) = q + \gamma b_x^- M( - b_x^+,q)/(b_x^+ - b_x^-), & \text{if } q > 0,\\ Q - \gamma w(x_{m}) \leq 0, & \text{if } q = 0. \end{array} \right\} \end{equation}

At this point, we have to distinguish between the cases $\alpha = 0$ and $\alpha > 0$, which give qualitatively different results.

3.5. Flux for the variable channel width case $\alpha > 0$

For $\alpha > 0$ and $p < 0$, the function $M(-p,q)$ defined in (2.12) is an increasing, concave function of $q$ with $M_q(-p,0) = \infty$. The right-hand side of (3.14)$_1$ therefore vanishes at $q = 0$, decreases for small $q$, reaches a global minimum and then increases thereafter (recall that $b_x^- > 0$ and $b_x^+ < 0$ for a shock at the seal). The right-hand side of (3.14)$_1$ is shown as a function of $q$ in figure 5(a) for fixed $b_x^+$ and a variety of values of $q_x^+$ as dashed and solid curves. The solution for $q$ can be read off the graph by identifying where the height of the curve reaches the prescribed value of base outflow rate $Q - \gamma w(x_m)$.

Figure 5. Values of base outflow rate $Q-\gamma w(x_{m})$ corresponding to a given flux $q$ as determined by (3.11), for $\gamma = 2$, $b_x^+ = -1$ and $b_x^- = 0$, $0.25$, $0.5$ …, $2$ (the end member cases $b_x^- = 0$ and $b_x^- = 2$ being labelled with arrows) for $\alpha = 0.5$ (a) and $\alpha = 0$ (b). Stable solutions are shown as solid lines, unstable as dashed lines. In panel (a), stable $q$ can be multi-valued for given $Q-\gamma w(x_{m})$, while in panel (b), there are combinations of $Q-\gamma w(x_{m})$ and $b_x^-$ for which no solution for $q$ exists (take for instance $Q-\gamma w(x_{m}) > 0$ and $b_x^- = 0.2$).

If the base outflow rate $Q - \gamma w(x_{m})$ is positive, there is therefore a single positive root for $q$ (a single value of $q$ for which the curve attains the prescribed value of $Q - \gamma w(x_m)$). That root $q$ increases with base outflow rate $Q-\gamma w(x_{m})$ and with storage capacity $\gamma$ (for fixed base outflow rate). By contrast, the solution for $q$ can become multi-valued for negative $Q - \gamma w(x_{m}) \leq 0$: no flow with $q = 0$ is then a valid solution of the original problem (3.11), since no flow implies the absence of dissipation-driven seal incision, and seal height will increase at or above the rate of lake filling. That does not mean that there cannot be any flow, however. In addition, there are two non-zero solutions if the negative base outflow $Q - \gamma w(x_{m}) < 0$ remains above a critical value (figure 5a). Incision of the seal by flowing water can then cause drainage of the lake at the right rate to maintain that rate of incision.

For the melt rate $M$ given by (2.12), that situation is possible when

(3.15)\begin{equation} 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)}, \end{equation}

where the critical value is the minimum of the right-hand side of (3.14)$_2$ with respect to $q$. For even more negative $Q - \gamma w(x_{m})$, $q = 0$ is the only solution.

A multi-valued solution for flux $q$ begs the question how outflow from the lake should be computed: How does the lake choose which solution branch to follow? In common with similar situations such as a glacier selecting a branch of a multi-valued sliding law during a surge cycle (Fowler 1987,  1989), or the behaviour of the van der Pol oscillator in the relaxation oscillation limit (Holmes 1995), we argue that multi-valuedness is simply the result of a process that occurs on a fast time scale having evolved to equilibrium, and the dynamics of that process needs to be resolved to pick the equilibrium that is reached.

An obvious fast time scale process is the adjustment of lake level. If $q > 0$, then lake level attains seal height at leading order (Appendix A), but more accurately, there is a small $O(\nu )$ difference between the two: the bigger the flux $q$, the more lake level exceeds the seal height by, even if that amount remains small compared with the $O(1)$ height of the seal itself. More specifically, solving a boundary layer problem around the seal (Appendix B.2 or, in more detail, in §§ 1.5–1.6 of the supplementary material) allows us to compute (although not in closed form) flux $q$ in terms of the difference between lake water level $h_0$ and seal height $b_{m}$, and in terms of the slopes $b_x^+$ and $b_x^-$ up- and downstream of the seal. That relationship can be expressed in the form

(3.16)\begin{equation} q = \mathcal{Q}_s( \nu^{ - 1}(h_0-b_{m}),b_{x}^-,b_{x}^+), \end{equation}

where $\mathcal {Q}_s$ is non-negative, vanishes when its first argument is negative or zero (meaning, water level is at or below seal height) and is otherwise $O(1)$ when its arguments are $O(1)$ (see Appendix B.2, figure 15b).

When applied to the mass balance of the lake, this prescription for flux leads to a regularized version of (2.17)

(3.17)\begin{equation} \gamma h_{0,t} = Q(t) - q, \quad q = \mathcal{Q}_s( \nu^{ - 1}(h_0-b_{m}),b_{x}^-,b_{x}^+). \end{equation}

Note that this replaces the cruder but structurally similar ordinary differential equation models for lake surface lowering in Raymond & Nolan (2000), Kingslake et al. (2015) and Ancey et al. (2019).

Let $h_0 = b_{m} + \nu h_1$, so $h_1$ is the appropriately rescaled water level elevation above the seal. Then, using (3.10) to re-write $h_{0,t} = \dot {b}_{m} + \nu h_{1,t}$, we obtain from (3.17)

(3.18)\begin{equation} \nu \gamma h_{1,t}= Q - \gamma w - \frac{\gamma b_x^- M( - b_x^-,\mathcal{Q}_s(h_1,b_x^-,b_x^+))}{b_x^+ - b_x^-} - \mathcal{Q}_s(h_1,b_x^-,b_x^+). \end{equation}

Reassuringly, we recover (3.14) for $q = \mathcal {Q}_s(h_1,b_x^-,b_x^+))$ at leading order in $\nu$. The flux $q$ implicitly defines the correction $h_1$, but this still does not resolve the multi-valuedness of the solution. Omitting the time derivative $\nu h_{1,t}$ is, however, a singular perturbation that neglects transient behaviour on the faster $O(\nu )$ time scale: rescaling time as $T = \nu ^{-1} t$ in (3.17) yields an ordinary differential equation for $h_1$ with $b_x^-$ and $b_x^+$ constant on that fast time scale. $h_1$ will evolve to a stable steady state in $T$ that satisfies either case in (3.14).

(There is a slight inconsistency here: (3.17)$_2$ is the result of the steady state boundary layer problem in Appendix B.2. When rescaling to the fast time scale $T$, that boundary layer problem actually becomes non-steady, so that $q$ is determined by (B2) with an additional term $B_T$ added to the left-hand side of (B2b), meaning the boundary layer is no longer necessarily in steady state, in which case flux cannot necessarily be expressed as a function of $h_1$ and the slopes $b_x^-$ and $b_x^+$ only. We leave an analysis of that situation to future work.)

Assuming that flux increases with water level above the seal, we have $\partial \mathcal {Q}_s/\partial h_1 > 0$ as indicated by the numerical solutions in figure 15(b). The stability of steady-state solutions to (3.18) is then easy to determine: when there are three solutions to (3.14), only the largest solution (for which $q$ increases with base outflow rate $Q - \gamma w(x_{m})$) and the zero solution are stable as indicated by solid lines in figure 5(a).

In common with analogous problems such as glacier surges or other relaxation oscillators, the relevant, stable solution branch of the original leading-order model (2.17) is chosen by continuity of $q$ in the original, slow time variable $t$ whenever such continuity is possible. (This is true at least provided there are no significant variations in water supply $Q$ on the short $\sim O(\nu )$ time scale over which the lake level correction $h_1$ 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_1$ 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}) = U b_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 6a). 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.

Figure 6. Contour plots of $q$ as a function of $(b_x^+,b_x^-)$ for steady-state upstream conditions $w(x_{m}) = b_x^-/U$. Logarithmically spaced contour intervals with five contours per decade, contour levels as indicated by the colour bars. The dashed contour in each case corresponds to $q = Q$, at which the shock is stationary, migrating forward for $q < Q$ and backward for $q > Q$. The solid black curve indicates the boundary of the region in which only the zero solution exists. Panel (a) shows $U = 1$, $\gamma = 1$, $Q = 1$ $\alpha = 0.5$, red dot-dashed line is the lower boundary of the region in which $q = 0$ is a solution. Panel (b) shows $U = 1$, $\gamma = 4$, $Q = 2$, $\alpha = 0$; the solid red curve is the boundary of the region in which no stable solution for $q$ exists.

3.6. Flux for the constant channel width case $\alpha = 0$

For $\alpha = 0$, the situation is qualitatively different: we have $M(-p,q) = -pq$ and hence (3.14) reads

(3.19)\begin{equation} q (b_x^+ - b_x^- \gamma b_x^- b_x^+ )/(b_x^+ - b_x^- ) = Q - \gamma w(x_{m}), \end{equation}

if $q > 0$, and $Q - \gamma w(x_m) \leq 0$ if $q = 0$. Stability again requires $q$ to increase with $Q-\gamma w(x_{m})$, so the coefficient of $q$ on the left of (3.19) must be positive when a stable non-zero solution exists. As a result, there is no multi-valuedness: $q$ vanishes if and only if the base outflow rate $Q-\gamma w(x_{m}) < 0$. Instead of multi-valuedness, there may, however, be no solution at all. This occurs when

(3.20a,b)\begin{equation} Q-\gamma w(x_{m}) > 0 \quad \text{and} \quad (b_x^+ - b_x^- \gamma b_x^- b_x^+ ) < 0. \end{equation}

We can again visualize the flux solution situation, plotting $q$ against $b_x^+$ and $b_x^-$ under the assumption that $w(x_{m}) = Ub_x^-$ (figure 6b). Unlike the case $\alpha > 0$, $q$ is indeed not multi-valued, but instead undefined in a ‘forbidden’ part of the $(b_x^+,b_x^-)$ plane if $Q > U$ as shown, reflecting the solvability condition (3.20a,b).

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 (3.19) vanishes, and $q$ undergoes runaway growth (as the red boundary of the forbidden region is approached from below in figure 6b). We will show in § 4.2 below that this runaway growth corresponds to abrupt lake drainage, with a short-lived but finite jump in water height across the seal that cannot be captured fully by our reduced model.

4. Results

We solve (2.13)–(2.17) numerically using the method of characteristics with a backward Euler step as described in Appendix E. We use the regularized flux prescription (3.17) 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 (3.17), we treat $\mathcal {Q}_s$ simply as a regularization rather than trying to emulate the function shown in figure 15(b). Consequently we drop the slopes $b_x^-$ and $b_x^+$ as arguments from $\mathcal {Q}_s$. In practice, we use $\mathcal {Q}_s(h_1) = [\max (h_1,0)]^2$, and put $\nu = 10^{-3}$.

Figures 7–10 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

(4.1)\begin{equation} w(x) = U\left\{ \overline{b_x} -2 b_1 \lambda (x-x_0) \exp\left[-\lambda\left(x-x_0\right)^2\right]\right\}, \end{equation}

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 7a i–a ii). 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(a ii), or as one of the black curves in figure 7(a i,a ii).

Figure 7. Solutions for $\alpha = 0.5$: $\gamma = 1$, $Q = 0.1962$ (a i–c i) and $\gamma = 4$, $Q = 1.570$ (a ii–c ii). Panels (b i,b ii) show contour plots of $b(x,t)$, with $t$ on the horizontal and $x$ on the vertical axis, contour intervals of 0.1 and levels given by the colour bar. Blue lines show water level in the lake and boundaries of ponded sections, black lines show the smooth lake seal, magenta the closest shock to the seal, excluding the seals of ponded sections downstream. Inset in (b ii) shows detail of shock migration: note that the shock first forms downstream of seal, and then migrates upstream to incise the seal as predicted for $\alpha > 0$ in § 3.3. Panels (a i,a ii) show the profiles indicated by black dashed lines in (b i,b ii), respectively, with blue indicating water surface in the lake or a ponded section. Panels (c i,c ii) show time series of $x_{m}(t)$ (blue) and $q(t)$ (black), using the same $t$-axis as (b ii,b ii).

Figure 8. Solutions for $\alpha = 0.5$: $\gamma = 4$, $Q = 0.7850$ (a i–e i) and $\gamma = 2$, $Q = 0,7850$ (a ii–e ii). Same plotting scheme as figure 7, panels (a–c) now show profiles at equal time steps during the intervals between the vertical dashed lines in panels (d i,d ii), those intervals being marked with the appropriate panel label (a i)–(c i) and (a ii)–(c ii). The black dashed curve in panel (a ii) is the unincised ice surface $s(x)$.

We use two different choices of shape exponent, $\alpha = 1/2$ (the semicircles and triangles of (2.1c) and (2.1d), with wetted perimeter increasing monotonically with channel cross-section), and $\alpha = 0$ (the fixed wetted-perimeter slot of (2.1e)). 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 Clarke 1982, for a discussion of lake hypsometries). Both of these assumptions are simplistic, but help understand the dynamics of surface lakes more clearly.

4.1. Lake drainage modes for steady water supply: $\alpha > 0$

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-inflow 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 (3.11). 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-inflow counterexample to the steady state of column 1, with $Q = 1.570$ and $\gamma = 4$. Here, a shock forms quickly: the inset in panel (b ii) shows that the shock (magenta) forms downstream of the smooth seal (black) as predicted at the end of § 3.3, and subsequently migrates upstream to breach the lake. This causes flux $q$ to increase as stored lake water is released. The downstream 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 a ii), 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 - c M_{-p}$ by (3.3). 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 § 6 below.

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 a i). 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 b i). Panel (d i) 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 c i,e i).

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 § 6 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.

A more systematic exploration of the parameter dependence of lake drainage styles is shown in figure 9, with each column corresponding to a fixed value of $Q$ and each row to a fixed value of $\gamma$. The figure indicates that inflow rate $Q$ alone determines whether the seal is breached: a critical value of $Q$ appears to separate 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$.

Figure 9. Solutions for $\alpha = 0.5$. Time series of $q(t)$ (black) and $x_{m}(t)$ (blue) (as shown in e.g. panel (c) of figure 7) for different combinations of $\gamma$ and $Q$. $\gamma = 0.5$ (a i–a v), $1$ (b i–b v), $2$ (c i–c v) $4$ (d i–d v) and $Q = 0.1962$ (a i,b i,c i,d i) $0.3525$ (a ii,b ii,c ii,d ii) $0.4371$ (a iii,b iii,c iii,d iii) $0.7850$ (a iv,b iv,c iv,d iv) and $1.570$ (a v,b v,c v,d v). The solutions in figures 7(a i–c i) and 7(a ii–c ii) and 8(a i–e i) and 8(a ii–e ii), are shown in panels (b i), (d v), (d iv) and (c iv), respectively. Note that the critical water input for seal breaching predicted by (5.7) is $Q_c = 0.3917$, between columns (a ii–d ii) and (a iii–d iii) here. This also marks the transition from steady outflow $q$ to outflow $q$ increasing after a seal breach in this figure.

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.

4.2. Lake drainage modes for steady water supply: $\alpha = 0$

The range of drainage styles observed for $\alpha = 0$ is more limited. At low water input $Q < U = 1$, the channel develops into a steady state in much the same way as shown in figure 7, column 1. At larger $Q > U$, a shock once more forms, although this time at the seal in agreement with §§ 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 10.

Figure 10. Solutions for $\alpha = 0$: $\gamma = 2$, $Q = 1.1$ using (3.17) with $\nu = 5\times 10^{-3}$ (a i–c i) and $\gamma = 2$, $Q = 2$ (a ii–c ii). Same plotting scheme a figure 7. In panel (c i), we show two solutions, one using (3.17) as in panels (a i)–(b i) ($q$ in black, $x_{m}$ in blue), and the other without using that regularization, ($q$ in red, $x_{m}$ in magenta). The latter solution fails to converge numerically after $t = 4.764$.

Column 1 shows a case with more moderate water input $Q = 1.1$ and $\gamma = 2$. Panel (c i) shows results for discharge $q(t)$ and seal position $x_{m}(t)$ from two computations: one without the regularization advocated in (3.17) (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 § 3.4 (see in particular figure 6b): 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 the 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 10, 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 any instances of oscillatory lake drainage for $\alpha = 0$, as detailed in the more systematic exploration of the effect of changing storage capacity $\gamma$ and water input $Q$ in § 5.1 of the supplementary material.

5. Criteria for lake drainage

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$, see figure 9). The shock migrates backwards, leading to at least partial lake drainage.

Below, we identify situations in which such shocks must form with parameter combinations for which there is no steady-state solution to (2.13)–(2.17) (see also § 4 of the supplementary material). Steady states are the natural consequence of the seal remaining intact: as discussed after (3.7), the constant upstream boundary condition $b(0,t) = b_{in}(0)$ ensures that

(5.1)\begin{equation} b_t = - \mathcal{H} = 0, \end{equation}

on any characteristic originating at the upstream end of the boundary. This remains true if that characteristic passes a steady-state smooth lake seal at $\bar {x}_m$, across which $\mathcal {H}$ is continuous (§ 3.3), and below which $\mathcal {H}_\tau = \mathcal {H}_q q_\tau = 0$ if the seal is steady and hence $q_\tau = 0$. The formation of steady states by characteristics entering the domain from above 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. In other words, a global steady-state results if the characteristics crossing a steady seal can fill the entire domain, while non-existence of steady states implies that such characteristics cannot propagate through the entire domain.

A critical flux $Q$ beyond which steady states fail to exist is easy to identify. To simplify matters, recall that we assume the uplift function to have a single root $w(\bar {x}_{m}) = 0$ at the steady seal location, with $w < 0$ for $x > \bar {x}_{m}$ (as is indeed the case for the uplift function in (4.1)). In steady state, there is then a single ponded section with $b_x > 0$ upstream of $\bar {x}_{m}$, so we can omit the ponding function $c$ from the definition of the Hamiltonian in steady state.

The fixed-width channel case $\alpha = 0$ of figure 10 differs qualitatively in terms of shock formation from the variable-width case $\alpha > 0$ of figures 7–9, and we have to treat the two separately. First, consider $\alpha = 0$. Then, $M(-b_x,q) = -H(-b_x)b_xq$ where $H$ is the usual Heaviside function. In steady state, (2.17) demands that $q = Q$, while the steady state (5.1) becomes (figure 11b)

(5.2)\begin{equation} \mathcal{H} = (U - Q H( - b_x) ) b_x - w = 0. \end{equation}

We can solve for $b_x$ everywhere when $Q < U$. By contrast, $\mathcal {H}$ cannot be zero in regions where $w < 0$ (that is, downstream of the seal) if

(5.3)\begin{equation} Q \geq U. \end{equation}

This is the criterion for lake drainage when $\alpha = 0$, and is consistent with the observation in § 3.3 that a smooth lake seal cannot persist if (5.3) holds.

Figure 11. The reduced Hamiltonian $\mathcal {H}_r = \mathcal {H}+w$ for (a) $\alpha = 0.5$ and (b) $\alpha = 0$, with $U = 1$, shown for $q = 0.5$, $0.75$, $1$, $2$. The red dots in panel (a) correspond to $p = p_c$, $\mathcal {H}_r = \mathcal {H}_c$. Steady states satisfy $\mathcal {H}_r = w$, with negative $w$ found downstream of the seal. Combinations of $p$ and $\mathcal {H}_r$ shown as dotted curves correspond to backward-propagating characteristics. We expect steady states to remain purely on the dashed/solid branches of the curves: backward-propagating characteristics in steady state require steady-state boundary conditions at the downstream end of the domain, and must meet forward-propagating characteristics at a stationary shock, an unlikely scenario.

Second, consider the variable-width channel case with $0 < \alpha < 1$. In that case, we can define a reduced Hamiltonian $\mathcal {H}_r$ through

(5.4)\begin{equation} \mathcal{H}_r(x,t,p,q) = Up + M( - p,q), \end{equation}

so $\mathcal {H} = \mathcal {H}_r - w$. The reduced Hamiltonian (see figure 11a) has a global minimum $\mathcal {H}_c$ with respect to $p$ at a value $p = p_c$ given by

(5.5a,b)\begin{equation} p_c(q) = - \frac{[(3-\alpha)U/3]^{(3-\alpha)/\alpha}}{q^{3(1-\alpha)/\alpha}}, \quad \mathcal{H}_c(q) = - \frac{\alpha(3-\alpha)^{(3-\alpha)/\alpha}U^{3/\alpha}}{3^{3/\alpha} q^{3(1-\alpha)/\alpha}}. \end{equation}

(More generally, such a minimum $\mathcal {H}_c$ can always be identified for generic melt rate functions $M(-p,q)$ that are strictly convex functions of slope $-p$ for downward slopes $p < 0$, with $M(0,q) = M_{-p}(0,q) = 0$.)

A steady state with $\mathcal {H} = \mathcal {H}_r - w \equiv 0$ exists if and only if $\inf (w) \geq \mathcal {H}_c$, or equally, we infer that lake drainage occurs if

(5.6)\begin{equation} \inf(w) < \mathcal{H}_c(Q). \end{equation}

If (5.6) is satisfied, the combined effect of downward motion $w$ of the ice and channel incision $M(-b_x,q)$ must overwhelm downstream advection $Ub_x$, no matter the channel slope $b_x$, and a steady state cannot be established. For the melt rate function given by (2.12), the criterion (5.6) can be re-written in the form

(5.7)\begin{equation} Q > \frac{\alpha^{\alpha/[3(1-\alpha)]}(3-\alpha)^{(3-\alpha)/[3(1-\alpha)]}}{3^{1/(1-\alpha)}}[-\inf(w)]^{-\alpha/[3(1-\alpha)]} U^{1/(1-\alpha)}, \end{equation}

which gives the desired critical flux for breaching the seal. While this differs from the criterion (5.3) for $\alpha = 0$, note that (5.7) reassuringly does reduce to $Q > U$ in the limit $\alpha \rightarrow 0$. Note also that (5.7) is consistent with our numerical results: the critical flux is $Q= 0.3917$ for the calculations in figure 9 and $Q = 1$ for the results reported in § 4.2 as well as in § 5.1 of the supplementary material.

6. Oscillatory lake drainage

Breaching of the seal need not lead to complete emptying of the lake: the lake can re-seal and re-fill temporarily instead (figure 8). Re-sealing results from changes in upstream slopes $b_x^-$ and $b_x^+$ at the seal during lake drainage, whose effect on $q$ is shown generically in figure 5. We have observed partial lake drainage only for $\alpha > 0$, as shown in figure 5(a). We superimpose ‘orbits’ of $(b_x^+,b_x^-)$ during different lake drainage events in figure 12 to track the effect of slopes and their role in re-sealing the lake.

Figure 12. ‘Orbits’ of $(b_x^+,b_x^-)$ at the shock that breaches the seal, superimposed on corresponding versions of the flux contour plot in figure 5(a), with contour lines for flux rendered in grey, These orbits are shown for the solutions in (a) figure 8(a i–e i), showing one of the periodic drainage cycles (b) figure 8(a ii–e ii) and (c) figure 7(a ii–c ii), showing the full solution for the latter two. The curves are colour coded by time as shown in each colour bar. The ‘orbit’ penetrates perceptibly into the zero flux ($q=0$) region at the top right of panel (a) because of the regularization (3.17) used in the computation of the time-dependent solution. The orbits terminating at $b_x^+ < 0$, $b_x^- = 0$ in panels (b,c) correspond to the shock reaching the bottom of the lake at the upstream end of the domain.

A steeper downstream slope $b_x^+ < 0$ leads to faster incision into the seal, and therefore to a greater rate of backward migration of the seal and hence of lake drainage at fixed upstream slope, so $q$ increases with decreasing $b_x^+$. The upstream slope $b_x^-$ has two conflicting effects: larger $b_x^-$ on the one hand slows the backward migration of the seal (through the denominator on the left of (3.14)) and corresponds to a greater rate of uplift, trying to re-seal the lake. On the other, for a given backward migration rate of the seal, a steeper upstream slope also corresponds to a greater rate of surface lowering and therefore volume loss from the lake at a given rate of seal migration. The latter effect dominates for steeper (more negative) downstream slopes $b_x^+$, the former for shallower $b_x^+$.

In the solutions we have reported above, termination of flow before the lake is fully empty generally hinges on two effects. First, while $b_x^-$ initially steepens after incision into the smooth seal, the upstream slope eventually flattens again after an inflection point in the unincised surface $s(x)$ is passed, and approaches zero as the seal point $x_{m}(t)$ approaches the bottom of the lake at $x = 0$, causing $q$ to decrease again (figure 12a,b).

Second, following an initial decrease, the downstream slope $b_x^+ < 0$ eventually increases (becoming less negative) during lake drainage, as already mentioned in § 4. The mechanism involved is the following: as incision into the seal occurs, $q$ initially increases. The increase in flux causes characteristic velocities downstream of the seal to become more negative by (3.3), so characteristics propagate upstream faster. As described in § 4, faster propagation of characteristics can cause a reduction in slopes: slope evolves as $p_\tau = w_x$ along characteristics, and $w_x$ is typically most negative around the original smooth seal location $\bar {x}_m$. The faster that characteristics move through this region of steepening because flux $q$ has increased, the less $p = b_x$ will steepen. As a result, characteristics that reach the shock at the seal $x_{m}(t)$ later during lake drainage do so with a less steep (that is, less negative) slope $b_x^+$.

Figure 12 shows that the increase in downstream slope $b_x^+$ (that is, reduction in magnitude $|b_x^+|$) is key in ensuring that flux is not only reduced in the later stages of lake drainage (as a flattening of $b_x^-$ already ensures), but actually vanishes entirely on reaching the boundary of the blank region of zero flux (marked $q = 0$ in the equivalent figure 6): compare panels (a,b) of figure 12 with panel (c), which shows the equivalent orbit for a lake that drains completely after the initial seal incision. Reaching that boundary in figure 12 implies that flow ceases abruptly, and the lake re-seals.

The case shown in figure 12(a) is additionally visualized in figure 13, where we show characteristics that reach the shock from downstream as multicoloured curves, the colouring indicating time (panel a) or slope $p = b_x$ along the characteristic (panel b). The increasingly rapid transit of characteristics past the point $x = \bar {x}_m$ (vertical dotted line marked $S_1$ in panel a) and the reduced steepening at later times during lake drainage is evident in panels (a) (later characteristics do not dip to larger negative values of $b_x$, and the colouring indicates only a short amount of time spent near $\bar {x}_m$) and (b) (later characteristics are steeper near $\bar {x}_m = 1.468$, indicating a faster passage, and retain a lighter blue colour indicating less steep slopes).

Figure 13. The solution in figure 8(a –e i) and figure 12(a). Panel (a) shows slope against position. The solid black curve is slope $b_x(x,t)$ against $x$ at $t = 41.1$, when lake level reaches the height of the smooth seal and lake discharge recommences. The dashed black line, partially obscured by the solid curve, is the slope $s_x$ of the unincised ice surface. Purple and red curves (solid and dashed) show the trajectory taken by downstream slope $b_x^+$ on the new shock that forms after flow commences (the upstream slope is $b_x^- = s_x$), the red arrow indicating the direction in which the shock traverses the curve as time $t$ increases. Purple indicates that the shock is downstream of the smooth seal at $\bar {x}_m$ (dotted vertical line), red indicates the shock has incised into the seal. Dashes (between points $C$ and $S_2$) indicate that there is zero flux $q$, while the solid portion of the curve between $A$ and $C$ corresponds to positive flux. The multi-coloured curves are characteristics that arrive at the shock at intervals of $\delta t = 0.1125$ while the seal is breached and water is flowing, coloured shading indicates time. Panel (b) shows same information but plotted as position against time, with coloured shading indicating the slope $b_x$ on the characteristics. The blue curve shows flux $q$ against time, plotted using the right-hand vertical axis tick marks. Here, $S_1$ marks the smooth seal where $w(\bar {x}_m) = Us_x(\bar {x}_m) = 0$ as indicated by the horizontal dotted line; $S_2$ marks the shock left by the previous drainage cycle. The point labels $A$–$D$ mark changes in the shock, from formation at $A$ to breaching the smooth seal at $B$, flow ceasing at $C$ to a new smooth seal forming as the shock passes the smooth seal location $x = \bar {x}_m$ at point $D$. Note that the dashed portion of the curve from $C$ to $S_2$ is a translated version of part of the black initial profile curve on which points $S_1$ and $S_2$ lie; this is no accident since both are characteristics with the same characteristic velocity $x_\tau = U$ and evolution equation $p_\tau = w_x$.

The effect of downstream flattening during seal incision becomes stronger if storage volume $\gamma$ is large or the inflow $Q$ is smaller (but still above the critical value for the initiation of drainage as discussed in § 5). Both larger $\gamma$ or smaller $Q$ lead to a bigger relative increase in flux $q$ during lake drainage, and hence to a stronger relative flattening of the downstream slope. This accounts for oscillatory drainage occurring at such parameter combinations in figure 9.

Spontaneous termination of lake drainage, however, need not lead to periodically recurring lake drainage, see for instance figures 8(d ii) and 12(b): consecutive filling and drainage cycles may have an increasing amplitude, leading to complete lake emptying eventually. This appears to be linked to rapid re-filling of the lake and re-activation of the same shock that caused the initial incision. Once the reactivated shock incises the smooth seal again, it may do so with a steeper downstream slope and incise further upstream (figure 12b).

Reactivation of the same shock is favoured by small lake storage capacity and larger fluxes $Q$, which allow the lake to refill rapidly. As a result, the shock that originally incised the smooth seal is not advected far enough downstream, and on reactivation reaches the re-formed smooth seal again. Periodic lake drainage by contrast results most easily if $\gamma$ is larger and inflow rates $Q$ are small but above the critical value for drainage.

In that case, lake re-filling takes longer and the shock that incised the seal on the previous drainage cycle is advected far enough downstream between cycles for it not to return to incise the seal again. This is illustrated in figure 13. The ice surface rebuilds to a local steady-state solution $Ub_x = w$ everywhere upstream of the advected shock by the time the lake refills and outflow of water recommences (the black curve in figure 13(a), with the advected shock being marked by $S_2$; the dashed black curve continues the steady state solution $Ub_x = w$ past the advected shock, where it now represents the unincised ice surface $s(x)$). When flow of water recommences, a new channel is incised and a new shock is formed in this previously steady part of the ice surface (the pink–red line originating at point $A$ in figure 13(a), see also the pink line in figure 8d i). This new shock migrates upstream, intersecting the rebuilt smooth seal $S_1$ at point $B$ in figure 13(a). Crucially, all characteristics that reach the new shock from downstream during the drainage cycle also originate upstream of the old shock (the coloured lines in figure 13(a), all of which start upstream of $S_2$). Once flow terminates again (point $C$), the new shock is in turn advected downstream, with a new smooth seal forming at point $D$. The new shock reaches the position $S_2$ of the previous shock at the start of the next cycle, which repeats the previous one exactly.

It is worth emphasizing that the flood termination mechanism described above involves a surface shape $s(x)$ that flattens upstream of the seal, as is likely to be generically the case for surface lakes on ice sheets that form due to a smooth local uplift anomaly. Lakes whose bottom does not flatten out do occur on mountain glaciers, for instance at glacier confluences (Werder, Loye & Funk 2009). We investigate this situation in § 5.3 of the supplementary material, prescribing an uplift velocity $w(x)$ such that $b_x$ tends to a positive constant far upstream of the smooth lake seal (that is, the lake surface does not flatten). Repeated floods with constant lake supply $Q$ are much less common, and follow a somewhat different mechanism: the same shock reactivates in each cycle but does not steepen from cycle to cycle in such a way as to cause complete lake drainage.

7. Discussion and conclusions

In this paper, we have derived and solved a reduced, ‘stream power’-type model (Whipple & Tucker 1999) for supraglacial stream incision (2.13), coupled to a model for lake drainage to determine the water flux $q$ (2.17), whose value depends on whether and how fast the stream is cutting into the lake seal. Note that, for completeness, the model is stated in dimensional form in § 1 of the supplementary material. At the most basic level, the model predicts that a lake drains if water input to the lake is sufficiently large to overcome the effect of forward advection of the channel by the flow of the ice: if the inflow criterion (5.7) is satisfied (again stated in dimensional form in § 1 of the supplementary material), then the incision of the outflow channel will cause the lake seal to be breached eventually by a backward-propagating shock. The criterion demonstrates that sufficiently large water supply, steep downward slopes on the far side of the seal (large $-\inf (w)$, where $w$ is the uplift velocity of the ice) and slow advection (small values of the horizontal velocity $U$) is key to lake drainage. In particular, forward advection of the channel is the critical difference between the supraglacial lake drainage case and other dam burst phenomena (e.g. Balmforth, von Hardenberg & Zammett 2009). Qualitatively, our lake drainage criterion (5.7) is at least consistent with the observation (Poinar & Andrews 2021) that non-draining lakes in Greenland are located at higher elevations (where water supply rates will be smaller, as are vertical velocities $w$) compared with ‘slowly draining lakes’, which may conceivably drain through surface channels rather than hydrofracture. Our model also suggests that, as surface melt rates and therefore rates of water supply $Q$ continue to increase, more lakes should eventually drain.

The model also predicts that initial incision into the seal need not lead to complete lake drainage. Instead, a flattening of both upstream and downstream slopes at the shock at the downstream end of the lake can lead to the lake re-sealing, with forward advection of the shock subsequently causing the original lake basin to re-form. The flattening of the downstream slope is facilitated not only by relatively slow water inflow rates to the lake but also, and perhaps counterintuitively, by a large lake storage capacity, with both facilitating a large relative increase in discharge during lake drainage and rapid retreat of the lake-terminating shock that ultimately causes the slopes downstream of the shock to flatten again (§ 6).

The dynamics of supraglacial lakes in our model ultimately permits four different outcomes: no incision of the seal (at inflow rates below the critical value given by condition (5.7)), a periodic cycle of the lake being breached and draining, followed by refilling (at large storage capacity and small above-critical water inflow), a sequence of lake drainage episodes of growing amplitude that progress until the lake fully empties (at intermediate storage volumes and water supply rates) and complete lake drainage at small storage volumes and large water supply rates. The possibility of oscillatory lake behaviour by overland drainage in particular has implications for the interpretation of lake observations by remote sensing, where the drainage mechanism may not be immediately apparent: it permits lakes to drain ‘unexpectedly’, that is, provided that seal incision is already underway, they may drain outside the melt season (e.g. Schaap et al. 2020; Benedek & Willis 2021), and for lakes to remain filled for multiple melt seasons until the seal is breached again (see also § 5.2 of the supplementary material). However, unlike the condition (5.7) for seal breaching, we are unable to give simple criteria for complete vs partial, oscillatory lake drainage; presumably, the boundaries between these phenomena in parameter space depends on the specifics of the uplift velocity $w(x)$.

Importantly, the results we have presented assume steady water supply, while 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.2 of the supplementary material, where we find that it leaves the qualitative behaviour of the system largely unchanged.

One shortcoming of our model relevant to the different drainage styles above is its one-dimensional nature. Implicit here is that, even if drainage terminates and the lake re-seals, subsequent overflowing of the reconstituted lake will re-activate the same channel as before. This is key to the drainage cycles with growing amplitude, leading to the lake emptying fully (column 2 of figure 8): the re-activation of the previous channel leads to subsequent drainage of the lake progressing further. If instead the previous channel is advected laterally as well as down-slope (Darnell et al. 2013), then a new channel may be formed each time and periodic drainage cycles may in fact be more common than our results indicate.

More broadly, it is worth revisiting the construction of our model. The glacial case is perhaps the simplest in which a ‘stream power’ model for channel incision can be justified from first principles: the product of ‘erosion’ by flowing water is simply more water, rather than sediment whose transport must then be accounted for (Fowler et al. 2007). Our results, however, do indicate that the model as stated is incomplete: the predictions of the model depend strongly on the choice of the exponent $\alpha$ that parameterizes the cross-sectional shape of the channel in our model (§ 2.1). For instance, for channels of fixed width (independently of their cross-sectional area) we have $\alpha = 0$. Unlike the case of channels with variable width ($\alpha > 0$), we have found no oscillatory lake drainage (§§ 4 and 6) when $\alpha = 0$. Instead, the lake can drain ‘abruptly’ in the sense that flux becomes large, incision becomes rapid and water level in the lake does not remain close to the height of the seal as stipulated by (2.17); in the model consisting of (2.13)–(2.17), this phenomenon manifests itself as flux becoming singular unless (2.17) is regularized (§§ 3.4–4).

To determine even the qualitative behaviour of lake drainage unambiguously, a more sophisticated model for channel evolution is therefore necessary, capable of predicting the shape of cross-sections self-consistently instead of imposing it as a constitutive relation. There is currently no particularly good template, although the work in Dallaston & Hewitt (2014) may be a good starting point. Closely linked to cross-sectional shape evolution is the need to be able to predict meandering (Karlstrom et al. 2013; Fernández & Parker 2021), which ultimately should modify our large scale model (2.13) through the introduction of an evolving tortuosity. Not only is a model for cross-sectional shape now required, but the secondary flows involved in meandering instabilities also need to be accounted for, which also occur at wavelengths comparable to channel width (Karlstrom et al. 2013). (That being said, it is worth remembering that even the more sophisticated subglacial drainage models in existence (e.g. Werder et al. 2013) do not attempt to account for evolving tortuosity.) Lastly, the ability to account not only for lateral instabilities driving meandering, but also for bedform formation and roll waves at supercritical Froude numbers (Fowler 2011, see also § 3.2 of the supplementary material) is also desirable, in order to be able to apply the model on steeper slopes or at large discharge rates.

There are numerous other shortcomings to our model as described in § 2.1, such as the neglect of melting due to heat exchange with the atmosphere and solar radiation, accumulation of snow in the channel and exchange of water with an underground firn aquifer. We conclude by pointing out that, these issues notwithstanding, our model provides a template for improving previous surface drainage models due to Raymond & Nolan (2000) and Kingslake et al. (2015). As with the prior, although slightly different work in Walder & Costa (1996) (which considers the widening rather than deepening of a pre-existing breach through the full thickness of an ice dam), the models for downward incision of a channel in Raymond & Nolan (2000) and Kingslake et al. (2015) are heavily parameterized and do not resolve position along the channel. In effect, they are ad hoc versions of the boundary layer problem in our Appendix B.2, aiming to compute the function $\mathcal {Q}_s$ of § 3.4 here: Raymond & Nolan (2000) equate the difference between lake level and seal height ($H_w(-\infty,t)$ in Appendix B.2) with the far-field water depth in the same boundary layer (our $\varSigma (\infty,t)^\beta$ in Appendix B.2), while Kingslake et al. (2015) questionably impose Bernoulli's law (valid in the inertia-dominated upstream far field of the boundary layer) at the same time as a balance between the downslope force of gravity and friction at the channel wall (valid in the friction-dominated downstream far field). The details of those calculations aside, it is unclear whether the precise form of the relationship between flux and water height above the lake seal are key to modelling a supraglacial outburst flood: our work suggests that it may often (except in the flux singularity case for fixed-width channels illustrated in figure 10a i–c i) suffice to require that lake level remains approximately at the seal, and to focus instead on the incision of the channel over longer length scales, which allows the channel slope at the shock-like lake seal to change as the outburst flood progresses, changing the rate of backward migration of the seal and hence the rate of lake drainage.