An analysis of instabilities and limit cycles in glacier-dammed reservoirs

Glacier lake outburst floods are common glacial hazards around the world. How big such floods can become (either in terms of peak discharge or in terms of total volume released) depends on how they are initiated: what causes the runaway enlargement of a subglacial or other conduit to start, and how big can the lake get before that point is reached? Here we investigate how the spontaneous channelization of a linked-cavity drainage system controls the onset of floods. In agreement with 5 previous work, we show that floods only occur in a band of water throughput rates, and identify stabilizing mechanisms that allow steady drainage of an ice-dammed reservoir. We also show how stable limit cycle solutions emerge from the instability, a show how and why the stability properties of a drainage system with spatially spread-out water storage differ from those where storage is localized in a single reservoir or ‘lake’.

widely adopted in the study of subglacial hydrology (Werder et al., 2013). We identifying instabilities that spontaneously lead to self-sustained oscillations, describing the mechanisms behind them and delineating the regions in parameter space where they occur. This is likely to be useful in diagnosing the behaviour of such models, regardless of their specific application to large outburst floods emanating from easily identifiable glacier 'lakes'.
To do so, we use a hierarchy of models, employing both, a spatially-extended, one dimensional drainage system model and 5 a lumped, box-type model that provides additional physical insight, as well as being the appropriate limit of the spatiallyextended model in the case of a long flow path. The paper is laid out as follows: in section 2, we develop the spatially extended and lumped models. In sections 3.1 and 3.2, we identify instabilities in the lumped model, and where they are suppressed.
'Instability' here refers to the system evolving away from a steady state if the latter is slightly perturbed. What that instability evolves into is investigated in detail in sections 3.3 and 3.4, where we show that stable limit cycles emerge in the lumped 10 model. In section 4, we show that the lumped model replicates the behaviour of the spatially extended behaviour well in the classical parameter regime that corresponds to the drainage of large glacier-dammed lakes. We also show that the extended model is unstable in more exotic parameter regimes, where the lumped model is not, and investigate the resulting dynamics.
These results are summarized in section 5.1, and we take the results for localized reservoirs as motivation to study stability of drainage systems with distributed storage in 5.2. Our results are summarized in section 6, and extensive additional information 15 are provided in the supplementary material, including the code used in computations.

A continuum model for outburst floods
We use the one-dimensional continuum model for drainage through a single conduit in Schoof et al. (2014), building on similar models used elsewhere (e.g. Ng, 1998;Hewitt and Fowler, 2008;Schuler and Fischer, 2009;Schoof, 2010;Hewitt 20 et al., 2012). Denoting conduit cross-section by S(x, t), effective pressure by N (x, t) and discharge by q(x, t), where x is downstream distance and t is time, we put on 0 < x < L, subject to boundary conditions where L is the length of the flow path, and we assume the closures Here c 1 , c 2 , c 3 , n, α, u b , h r and S 0 are positive constants with α > 1, while q in is also prescribed as a function of time and V p > 0 is a function of N . Ψ 0 is a geometrically determined background gradient, given in terms of ice surface elevation s(x) and bed elevation b(x) through Here ρ w and ρ i are the densities of ice and water, respectively, and g is acceleration due to gravity. Note that effective pressure N is linked to water pressures p w , which is the primary observable in the field, through where p i is overburden (or more precisely, normal stress in the ice at the bed, averaged over a length scale much larger than 10 that occupied by the channel). Typically (including in (1h)), p i is assumed to be cryostatic, equal to ρ i g(s − b).
Physically, the model represents a conduit whose size evolves due to a combination of dissipation-driven wall melting at a rate c 1 qΨ, opening due to ice sliding over bed roughness at rate u b h r (1 − S/S 0 ), and creep closure at rate c 2 SN n , with discharge in the conduit given by a Darcy-Weisbach or Manning friction law as c 3 S α |Ψ| −1/2 Ψ through (1g). u b is sliding velocity, h r is bed roughness and S 0 is a cut-off cavity size at which bed roughness is drowned out . S 0 is 15 typically a regularizing parameter that prevents conduits from becoming excessively large in regions where effective pressures are low, such as near the glacier terminus, but has little effect on conduit sizes along most of the flow path; this corresponds to the limit of large S 0 .
Vanishing effective pressure at the end of the flow path x = L simply reflects the assumption that overburden and water pressure vanish simultaneously (where we have neglected atmospheric pressure as a gauge throughout). The upstream boundary 20 condition at x = 0 instead reflects water storage in a reservoir at the head of the conduit: (1e) represents conservation of mass in the reservoir, with inflow rate q in , outflow at rate q(0, t) through the conduit modelled by (1a)-(1d) and with water storage V in the reservoir a function of effective pressure at x = 0 such that dV / dN = −V p . Note that we deliberately use 'reservoir' rather than 'lake' here since we will be interested not only in large lake, but also in more modest-sized reservoirs as in Schoof et al. (2014). We will treat V p and q in as constants in most of what follows. Conservation of mass along the conduit (1b) by 25 contrast assumes that water storage is negligible along the flow path. We show that the contribution of wall melting to mass flux is in general negligible in terrestrial drainage systems in section 2.2 of the supplementary material.
More generally, (1e) is justified as follows (Clarke, 2003). Water volume in the reservoir can be related to water depth h w and hence to water pressure p w0 at the bottom of the reservoir as dV / dh w = A, where A is the surface area that corresponds to the filling level h w . With p w0 = ρ w gh w , dV / dp w0 = A/(ρ w g), In general, A(h) is given by the shape of the reservoir. If we put N = p i − p w0 at the head of the conduit, then h w and therefore A are functions of N at x = 0. Equation (1e) follows with V p = A(h w )/(ρ w g); constant V p therefore corresponds to a lake with vertical sides, whose surface area is independent of filling level.
There is an important simplification we have made, irrespective of the particular choice of V p : Our model effectively imposes 5 an ice cliff at the conduit inlet x = 0, with the cliff typically above the flotation thickness. This ensures that N (0, t) can change over time without the upstream end of the conduit migrating. This contrasts with a glacier that partially floats on the lake, in which case the upstream end of the conduit is always at N = 0, but migrates as the lake fills or empties, and the location of that upstream end rather than effective pressure is related to lake volume. The filling level in the lake no longer dictates the effective pressure at the conduit inlet, but how much of the glacier has floated and therefore where that inlet has migrated to 10 (see section 2.4 of the supplementary material).
The assumption of an ice cliff, in addition, allows us to assume that Ψ 0 > 0 along the entire flow path, without a 'seal' at a finite distance from the reservoir at which Ψ 0 changes sign (Fowler, 1999). With a seal, there is a finite region of negative Ψ 0 between reservoir and seal location, with water flow out of the reservoir possible only when local effective pressure gradients ∂N/∂x are large enough to overcome the effects of the seal to make the total hydraulic gradient Ψ positive in (1d). In the 15 absence of a seal and without an ice cliff (so N = 0 at the edge of the reservoir), normal stress at the glacier bed just outside of the reservoir would then be below the flotation pressure, and there would be nothing to dam the reservoir (see section 2.4 of the supplementary material).
The two situations (an ice cliff and a seal) can be reconciled in the sense that a seal very close to the edge of the reservoir corresponds to a short, steep ice surface slope rather than an actual cliff. This results in a large, negative Ψ 0 between reservoir 20 and seal over a short distance, which can be balanced by an equally large ∂N/∂x over the same short distance, leading to a non-zero effective pressure at the seal itself, only a short distance from the lake.
Nevertheless, our simplifying assumption is relevant as the mechanism for initiating periodically recurring floods is fundamentally different from that in Fowler (1999): his model relies on a geometrical seal with negative Ψ 0 near the reservoir, changing sign at the seal location, and an englacial water supply to the conduit that dictates the gradient ∂N/∂x while the 25 lake is filling. As described in the introduction, the onset of the outburst flood corresponds to the instant at which Ψ at the the conduit inlet x = 0 changes from negative to positive, and allows water to flow out of the reservoir. This leads to runaway enlargement of the conduit through dissipation-driven melting.
By contrast, our model assumes that the reservoir always experiences some amount of leakage, and requires no englaciallysupplied conduit: the water supply to the drainage system in our model is purely through the inflow term q in into the reservoir.

