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 Reference Fowler2011, 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 Reference Fowler2011)

(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. Reference Evatt, Fowler, Clark and Hulton2006; Ogier et al. Reference Ogier, Werder, Huss, Kull, Hodel and Farinotti2021); 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 Reference Dallaston and Hewitt2014) 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 Reference Schoof2002), 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 Reference Jarosch and Gudmundsson2012) as well as snow accumulation at the base of the channel during winter, but also that $s$ is in steady state (Schoof Reference Schoof2002), 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 Reference Karlstrom, Gajjar and Manga2013; Fernández & Parker Reference Fernández and Parker2021). 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 Reference Buzzard, Feltham and Flocco2018). 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 Reference Raymond and Nolan2000). 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 Reference Fowler2011). 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 Reference Holmes1995). 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 (Reference Luke1972), Fowler, Kopteva & Oakley (Reference Fowler, Kopteva and Oakley2007) and Kwang & Parker (Reference Kwang and Parker2017).

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.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 Reference Luke1972),

(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 Reference Courant and Hilbert1989, § 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 Reference Luke1972). 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 Reference Luke1972; Royden & Perron Reference Royden and Perron2007).

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 Reference Royden and Perron2007)

(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 Reference Fowler1987,  Reference Fowler1989), or the behaviour of the van der Pol oscillator in the relaxation oscillation limit (Holmes Reference Holmes1995), 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 (Reference Raymond and Nolan2000), Kingslake et al. (Reference Kingslake, Ng and Sole2015) and Ancey et al. (Reference Ancey, Bardou, Funk, Huss, Werder and Trewhela2019).

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.