30
Leakage out of the reservoir is facilitated by the linked-cavity behaviour of the conduit between floods, associated with the cavity opening term v o that keeps the conduit open even when flow rates are small. This is absent from the models in Fowler (1999) and Ng (1998). As the lake fills, these cavities enlarge because the effective pressure is reduced, and the conduit then grows through the same melt-driven channel enlargement as in Fowler (1999) and Ng (1998), but the initiation of that enlargement differs.

A lumped model
The model (1) can be simplified if we assume that Ψ 0 not only does not change sign along the flow path, but can be treated as a constant, and if we assume that the effective pressure gradient ∂N/∂x along the flow path can be treated as negligible (see also Fowler, 1999;Ng, 2000, and sections 2.2 and 2.3 of the supplementary material).
Under these assumptions, we can relate the flux q along the flow path (which is independent of position by (1b)) purely to 5 the conduit size S and hydraulic gradient Ψ 0 at the head of the channel, leading to a system of ordinary rather than partial differential equations: where a dot signifies an ordinary derivative with respect to time, and we continue to assume the closure relations (1g); by the argument above, we then strictly speaking have to put Ψ = Ψ 0 . We generalize this reduced model slightly to account qualitatively, if not quantitatively, for the effect of pressure gradients along the flow path by putting In what follows, we proceed first with an analysis of the simpler, 'lumped' model (2), to which we can bring to bear the theory of finite-dimensional dynamical systems (Wiggins, 2003), leading a number of semi-analytical results. Subsequently, we study the full, spatially-extended model (1), for which we are limited to comparatively expensive numerical methods that we can guide by our analysis of the simpler, lumped model. 3.1 Nye's jökulhlaup instability in the reduced model Nye's (1976) original theory of jökulhlaups centers on the idea that steady flow in a channel is unstable if there is a water reservoir that keeps water pressure approximately constant. This instability results in the runway growth of a channel, eventually allowing the reservoir to drain in an outburst flood. The assumption of a constant effective pressure is of course a simplification of the more complete model (2), based on the notion that the drainage of the lake is generally too slow to lead to changes in 25 effective pressure that could stabilize the channel against growth.

An analysis of the lumped model
The basis of this instability is relatively easy to capture mathematically. To set the scene, we present a simplified version first before tackling an analysis of the complete mdel (2). If we assume a classical 'channel' in the sense of Röthlisberger (1972) and Nye (1976) and therefore set the cavity opening term v o to zero in (2a), and use the remaining relations in (1g), the conduit with a steady state conduit sizeS given by c 2S |N | n−1N = c 1 c 3S α |Ψ| 3/2 . At fixed effective pressure N , and hence at fixed Ψ, this steady state is unstable: an increase S away from the steady state sizeS will lead to a further growth in he conduit, as where, with α > 1, the right-hand side is positive if S is, signifying that S will continue to grow in a positive feedback.
In more abstract terms, if Ψ and N are fixed, (2a) can be expected to lead to unstable growth of conduits away from an equilibrium stateS if (see also pages 5-6 of the supplementary material to Schoof (2010)) This is also precisely the condition that defines a conduit as being 'channel-like' in Schoof (2010), in the sense that two neighbouring conduits will compete for water with one eventually growing at the expense of the other if both are 'channellike'. A channel-like conduit is also distinguishable by the fact that the effective pressure required to balance opening processes (dissipation-driven wall melting and cavity opening due to flow over bed toughness) by creep closure increases with discharge 15 q in in the channel; the opposite is true for a 'cavity-like' conduit.
Key to the instability mechanism above was the notion that effective pressure N is kept constant as the conduit evolves, which is not actually the case: instead, the lake simply buffers changes in N from happening rapidly, since a change in N can only occur after an imbalance develops between inflow q in and outflow q. Next, we investigate in more detail how water storage and the finite size L of the system control whether Nye's instability does occur in the reduced model of section 2.2.

20
Our main result will be that this is not unconditionally the case, and that some reservoirs can drain steadily without outburst floods.
Note that v o and v c satisfy q also has the same sign as Ψ, and satisfies q(0, Ψ) = q(S, 0) = 0. These are the minimal assumptions we make on these functions, allowing us to generalize from the specific forms in (1g) in analyzing Nye's instability.
The primary dependent variables in the model (2) are S and N . For constant water input q in , the model admits a steady state (S,N ) given implicitly by If q in > 0, it follows thatΨ > 0, and from (8a) and the properties of v c , we also haveN > 0. For future convenience, we also write q(S.N ) =q. Note that, for the specific choices in (1g), a steady state exists for every positive q in . EliminatingS andΨ, we find a problem forN alone: The left-hand side is a monotonically decreasing function ofN for 0 ≤N < Ψ 0 L, positive atN = 0 and tending to −∞ as and puttingV p = V p (N ) gives the following leading-order form of (2): As before, we want to know whether S grows over time. Looking for solutions of the form S = S 0 exp(λt), N = N 0 exp(λt), Setting the determinant of the matrix on the left to zero leads to a polynomial for λ, where the coefficients take the form But, from our assumptions on the various functions involved, we see that a 2 > 0 (recall that v o,S ≤ 0 from (6)), while a 1 can be either sign.
The characteristic quadratic (10) has solutions Since a 2 > 0, we have a 2 1 − 4a 2 < a 2 1 and two possible types of solution: either a 2 1 − 4a 2 > 0 and we have two real roots, both 10 of which have the same sign as a 1 . Alternatively, we have a 2 1 − 4a 2 < 0 and a complex conjugate pair of roots, both of which have real part a 1 . In either case, we see that the system is linearly unstable if and only if a 1 > 0, or We have deliberately written the left-hand side of (13) as the difference of two terms, a potentially destabilizing term c 1 q SΨ + v o,S − v c,S and a stabilizing term −V −1 p q Ψ L −1 . The first term can be recognized as the growth rate sensitivity we 15 previously identified as being at the heart of Nye's instability in (5), as well as an indicator of whether the steady-state conduit is 'channel-like' in the terminology of Schoof (2010).
The second term in (13) is a stabilizing term that is inversely related to storage capacity. Its physical origin is the following: the sensitivity of channel growth to perturbations in conduit size may be positive, potentially leading to unstable conduit growth.
However, growth of the conduit also allows water to drain out of the system, which will increase the effective pressure N . As 20 the effective pressure is increased, the hydraulic gradient Ψ = Ψ 0 − N/L is reduced. This leads to less turbulent dissipation in the conduit as the conduit grows, and increased N further leads to faster creep closure of the channel. Both of these will suppress further growth of the conduit.
How strong this stabilizing effect is will depend on the storage capacity of the system. If the storage capacity is large, then the growth of the channel will have a minimal impact on effective pressure (essentially, the amount of water that drains out 25 due to the widened conduit is insufficient to affect water levels and therefore water pressure in the storage system sufficiently).
The stabilizing effect of the second term is therefore small when storage capacityV p is large. Similarly, if the system size L is large, then the role of effective pressure in controlling hydraulic gradients is small, and Ψ is always close to the background hydraulic gradient. The stabilizing effect is then again small, as drawing water out of the storage system does not significantly affect hydraulic gradients and turbulent dissipation.

30
In summary, Nye's instability will occur if the conduit is channel-like and water storage in the system is sufficiently large, and if the system length L is big enough. Note that with the choices in (1g), the system is always 178 Pa m −1 Table 1. Parameter values common to all calculations except those in figure 11, for which Ψ0 = 1630 Pa m −1 (corresponding to a much steeper 17:100 slope), u b hr = 1.05 × 10 −07 m 2 s −1 and L = 5 km.
conduit is always a channel) and L = ∞ (so effective pressure does not alter the hydraulic gradient). This is in agreement with Ng (1998).

Stability boundaries
So far, we have only established abstract conditions under which steady drainage through the conduit is unstable and will therefore not persist, potentially leading to periodic outburst floods. We can go further and determine explicitly the regions of 5 parameter space in which this stability occurs. While we present our results here in dimensional form, note that it is possible to reduce the spatially extended and lumped models (1) and (2) to a four-dimensional parameter space by non-dimensionalizing them (sections 2.2 and 2.3 of the supplementary material). These parameters are dimensionless versions of the inflow rate q in , storage capacity V p , system length L and conduit cut-off size S 0 . Recall that S 0 is intended to have minimal impact on conduit evolution away from the glacier margin, and we are therefore interested in the limit of large S 0 . As a result, we restrict 10 ourselves to the three-dimensional parameter space spanned by (q in ,V p , L) and set S 0 = ∞ in the lumped model (2). The remaining fixed parameter values we have used are given in table 1. The low value of Ψ 0 stated corresponds approximately to a 1:50 surface slope.
It can be shown analytically that there is at most a finite range of values of q in for which the instability occurs for givenV p and L (section 3.1 of the supplementary material): when q in is too small, the conduit is cavity-like as opposed to channel-like, 15 so (5) is not satisfied. By contrast, when q in is too large, the stabilizing termV −1 p q ψ L −1 in (13) dominates. The intermediate range of inflow values q in that makes the system unstable is not guaranteed to exist: a sufficiently large reservoir sizeV p and a long flow path length L are required to prevent the stabilizing term from dominating. The range of unstable q in values also increases asV p and L increase. This is confirmed by direct numerical computations of the stability boundaries. These are the locations in parameter space 20 where the real part of the growth rate λ in (12)  10 km (red), 50 km (blue) and 250 km (black); the latter is not intended to be physically realistic but to reflect the limit of large L. The Hopf bifurcation at the stability boundary is supercritical where the curve is solid, subcritical where dashed (see section 3.3. The dot-dashed curves show the asymptotic stability boundaries (14) and (15). stability boundaries in the (q in ,V p )-plane for different values of L. Note thatV p is displayed in units of km 2 , normalizing by ρ w g; in other words, what is really plotted isV p /(ρ w g), the surface area of the lake (see section 2). Similarly, we will plot N in units of metres, normalizing by ρ w g: changes in N plotted then correspond directly to changes in water level in the lake.
In each case, there is a region of instability above some critical value ofV p , and for q in in some intermediate range, where largerV p is required and the unstable inflow range is shrunk for shorter flow path lengths L.

5
For large enoughV p and L, we can determine analytical expressions for the stability boundaries. The lower critical value of q in at which the drainage system first becomes unstable corresponds to the switch from a cavity-like to a channel-like conduit, and this occurs for large L at (see Schoof, 2010, and section 3.1 of the supplementary material) The upper critical value of q in at which the system stabilizes again can similarly computed in the limit of large L by omitting 10 the cavity opening term v 0 and reducing conduit evolution (2a) to a balance between the first (dissipation-driven melting) and third (creep-closure) terms on the right hand side. This yields These limiting forms are also shown in figure 1, where we see that they are most accurate for large L as expected.

Hopf bifurcations and limit cycles
The analysis above has been purely linear, identifying parameter regimes in which steady drainage is unstable. What the 5 analysis does not allow is to say is what happens when the perturbations grow in size to the point where the linearization fails.
A key aspect of many outburst floods is that they are a recurring phenomenon; in terms of our reduced model (2) with two dynamical degrees of freedom and steady forcing, such recurrence must correspond to a stable periodic oscillation in the absence of time-dependent forcing (see also Kingslake, 2013, for time-dependent forcing that leads to chaotic solutions). As was underlined by Ng (1998) and Fowler (1999), the existence of such a limit cycle does not simply follow from the instability 10 itself: we need to ensure that the evolution away from the steady state leads to bounded growth of the instability once it reaches a finite amplitude (as opposed to the infinitesimal perturbation assumed by the linearization in section 3.1). This cannot be done in the confines of the linearization of section 3.1 alone.
It is straightforward to demonstrate bounded growth in our model computationally. Figure 2 shows a sample calculation of a periodic solution, with the classical attributes of a outburst flood cycle: effective pressure slowly decreases during the interval 15 between outburst floods, when conduit size is small. This is simply the lake refilling. As N approaches its minimum, the conduit size S starts to grow rapidly, initiating the outburst flood: effective pressure is no longer large enough to keep conduit size small, and enough water can flow to start enlarging the conduit in the runaway growth envisioned by Nye (1976). The lake then drains, rapidly increasing N . Once effective pressure gets large enough, creep closure becomes dominant, causing S to shrink again, and the cycle repeats. By contrast with Ng (1998) and Fowler (1999), key to this periodic behaviour is that the 20 conduit cannot become arbitrarily small between floods, since it is kept open by ice flow over bed roughness.
In fact, S in our model actually starts to increase immediately after flood termination, as the refilling of the reservoir leads to decreasing effective pressure allowing the now cavity-like conduit to grow. As explained above, water flow through them will eventually lead to re-initiation of dissipation-driven melting and enlargement of the conduit. By contrast, once the amplitude of the floods has become large enough, conduit size in the pure channel model of (Ng, 1998) and (Fowler, 1999) keeps shrinking 25 during the refilling phase until the effective pressure changes sign to negative, and this underpins the ever-increasing flood amplitudes in their model.
An alternative and possibly preferable way to visualize the periodic solution is in a phase plane, plotting S against N as the system evolves; a periodic solution then corresponds to a closed orbit. Figure 3 shows the phase plane equivalent of (2), with the dashed lines corresponding to nullclines (the curves on which eitherṠ = 0 orṄ = 0. During the refilling phase, the point 30 (N (t), S(t)) closely follows the S-nullcline, with S steadily increasing as explained above. periodic solutions and trace how they change under parameter changes using an arc length continuation method (section 3.4 of the supplementary material), and to determine simultaneously whether the periodic solutions are stable: this is not guaranteed, and some periodic solutions are never attained by forward integration of (2) because small perturbations will cause the system to evolve away from them.
As in section 3.2, we focus on how solutions change as the water input to the reservoir q in is varied. Figure 4 shows how the As q in crosses the lower critical value at which instability first occurs (associated with the change in the steady state conduit 15 from cavity-to channel-like behaviour), a limit cycle of small amplitude is formed, with that amplitude growing progressively as q in increases. This is consistent with a supercritical Hopf bifurcation: in general, the change from stability to instability in our model corresponds to the eigenvalue λ in (12) attaining the purely imaginary value ±i √ a 2 , and a standard result from the theory of dynamical systems is that a local, small-amplitude periodic solution exists near the bifurcation (Wiggins, 2003). What is non-trivial to determine a priori is whether that periodic oscillation exists on the stable or unstable side of the bifurcation (in 20 this, for values of q in less than or larger than the critical value, respectively). The stability of the periodic solution complements that of the corresponding steady state: if the periodic solution appears on the side of the bifurcation where the steady state has become unstable (a case termed a supercritical Hopf bifurcation), then the periodic solution is stable. This is the case for the lower critical value of q in in all panels of figure 4.
As the upper critical value of q in is approached, the stability analysis of section 3.1 predicts that the steady state returns to 25 stability. As this happens, we see that the amplitude of oscillations only shrinks continuously back to zero for the smallest value ofV p considered (in which case the upper critical value corresponds to another supercritical Hopf bifurcation). In panels a-c, the amplitude of stable periodic solutions continues to increase until values of q in near the upper critical value are reached. The amplitude then decreases slightly with a further increases in q in , before the periodic solution ceases to exist at a third critical value of q in that is larger than the threshold at which the steady state has become stable again, as shown in more detal in the its supercritical cousin that we encountered above). Panel d (see inset) represents an exotic exception to this, where there are two stable periodic solutions for a small interval of q in .
The nature of the Hopf bifurcations can be determined more quickly than by the numerical continuation method used above, using a weakly nonlinear stability analysis in the vicinity of the stability boundary (see section 3.3 of the supplementary material). This allows us to map out in figure 1 where that boundary corresponds to a supercritical Hopf bifurcation, in whose 5 vicinity a small-amplitude stable periodic solution emerges (stability boundary indicated as a solid curve), and where the Hopf bifurcation is subcritical and the stable periodic solution has a large amplitude. As in the few samples shown in figure 4, we see that the lower critical value of q in is always supercritical, while the upper critical value is generally subcritical, except at low V p .
In practical terms, figure 4 shows that for every unstable steady state, there is a corresponding stable periodic solution, so 10 the system evolves away from the steady state into an oscillation of finite amplitude as shown in figures 2 and 3. In practice, we may wish to know not only how the amplitude of oscillations (that is, of variations in N during the flood cycle) depends on water supply to the reservoir, but also what the recurrence period of the floods is. Figure 5 shows the period of the stable periodic solutions shown in figure 4 as a function of water supply rate q in , using the same colouring scheme to distinguish different values of V p . Perhaps unsurprisingly, large inflow rates q in correspond to more rapid flood cycles, and large reservoir 15 volumes correspond to floods repeating more slowly, albeit with a larger amplitude.
There is a significant caveat here: in many cases, the limit cycle solutions we have computed predict that N becomes negative during the cycle of reservoir filling and draining. This is not something that the model (2) It is important to know where in parameter space the outburst flood mechanism will change away from conduit growth due to cavities becoming channel-like at the end of the refilling phase, and instead involve water separating ice from bed at vanishing effective pressure. The boundary between these two regimes should correspond to the location in parameter space 5 where the minimum effective pressure during the flood cycle is zero. Figure 6 shows that parameter regime boundary, computed numerically (section 3.4 of the supplementary material) and superimposed on the stability boundary plot of figure 1. Clearly, sheet-flow-initiated floods (in which the glacier starts to float at the start of the outburst flood) are favoured at high water input rates and smaller reservoir volumes, where the conduit is less able to adjust to the rapid refilling of the reservoir.

Asymptotic solutions 10
We can also address limit cycle solutions through asymptotic methods in some parametric limits in our model. The most relevant limit for 'real' glacier-dammed lakes is likely to be that of a relatively large reservoir that is filled relatively slowly, but where water supply is not so small as to allow the conduit to be cavity-like in steady state (in which case the reservoir would be drained steadily, without a flood cycle): This is the case of large V p and moderate but not small q in , and is described in appendix A and, in detail, in section 4 of the supplementary material. This region of parameter space lies in the upper left of the unstable region of figure 6, between the near-vertical solid black line and the solid red curve.
In brief, the asymptotic solution confirms that there is a periodic flood cycle that the system very quickly settles into, and 5 that the flood cycle consists of three distinct stages. During the main flood stage, the evolution of the conduit is rapid and dominated by dissipation-driven melt c 1 qΨ and creep closure v c (S, N ) in (2a), and water input q in to the lake is much smaller than outflow q. This is followed by a long refilling phase in which the conduit has shrunk dramatically and therefore behaves as a cavity. Its size is dictated by a quasi-equilibrium in which the cavity opening term v o (S) balances the creep closure term v c (S, N ), leading to a slow opening of the conduit as the reservoir fills and N consequently dropped. Outflow q from the 10 reservoir during this phase is insignificant, and the mass balance of the lake is dominated by inflow q in . The refilling phase is terminated by a flood initiation phase whose length is intermediate between the main refilling phase and the rapid flood phase.
During this initiation phase, the reservoir is still filling, but the conduit has enlarged sufficiently that dissipation-driven wall melting c 1 qΨ starts to be significant.
Qualitatively, this solution is illustrated by the limit cycle in figures 2 and 3 and panel a of figure 8. The asymptotic solution 15 only becomes quantitatively accurate when V p is large enough to be physically unrealistic for real glacier-dammed lakes (section 4.4 of the supplementary material); this is because the relatively large exponent n = 3 in Glen's law makes the creep closure term quite sensitive to changes in N when N is small, and the initiation phase of the floods is affected significantly by this.
One of the predictions of the asymptotic solution is that the amplitude of effective pressure oscillations should be insensitive to refilling rate q in , except close to the Hopf bifurcation at which oscillations are initiated. As a result, the flood recurrence period (which is essentially the time taken to refill the reservoir in the limit where reservoir drainage is fast) should simply be inversely proportional to q in . Specifically, the result states that The same asymptotic solution also provides an estimate for the inflow rate q in at which zero effective pressure is first reached during the flood cycle, and our model ceases to be physically realistic as described above. The estimate is given by an analysis of the flood initiation phase (section 4.4 of the supplementary material) as This is superimposed on figure 6. The formula above is in general an underestimate of the value of q in at which zero effective pressures are reached, but clearly gives the correct scaling of how that value relates to V p .
We are also able to construct a second asymptotic solution for the opposite case of a large water supply rate q in (section 5 5 of the supplementary material). This predicts rapid oscillations in the reservoir level (or effective pressure) whose amplitude slowly evolves to a steady value. These rapid oscillations however invariably involve effective pressure changing from negative to positive: at leading order the growth and shrinkage of the conduit is controlled by the creep 'closure' term, which becomes a creep opening term during about half of the cycle as in panel b of figure 8, and dissipative melting is a higher order correction that ultimately dictates the amplitude that the rapid oscillations settle into over time. As discussed above, such solutions are 10 clearly not physically viable, and the model must be amended to take account of ice-bed separation; the asymptotic solution here serves merely to confirm that negative effective pressures are a robust prediction of our model for these large inflow rates, when the possibility of ice-bed separation is not taken into account.
4 The spatially extended model The analysis described above provides a comprehensive picture of the qualitative attributes of the lumped model (2). Here, we 15 consider how well that lumped model represents the behaviour of the its more complete, spatially-extended counterpart (1).
We begin by recreating the stability boundary diagrams of figure 1. The method by which the latter were computed explicitly (sections 3.1 and 3.2 of the supplementary material) cannot be applied directly to the extended model (1), since we have no closed-form solution to a linearized version of the model analogous to (12). Instead, we grid the (q in ,V p ) parameter space used in figure 1 finely. For each (q in ,V p ) pair, we discretize (1), compute steady states and perform a linear stability analysis 20 numerically as described in Schoof et al. (2014); this allows us to delineate regions of instability, although in a less sophisticated way.
Results are shown in figure 7. It is clear that the lumped model consistently underestimates the range of parameter values over which steady drainage is unstable. The onset of instability at the transition from cavity-to channel-like conduit behaviour appears to remain robust except in the case of small domain lengths L (panel c), for which the system appears to be unstable for 25 combinations of small storage capacitiesV p and inflow rates q in . Where the lumped model typically underestimates instability is at low storage capacities, and at large flow path lengths L. The latter is particularly significant since we have previously attributed stabilization to the effect incipient reservoir drainage reducing the hydraulic gradient along the flow path and therefore reducing flow through the conduit. While stabilization at large water input rates q in does eventually occur, this occurs at values that can be several orders of magnitude larger than those predicted by the lumped model, especially for the case of a large flow path lengths L. Furthermore, for a given storage capacityV p , the lumped model predicts that there is a single interval of inflow rate values q in over which instability occurs. The spatially extended model by contrast has two or more such intervals for most values ofV p , with a narrow region of stability between (the diagonal bands of solid black diamonds in panels a and b of figure 7).

5
To understand this discrepancy better. we have solved for the nonlinear evolution of the draiange system as described by the spatially extended model (1) for the parameter values indicated by magenta circles in panel a of figure 7.
The there smallest values of q in all correspond to unstable steady states in the lumped model (2), and we can compare spatially extended and lumped solutions. Figure 8 shows results. While the spatially extended solution is an infinite-dimensional dynamical system and strictly speaking cannot be visualized using a phase plane, we can overlay plots of conduit size S(0, t) 10 at the upstream end of the conduit against effective pressure N (0, t) at the same location onto the S-N phase plot for the lumped model. This is shown in the top row of panels. In all three cases, the extended models settles into a limit cycle, and for small and intermediate q in , we find good agreement between extended and lumped models; this only breaks down as the upper critical value of q in at which the lumped model stabilizes is approached. The lower row of panels in figure 8 shows snapshots of N (x, t) against x for the correspnding limit cycle shown in the top panel of the same column. For the smaller two values of q in , we see that pressure gradients ∂N/∂x are in general moderate away from the glacier terminus and therefore do not contribute significantly to the hydraulic gradient Ψ: this was the basis for the reduction of the spatially extended model (1) to the lumped form (2), and explains the good agreement. For larger q in , effective pressure N (x, t) along the conduit starts to develop wave-like structures, causing the approximation of negligible 5 pressure gradients to break down, and the discrepancy between lumped and extended model grows.
The extended model remains unstable beyond the stability boundary model of the lumped model, but our numerical solutions no longer support the conclusion that the system necessarily settles into a limit cycle. Figure 9 shows the evolution of the system for V p = 408 m 3 Pa −1 s −1 (a lake with a surface area of 4 km 2 ) and q in = 2.18×10 3 m 3 s −1 as in the magenta marker labelled 'D' in figure 7a (the fact that this is an inflow rate of biblical proportions is not as relevant as it may at first seem, as we path over time, but without the growth becoming bounded. Panels c and d show snapshots over the last two oscillations in the computation. Note that the vertical scale in panel c is logarithmic. What we see is that conduit size develops an aneurysmlike, massively enlarged feature near the terminus, blocked by a narrow constriction that requires an extremely large pressure gradient to overcome (note that panel c uses a logarithmic verticla scale) . The simulation eventually terminates when the solver fails to converge, and it is unclear whether this is merely a computational problem or indicates that the true continuum solution 5 becomes pathological or ceases to exist.
For even larger q in , a narrow band of stable values is passed and a different instability ensues as shown in figure 10. This instability appears to lead to bounded growth, again marked by rapid oscillations whose amplitude slowly grows, and significant pressure gradients along the flow path.
Once again, a key observation is that the spatially extended model predicts negative effective pressures long before any 10 kind of limit cycle is reached for all but one of the solutions shown, and that the discrepancies between lumped and spatially extended models occur only where this is the case. Moreover, they occur for water supply rates to the reservoir that far exceed values that would be plausible for typical glacier-dammed lake systems: the lumped model appears to be robust for the latter, including its prediction of where overpressurization of the drainage system with negative N first occurs.
Nonetheless, as we will discuss shortly, this does not render some of the more exotic instabilities shown in figures 8-10 15 irrelevant: The parameter values we have chosen here were deliberately chosen to reflect typically large reservoirs dammed by low-angled and large glaciers. Much smaller reservoirs in shorter more steeply-angled glaciers are liable to give rise to some of the instabilities that lie outside the reach of standard glacier-dammed lakes.

Glacier-dammed lakes and small reservoirs
We have shown that a drainage system that switches between cavity-like and channel-like behaviour spontaneously is capable 5 of supporting periodic outburst floods from a glacier-dammed reservoir. At issue here is the initiation of the flood: as discussed by Ng (1998) and Fowler (1999), if the conduit consists purely of a Röthlisberger-type channel, kept open only by melting due to heat dissipated in the turbulent flow of water through the channel, then the initiation of floods can be delayed more and more, without a limit cycle emerging. Both of these authors proposed that englacially-routed water supplied to the channel should keep it open at a minimum size, and the migration of the flow divide within the channel as the lake fills then becomes 10 the control on flood initiation. Our theory differs from theirs in the sense that no englacial water supply is necessary and the lake can continuously leak water through a linked-cavity-type drainage system that spontaneously becomes 'channel-like' and undergoes runaway enlargement through dissipation-driven melting as the lake fills.
At the heart of the oscillatory behaviour of glacier-dammed lakes is exactly that instability of a channel, which prevents steady discharge of the water supplied to the reservoir (Nye, 1976;Ng, 1998;Fowler, 1999). As described in Schoof et al.

15
(2014) and section 3.1 above, this instability does not occur when the conduit draining the reservoir acts as a set of linked cavities: when these are capable of steadily draining the water supplied to the reservoir, no outburst floods occur.
By contrast, for typical glacier-dammed lakes, in which moderate inflows q in exceed the drainage capacity of a linked cavity system but the the large storage capacity V p ensures that the lake fills slowly compared with the time scale over which a basal conduit can evolve, the cycle of filling and draining follows the characteristic sequence of a long filling period in which outflow 20 from the lake is negligible, followed by a brief onset period in which the conduit starts to experience significant melt-driven enlargement but the lake level continues to rise, and an even shorter outburst flood in which inflow to the lake is dwarfed by drainage through the subglacial conduit. For a given lake size V p , larger inflow rates will slowly increase the amplitude of the lake level fluctuation during the flood cycle (figure 4), while significantly shortening the length of the flood cycle (figure 5).
As in Ng (1998) andFowler (1999), our flood initiation mechanism can also be too slow to respond to water input and lead 25 to our models (1) and (2) predicting negative effective pressures in the reservoir: in that case, the glacier ought to float and flood initiation is likely to take the form a sheet flow that subsequently channelizes, as explored in Schoof et al. (2012) and Hewitt et al. (2012). This effect is not included in our models here, but we are at least able to give an asymptotically valid (in the limit of very large reservoir sizes V p ) criterion for the water supply rate at which the change in flood initiation between unstable conduit growth and partial glacier flotation occurs (equation (17)). Adapting our model to account for this alternative 30 flood initiation mechanism is the obvious next step to take.
We have also investigated a mechanism previously identified in Schoof et al. (2014), by which flow out of a reservoir can be stabilized when the storage capacity in the reservoir is relatively small, so that typical drainage rates (comparable to the inflow rate q in ) lead to adjustment of effective pressure in the lake much faster than the conduit can evolve due to wall melting. The rapid response of reservoir water levels can have a large enough effect on hydraulic gradients to stabilize the flow. This is true at least for a short enough drainage system. Importantly, this happens at water throughput rates that are completely unrealistic for typical, large glacier lakes (and even these rates are underestimated by a simplified, 'lumped' model as described in section 4).

5
The reason why this stabilization mechanism is relevant is that it may explain why much smaller water reservoirs that are typically not recognized as lakes but nonetheless provide storage capacity (such as large moulins) do not invariably generate outburst-flood type behaviour: the stability diagram of figure 1 was generated for parameter values chosen to represent large lakes dammed by a long glacier with a small surface slope, but can easily be rescaled to represent other glacier geometries, reservoir storage capacities and water throughput. The relevant scaling is explored in section 2 of the supplementary material.

10
Here, we confine ourselves to a simple numerical demonstration: figure 11 recalculates the stability boundary for steeper, shorter valley glacier with a relatively small reservoir sizes. We use the same parameter values as in table 1, except for Ψ 0 , u b h r andL. The 'upper' stability boundary (the larger critical value of q in where the lumped model (2) becomes stable again) is much more accessible for reasonable flow rates, and some of the more exotic manifestations of instability explored in section 2.1 are more likely to be within reach of typical water input rates.

15
With this as motivation, we also briefly discuss a second instability. This diverges from the classical picture of a glacier outburst flood in that it does not require channel-like behaviour but can occur in a purely cavity-based drainage system. We analyze this instability briefly below, since the work already done in the present paper provides the necessary framework.

Distributed storage: a modification of Nye's instability
In section 2, we formulated a model for a single reservoir, for which Nye's instability predicts the onset of oscillations when 20 the conduit becomes channel-like. In Schoof et al. (2014, section 5), the case of spatially spread-out storage capacity, for instance in the form of numerous basal crevasses arrayed along the flow path, was considered in addition to a single reservoir, A numerical linear stability analysis was used to demonstrate that instability can also occur for cavity-like conduits (which are generally stable for a single reservoir) for the case of such distributed water storage. As in the case of a single reservoir, Schoof et al. (2012) find that there is a finite range of values q in for which instability then occurs, only that this extends to lower values of q in , where the conduit is cavity-like.

5
Here we build on the analysis in section 3.1 to shed further light on their resultts. The model used in Schoof et al. (2012) was where v(N ) is a decreasing function of N that describes storage of water per unit length of the conduit, and the rest of the notation used replicates that of the single-reservoir model (1). Below, we will denote by the equivalent of V p in (1). This model also requires boundary conditions; an obvious choice is a prescribed flux q = q in at the 15 inflow x = 0, with no actual reservoir there, and again N = 0 at the terminus x = L.
In order to take advantage of the work already done in section 3.1, we adopt an analytical approach to understanding the instabilities of this modified drainage model, complementing the numerical stability analysis of Schoof et al. (2014) and providing additional insight into the feedbacks involved.
With a finite domain size and the proposed boundary conditions above, a steady state solution to (18) will in general have 20 spatial structure as in Schoof et al. (2014), and a linearization around that steady state will lead to a boundary value problem with non-constant coefficients that is not amenable to a closed-form solution. To avoid this, we concentrate here on shorter length scales, and assume that we can use (18) with periodic boundary conditions at these scales. This yields a spatially uniform steady state solution defined implicitly by with a constraint onN through the prescribed volume of water in the system (which is preserved due to the assumed periodic boundary conditions). Note that this is closely analogous to (8) but simpler, as the hydraulic gradient here does not contain the gradient term retained in (8c).
Linearizing as N =N + N exp(ikx + λt), S =S + S exp(ikx + λt), we find The eigenvalue λ satisfies a quadratic of the form (10), with solution (12), but now with a 1 and a 2 given by The only term that differs intrinsically from (11) is a 2 , which now has an imaginary part. We see that the real part of a 2 is still invariable positive, so the real part of a 1 − 4a 2 is smaller than a 1 . However, we can no longer conclude that the real part of λ has the same sign as a 1 except in appropriate limits.

10
For short wavelengths (large k), the system is stable with (λ) < 0: the two limiting forms of the two eigenvalues are and the assumptions about q ψ , q S , v o,S and v c,S in (6)-(7) ensure that both of these expressions are negative. Limited storage over short length scales prevents large variations in water discharge and prevents short-wavelength disturbances in conduit size and effective pressure from growing.
Conversely, at sufficiently large wavelengths (small k), we have a 1 ∼ (c 1 q S + v o,S − v c,S ), a 2 ∼ O(k), and a channel-like conduit will be unstable with one eigenvalue behaving as λ ∼ (c 1 q S + v o,S − v c,S ), while a cavity-like conduit is stable at 15 long wavelengths. This is simply Nye's instability at work. Clearly, then, there is a critical wavelength at which a channellike conduit with distributed water storage will become unstable (since it is stable at short wavelengths and unstable at long wavelengths), and the critical wavelength is determined by the various parameters in the model. How large that wavelength is matters, because we have assumed in setting up our stability analysis that we are looking at only a relatively small part of the domain, with periodic boundary conditions. If the critical wavelength approaches the full domain size L, our simplified 20 analysis breaks down, and the instability is no longer guaranteed to occur. This is consistent with figure 7 of Schoof et al. (2014), where channel-like conduits eventually become stable for large enough water supply rates.
Notably, the critical wavelength for instability of channel-like conduits will depend on storage capacityv p . The more storage capacity there is, the smaller all the terms containing k are in (22) for a given wavenumber k, and the more likely the first term in the definition of a 1 is to dominate the eigenvalue λ, leading to instability for a channel-like system.

25
To increase the size of the potentially stabilizing terms therefore requires larger wavenumbers k, and therefore a larger range of short wavelengths is likely to be unstable. Although we are not able to address a system of finite length directly with this approach, our results indicate that systems of limited size may remain stable if storage capacity is sufficiently limited, and this is confirmed by Schoof et al. (2014).
that propagates, in this case downstream as the imaginary part of λ is then negative. It is also of the same size as the real part, so propagation is not slow. This unstable wave is not the result of Nye's instability as the conduit is cavity-like. Instead, it is the result of an interaction between the dependence of the conduit closing rate on effective pressure and the dependence of water drainage (which affects effective pressure through water storage) on conduit size: the dominant balance that underpins it is where the perturbed conduit evolution does not include the melting term c 1 qΨ at all. Equations (24) can be combined into the single 'diffusion' equationv p ∂ 2 N /∂t 2 − q s v −1 c,N ∂N /∂x ∼ 0, except that the roles of time t and space x have been reversed. We can identify the positive feedback causing growth as a the result of a phase lag, where S leads N in phase and thus ensures that ∂S /∂x is positive where N has a maximum, therefore ensuring that ∂N /∂t is also positive.
Naturally, this sketch is incomplete, as the 'diffusion' problem with the roles of x and t reversed is not well-posed. The 5 growth rate in (23) grows unboundedly as wavelength k −1 approaches zero, which is a clear sign something is amiss. The discussion above was merely designed to identify a positive feedback that can drive growth; a negative feedback is necessary to ensure there is a fastest growing wavelength and that short wavelengths are damped as already discussed.

Conclusions
As demonstrated previously using a much more restricted sweep of parameter space in Schoof et al. (2014), we have shown that 10 drainage systems capable of switching spontaneously between channel-and cavity-like behaviour are stable in the presence of a localized water reservoir at low and high water throughput, with an unstable intermediate range of water fluxes. In that range, spontaneous oscillations in reservoir level will occur, driven by Nye's (1976) instability mechanism driving outburst floods.
These outburst floods turn out to be regular, periodic oscillations in water level at least at low-to-moderate water input, where a simplified, 'lumped' model of reservoir drainage generally reproduces the results of a more sophisticated, spatially-extended 15 drainage model. At high water throughput rates, our results are more equivocal as to the emergence of such limit cycles in the model used; in any case, the model necessarily breaks down physically (as opposed to mathematically) because it predicts that the oscillations will invariably reach negative effective pressures and therefore flotation of the ice dam at such large throughput rates.
It is worth pointing out that part of our focus has been on the emergence of limit cycle solutions in order to identify how 20 the flood initiation mechanism can prevent flood magnitude from increasing progressively from cycle to cycle as observed in Ng (1998) andFowler (1999), and to present an alternative mechanism for suppressing the continued growth in amplitude from that considered by Fowler (1999), whose lake is necessarily 'sealed' between floods. That mechanism is the ability of a cavity-like drainage system that remains during the reservoir recharge phase to switch to a channlized drainage mode.
The fact that we illustrate this by showing evolution towards a limit cycle is not at odds with the chaotic behaviour observed 25 by Kingslake (2015) using a variant of Fowler's (1999) model: the chaotic behaviour is intrinsically the result of a time-varying water input to the reservoir, which we have not studied in this paper. It is entirely possible, and in fact likely, that our model will also behave chaotically under such time-varying water input rates q in (t). The point of demonstating that limit cycles emerge was rather to underline that the growth of Nye's (1976) instability is bounded in our model without having to appeal to additional physics.

30
The lower cut-off to the drainage instability that leads to outburst floods corresponds to the drainage system switching to a cavity-like state under steady flow conditions when water input to the reservoir can be drained steadily by those cavities.
The mechanism for the cut-off at high water throughput rates is harder to identify. The lumped model predicts that the high sensitivity of water level in the reservoir to the evolution of the draining conduit will induce water pressure gradients that reduce flow as the conduit grows, and therefore suppress its enlargement due to heat dissipation. The threshold at which this happens is however significantly underestimated by lumped model relatively to its spatially-extended counterpart, which develops more wave-like instabilities at higher water throughput rates.
In closing, we have also investigated how wave-like instabilities can occur when the water reservoir is not localized but 5 spread out or 'distributed' along the flow path (for instance, in the form of many small reservoirs like basal crevasses). This type of instability was first observed in Schoof et al. (2014), and persists even where water throughput is insufficient to lead to channel-like behaviour. Adapting the stability analysis performed on a model with localized storage, we have shown that Nye's instability persists, but also that a second instability mechanism emerges, in which a phase shift between water storage and flux arises that causes water to accumulate in regions where effective pressure is already low.

10
Future work is likely to focus on capturing the role of overpressurization of the drainage system in initiating and mediating the instability driving outburst floods, since flood initiation at water pressures below flotation is confined to a relatively small part of parameter space, and the model predicts that reaching zero effective pressure and initiation by partial flotation of the ice dam is likely to be common.  . We assume that the exponents n and α satisfy n + 1 > α > 1, in which case the asymptotic solution we develop is valid when δ δ (α−1)(n+1)/(αn) 1 and ν 1.
The main flood phase is described by omitting terms of O(δ) and O( )

15
where matching with the initiation phase described later requires that (N, S) → (0, 0) as t → −∞. The transformation P = N/S allows the system (A2) to be re-cast in such a way as to demonstrate that the orbit into the fixed point (0, 0) is unique, so there is only one flood phase solution. The orbit terminates at a finite N =Ñ f as t → ∞; this then sets the amplitude of the lake level fluctuations that give the asymptotic formula for the flood cycle period (16).
The refilling phase is described by the rescalingÑ = N ,S = δ −1 S andt = (t − t f ), where t f is the time of the last flood. 20 At leading order, dÑ dt = −1; conduit size is quasi-steady and cavity-like, while lake level evolves purely because of inflow.