the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
An analysis of instabilities and limit cycles in glacierdammed 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 the flood, and how big can the lake get before that point is reached? Here we investigate how the spontaneous channelization of a linkedcavity drainage system can control the onset of floods. In agreement with previous work, we show that floods only occur in a band of water throughput rates in which steady reservoir drainage is unstable, and we identify stabilizing mechanisms that allow steady drainage of an icedammed reservoir. We also show how stable limit cycle solutions emerge from the instability and identify parameter regimes in which the resulting floods cause flotation of the ice dam. These floods are likely to be initiated by flotation rather than the unstable enlargement of a distributed drainage system.
 Article
(4815 KB)  Fulltext XML

Supplement
(889 KB)  BibTeX
 EndNote
Glacier lake outburst floods or “jökulhlaups” are a glacial hazard in many parts of the world (e.g. Björnsson, 1988; Clague et al., 2012). In addition, they provide a window into subglacial drainage systems: most outburst floods involve the opening and closing of a subglacial conduit, driven by melting of its walls through heat dissipation in the turbulent flow of water and by viscous creep closure of the ice. In the early stages of the outburst flood, wall melting dominates through a positive feedback in which conduit enlargement leads to faster water flow and therefore more dissipation of heat. Later in the flood, creep closure accelerates as the lake level drops, eventually causing the conduit to close again and terminating the flood. Alternatively, the flood can terminate because the lake has run dry. The conduit then becomes partially air filled and closes with minimal flow going through it.
While the mechanics of the main discharge phase of an outburst flood are relatively well understood (Nye, 1976; Clarke, 1982, 2003; Fowler and Ng, 1996; Ng, 1998; Fowler, 1999; Kingslake and Ng, 2013; Kingslake, 2013), the mechanisms that initiate the flood are still poorly understood, even though they dictate the level to which the lake is able to fill and therefore the magnitude and timing of the flood. In other words, the challenge is not to explain how a single flood progresses once started, but how it starts and, more to the point, what makes repeated floods occur cyclically, as a selfsustaining oscillation in the drainage system (Fowler, 1999; Kingslake, 2015). The main discharge phase of the flood is generally only part of a larger cycle, which also comprises a slow recharge phase in which outflow from the lake is small or absent altogether, and the transition between the recharge and discharge phases, where both lake inflow and outflow are of comparable magnitude.
One possibility for flood initiation is that the lake simply fills to the level at which the ice dam starts to float, and a sheet flow emerges between ice and bed that subsequently channelizes (Flowers et al., 2004). This may occur irregularly in some lakes such as Grímsvötn due to exceptionally large inflow rates to the lake, for instance during volcanic eruptions, or as part of a repeating flood cycle in others (Bigelow et al., 2020). Some lakes are however known to initiate outburst floods before they reach that flotation level.
Motivated by the subglacial lake at Grímsvötn in Iceland, Ng (1998) and Fowler (1999) consider how outburst floods are initiated in a drainage system that consists of a Röthlisberger (R) channel (Röthlisberger, 1972). In the simplest version of their model, where the channel is fed purely by the lake, they find that the amplitude of floods grows from cycle to cycle and that water pressures exceeding ice overburden are eventually reached, meaning that floods should start through the flotation of the ice dam after a few cycles. At the heart of this behaviour is the fact that the R channel can shrink to progressively smaller sizes between successive floods, making it harder to reinitiate drainage.
Grímsvötn typically starts its flood when water levels are below flotation (Björnsson, 1988), without successive floods growing in amplitude, except when the flood results from a large increase in water input due to volcanic activity (Gudmundsson et al., 1997). To explain this behaviour, Fowler (1999) considers the effect of a water supply along the length of the R channel on its evolution. Such a water supply can maintain a minimum channel size even between outburst floods, with flow of water being directed partly into the lake and partly downglacier to the margin, provided the glacier has a geometrical “seal”. As the lake fills, the flow divide inside the channel migrates towards the lake, and the flood begins when the divide reaches the edge of the lake.
While this mechanism successfully explains how limit cycles (stable, periodic oscillations in lake level) can emerge in the model, it also predicts that no water can leave the lake between floods. Tracer experiments conducted at Salmon Glacier in Canada (Fisher, 1973) demonstrate that lakes can leak continuously throughout their recharge phase. Summit Lake, dammed by the Salmon Glacier, also has a history of flood initiation at lake levels below the flotation level (Post and Mayo, 1971). At Gornersee in Switzerland, observations of water lake levels and inferences of lake water balance based on meteorological measurements over two consecutive years also suggest leakage and flood initiation below flotation level during one year and flood initiation at or above flotation level during the previous year, potentially in the absence of any preflood leakage (Huss et al., 2007). Observations at Hidden Creek Lake (Anderson et al., 2003) by contrast suggest no leakage prior to flood initiation, though that conclusion is again based on water balance estimates rather than direct tracer experiments.
In this paper, we consider an alternative mechanism by which floods can initiate in a recurring fashion without the need for a sealed lake, but with continuous leakage throughout the floodingandrefilling cycle. We show that a conduit that is able to switch spontaneously between the behaviour typical of an R channel and the behaviour of a linkedcavity system (Kessler and Anderson, 2004; Schoof, 2010; Hewitt et al., 2012; Hewitt, 2013; Werder et al., 2013) can sustain limit cycles without leading to flotation of the glacier at flood initiation, and without requiring the flowdivide migration appealed to in Fowler (1999). Note that many elements of the model studied here were included in the earlier study of outburst flooding by Kessler and Anderson (2004), who were able to replicate leakage from a glacierdammed reservoir before the onset of the outburst flood but did not attempt to reproduce recurring floods. It is unclear if their model would have been able to do so since, in common with the model due to Fowler (1999), the main trunk of their drainage system was set up as a pure R channel.
The behaviour studied in the present paper was first documented in a model by Schoof et al. (2014). In this paper, we approach the problem from a more theoretical perspective, focusing on the following: first, we delineate when a “lake” or other storage reservoir can be drained steadily, identify conditions under which steady drainage becomes unstable, and identify how that instability leads to a limit cycle. Next, we investigate how the amplitude and period of floods depends on reservoir size and inflow rate. As a corollary, we determine at which point (in parameter space) the model predicts that flotation does occur and a different flood initiation mechanism is likely to take over. In closing, we also briefly comment on the difference in flood dynamics caused by distributing water storage along the flow path rather than having a single water reservoir (Schoof et al., 2014). We outline how the classical runaway melt mechanism for outburst floods first described by Nye (1976) can be replaced as a mechanism for selfsustaining water pressure oscillations by a spatial phase shift between conduit size and water pressure along the flow path.
The original motivation for the work in Schoof et al. (2014) was to demonstrate that behaviour akin to outburst floods can occur in systems with limited or even distributed water storage, manifesting itself in the form of unforced water pressure oscillations. Such limited storage could in principle result from moulins or any other vertical shaft that can fill progressively with water as water pressure rises. Viewed from that alternative perspective, the present paper analyses a drainage model that has become widely adopted in the study of subglacial hydrology (Werder et al., 2013). We identify instabilities that spontaneously lead to selfsustained 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.
We use a hierarchy of model versions, employing both a spatially extended, onedimensional drainage system model and a “lumped”, boxtype model that provides additional physical insight, as well as being the appropriate limit of the spatially extended model in the case of a long flow path. The paper is laid out as follows: in Sect. 2, we develop the spatially extended and lumped models. In Sect. 3.1 and 3.2, we identify instabilities in the lumped model. The evolution of that instability at finite amplitude is investigated in detail in Sect. 3.3 and 3.4, where we show that stable limit cycles emerge in the lumped model. In Sect. 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 glacierdammed 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 Sect. 5.1, and we briefly touch on the dynamics of systems with distributed water storage in Sect. 5.2, relegating detail to the Supplement.
2.1 A continuum model for outburst floods
We use the onedimensional continuum model for drainage through a single conduit stated 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 et al., 2012). Note that we use the term “conduit” here to refer to a generic drainage element that can evolve dynamically to behave as an R channel, in which dissipation and creep closure are the dominant mechanisms by which channel size changes, or as a cavity, in which opening due to sliding over bedrock and creep closure are dominant (Kessler and Anderson, 2004; Schoof, 2010).
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 $\mathrm{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 described in further detail below, 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 hydraulic 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 that occupied by the conduit). Typically (including in Eq. 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 dissipationdriven wall melting at a rate c_{1}qΨ, opening due to ice sliding over bed roughness at rate ${u}_{\mathrm{b}}{h}_{\mathrm{r}}(\mathrm{1}S/{S}_{\mathrm{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}_{\mathrm{3}}{S}^{\mathit{\alpha}}\mathrm{\Psi}{}^{\mathrm{1}/\mathrm{2}}\mathrm{\Psi}$ through Eq. (1g). u_{b} is sliding velocity, h_{r} is bed roughness and S_{0} is a cutoff cavity size at which bed roughness is drowned out (Schoof et al., 2012). S_{0} is 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}. Equation (1a) also ignores the effect of water pressure on the melting point, which affects the fraction of the dissipated heat available for melting of the conduit walls (Röthlisberger, 1972; Werder, 2016). We have also neglected the effect of water storage in the conduits and of dissipationdriven melting in the mass balance relation (1b); the reason for this is that a simple scaling exercise shows that these terms are invariably small for a terrestrial glacial system (Supplement Sect. S2.2).
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). A computationally preferable assumption to N=0 may be to assume that there is an ice cliff at the glacier terminus, so that vanishing water pressure still corresponds to a finite N and the effect of dissipationdriven opening does not have to be offset by an unrealistic closing of the conduit through the sliding term, which becomes negative when S>S_{0}; our results are however not sensitive to this assumption.
The upstream boundary condition at x=0 reflects water storage in a reservoir at the head of the conduit: Eq. (1e) represents conservation of mass in the reservoir, with inflow rate q_{in}, with outflow at rate q(0,t) through the conduit modelled by Eqs. (1a)–(1d) and with water storage V in the reservoir a function of effective pressure at x=0 such that $\mathrm{d}V/\mathrm{d}N={V}_{p}$ (see e.g. Clarke, 2003). Note that we deliberately use “reservoir” rather than “lake” here since we will be interested not only in large lakes, but also in more modestsized reservoirs as in Schoof et al. (2014). We will treat V_{p} and q_{in} as constants in most of what follows, although realistic lake shapes may argue for particular forms of the function V_{p}(N) (Clarke, 2003).
For simplicity, we have neglected the possibility of a partially air filled conduit and the effects of overpressurization, where effective pressure becomes negative (Schoof et al., 2012; Hewitt et al., 2012). Overpressurization is particularly relevant to the initiation mechanism for subglacial floods where unstable enlargement of the conduit is not fast enough for the flood to occur before the lake reaches N=0 and overfills. A partially air filled conduit may form by contrast at the upstream end of the flood path if the lake fully drains at the end of the outburst floods, and it could form more permanently at the downstream end as shown in Hewitt et al. (2012). We defer consideration of both processes to future work.
There is an important simplification we have made, irrespective of the particular choice of V_{p}: our model effectively imposes an ice cliff at the conduit inlet x=0, with the cliff height 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. In that case, the upstream conduit end migrates as the lake fills or empties, and its location is related to lake volume (see Sect. S2.4 of the Supplement).
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 $\partial N/\partial x$ are large enough to overcome the effects of the seal to make the total hydraulic gradient Ψ positive in Eq. (1d). In the 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 Sect. S2.4 of the Supplement).
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 and seal over a short distance, which can be balanced by an equally large $\partial N/\partial x$ over the same short distance, leading to a nonzero 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 $\partial N/\partial x$ while the lake is filling. As described in the introduction, the onset of the outburst flood corresponds to the instant at which Ψ at 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 dissipationdriven melting.
By contrast, our model assumes that the reservoir always experiences some amount of leakage and requires no englacially supplied conduit: the water supply to the drainage system in our model is purely through the inflow term q_{in} into the reservoir. Leakage out of the reservoir is facilitated by the linkedcavity 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 meltdriven conduit enlargement as in Fowler (1999) and Ng (1998), but the initiation of that enlargement differs.
2.2 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 $\partial N/\partial x$ along the flow path can be treated as negligible (see also Fowler, 1999; Ng, 2000, and Sect. S2.2 and S2.3 of the Supplement) or approximated by a simple divided difference.
Under these assumptions, we can relate the flux q along the flow path (which is independent of position by Eq. 1b) purely to the conduit size S and hydraulic gradient Ψ_{0} at the head of the conduit, 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
which corresponds to a crude divided difference approximation of the actual gradient $\partial N/\partial x$ by $\left[N\right(L,t)N(\mathrm{0},t\left)\right]/L$.
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 finitedimensional dynamical systems (Strogatz, 1994; Wiggins, 2003), leading a number of semianalytical 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 centres on the idea that steady flow in an R 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 effective pressure that could stabilize the channel against growth.
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 model (2). If we assume a classical R channel in the sense of Röthlisberger (1972) and Nye (1976) and therefore set the cavity opening term v_{o} to zero in Eq. (2a), and use the remaining relations in Eq. (1g), the conduit evolution Eq. (2a) becomes
with a steadystate conduit size $\overline{S}$ given by ${c}_{\mathrm{2}}\overline{S}\overline{N}{}^{n\mathrm{1}}\overline{N}={c}_{\mathrm{1}}{c}_{\mathrm{3}}{\overline{S}}^{\mathit{\alpha}}\mathrm{\Psi}{}^{\mathrm{3}/\mathrm{2}}$. At fixed effective pressure N, and hence at fixed Ψ, this steady state is unstable: an increase S^{′} away from the steadystate size $\overline{S}$ will lead to a further growth in the conduit, as we have $(\overline{S}+{S}^{\prime}{)}^{\mathit{\alpha}}\approx {\overline{S}}^{\mathit{\alpha}}+\mathit{\alpha}{\overline{S}}^{\mathit{\alpha}\mathrm{1}}{S}^{\prime}$ and so
where, with α>1, the righthand 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, Eq. (2a) can be expected to lead to unstable growth of conduits away from an equilibrium state $\overline{S}$ if (see also pp. 5–6 of the Supplement to Schoof, 2010)
This is also precisely the condition that defines a conduit as being channellike 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.
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. 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 Sect. 2.2.
Note that v_{o}, v_{c} and q satisfy
with v_{o}(S)>0 bounded as S→0, ${v}_{\mathrm{c}}(S,\mathrm{0})=\mathrm{0}$ and ${v}_{\mathrm{c}}(\mathrm{0},N)=\mathrm{0}$; q also has the same sign as Ψ and satisfies $q(\mathrm{0},\mathrm{\Psi})=q(S,\mathrm{0})=\mathrm{0}$. These are the minimal assumptions we make on these functions, allowing us to generalize from the specific forms in Eq. (1g) in analysing 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 $(\overline{S},\overline{N})$ given implicitly by
where we assume positive lake inflow q_{in}>0. It can be shown (see Sect. S3 of the Supplement) that there is a unique steadystate solution with positive $\overline{N}$, $\overline{S}$ and $\overline{\mathrm{\Psi}}$. To establish whether the solution is stable, we can linearize around the steady state as
This yields the following eigenvalue problem
where
Setting the determinant of the matrix on the left to zero, we obtain a quadratic with solution
where
From our assumptions on the various functions involved, we see that a_{2}>0 (recall that ${v}_{\mathrm{o},S}\le \mathrm{0}$ from the inequality Eq. 6), while a_{1} can be either sign. Since a_{2}>0, we have ${a}_{\mathrm{1}}^{\mathrm{2}}\mathrm{4}{a}_{\mathrm{2}}<{a}_{\mathrm{1}}^{\mathrm{2}}$ and two possible types of solution: either (i) ${a}_{\mathrm{1}}^{\mathrm{2}}\mathrm{4}{a}_{\mathrm{2}}>\mathrm{0}$, in which case we have two real roots, both of which have the same sign as a_{1}, or (ii) ${a}_{\mathrm{1}}^{\mathrm{2}}\mathrm{4}{a}_{\mathrm{2}}<\mathrm{0}$, resulting in 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 lefthand side of Eq. (11) as the difference of two terms, a potentially destabilizing term ${c}_{\mathrm{1}}{q}_{S}\overline{\mathrm{\Psi}}+{v}_{\mathrm{o},S}{v}_{\mathrm{c},S}$ and a stabilizing term ${\overline{V}}_{p}^{\mathrm{1}}{q}_{\mathrm{\Psi}}{L}^{\mathrm{1}}$. The first term can be recognized as the growth rate sensitivity we previously identified as being at the heart of Nye's instability in Eq. (5), as well as an indicator of whether the steadystate conduit is channellike in the terminology of Schoof (2010).
The second term in the relation (11) is a stabilizing term that is inversely related to storage capacity. Its physical origin is the following: the sensitivity of conduit 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 the effective pressure is increased, the hydraulic gradient $\mathrm{\Psi}={\mathrm{\Psi}}_{\mathrm{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 conduit. 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 and on the length L of the flow path: a large storage capacity ${\overline{V}}_{p}$ or a long flow path L leads to a reduced stabilizing effect. In practice, we are likely to see stabilization for small storage elements, such as moulins and individual crevasses.
3.2 Stability boundaries
We have established that Nye's instability will occur if the conduit is channellike and water storage in the system is sufficiently large, as well as if the system length L is big enough. Note that, with the choices in Eq. (1g), the system is always unstable if v_{o}=0 (so the conduit is always a channel) and L=∞ (so effective pressure does not alter the hydraulic gradient). This is in agreement with Ng (1998).
We can go further and determine explicitly the regions of 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 fourdimensional parameter space by nondimensionalizing them (Sect. S2.2 and S2.3 of the Supplement). These parameters are dimensionless versions of the inflow rate q_{in}, storage capacity V_{p}, system length L and conduit cutoff 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 ourselves to the threedimensional parameter space spanned by $({q}_{\mathrm{in}},{\overline{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.
Figure 2 shows stability boundaries in the $({q}_{\mathrm{in}},{\overline{V}}_{p})$ plane for different values of L. These are the locations in parameter space where the real part of the growth rate λ in Eq. (9) is zero (see Sect. S3.2 of the Supplement for details); invariably, the eigenvalues λ then form a purely imaginary conjugate pair, and we have a socalled Hopf bifurcation (see Sect. 3.3 below). Note that ${\overline{V}}_{p}$ is displayed in units of square kilometres (km^{2}), normalizing by ρ_{w}g; in other words, what is really plotted is ${\overline{V}}_{p}/\left({\mathit{\rho}}_{\mathrm{w}}g\right)$, the surface area of the lake (see Sect. 2). Similarly, we will plot N in units of metres (m) in subsequent figures, normalizing by ρ_{w}g: changes in N plotted then correspond directly to changes in water level in the lake.
In each case shown in Fig. 2, there is a region of instability above some critical value of ${\overline{V}}_{p}$, and for q_{in} in some intermediate range, where larger ${\overline{V}}_{p}$ is required and the unstable inflow range is shrunk for shorter flow path lengths L. This is consistent with theoretical results (Sect. S3.1 of the Supplement) showing that the system is always stable for either sufficiently small or sufficiently large q_{in} but can be unstable at intermediate inflow rates.
For large enough ${\overline{V}}_{p}$ and L, we can determine approximate 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 cavitylike to a channellike conduit, and this occurs for large L at (see Schoof, 2010, and Sect. S3.1 of the Supplement)
Note that for finite h ${\overline{V}}_{p}$ and L, there is generally at least a small region of stability in which q_{in} is already large enough for the conduit to be channellike (the space between the dotted–dashed vertical line and the solid curves in Fig. 2). This region exists because the first term in Eq. (11) is positive but not large enough to overcome the second, stabilizing term here.
The upper critical value of q_{in} at which the system stabilizes again can be similarly computed in the limit of large L by omitting the cavity opening term v_{o} and reducing the conduit evolution Eq. (2a) to a balance between the first (dissipationdriven melting) and third (creepclosure) terms on the righthand side. This yields
These limiting forms are also shown in Fig. 2, where we see that they are most accurate for large L as expected.
3.3 Hopf bifurcations and limit cycles
The analysis above has been purely linear, identifying parameter regimes in which steady drainage is unstable. That linearization will eventually fail where unstable growth is predicted, and a nonlinear model is needed to predict what the instability grows into. A key aspect of many outburst floods is that they are a recurring phenomenon. For the reduced model (2) with two dynamical degrees of freedom and steady forcing, such a recurrence must correspond to a stable periodic oscillation in the absence of timedependent forcing (see also Kingslake, 2013, for timedependent 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 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.
It is straightforward to demonstrate bounded growth in our model computationally. Figure 3 shows a sample calculation of a periodic solution, with the classical attributes of an outburst flood cycle: effective pressure slowly decreases during the interval 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 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 almost immediately after flood termination, as the refilling of the reservoir leads to decreasing effective pressure allowing the now cavitylike conduit to grow. As explained above, water flow through them will eventually lead to reinitiation of dissipationdriven 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 during the refilling phase until the effective pressure changes sign to negative, and this underpins the everincreasing flood amplitudes in their model.
An alternative visualization of the physics involved is a phase plane, plotting S against N as the system evolves; a periodic solution then corresponds to a closed orbit. Figure 4 shows the phase plane equivalent of Fig. 3, with the dashed lines corresponding to nullclines (the curves on which either $\dot{S}=\mathrm{0}$ or $\dot{N}=\mathrm{0}$). During the refilling phase, the point (N(t),S(t)) closely follows the S nullcline, with S steadily increasing as explained above.
Figures 3 and 4 show only one example. It is possible to solve directly for such periodic solutions and trace how they change under parameter changes using an arc length continuation method (Sect. S3.4 of the Supplement), as well as to determine simultaneously whether the periodic solutions are stable: some periodic solutions are in fact never attained by forward integration of the model (2) because they are themselves unstable. Small perturbations will therefore cause the system to evolve away from them.
As in Sect. 3.2, we focus on how solutions change as the water input to the reservoir q_{in} is varied. Figure 5 shows how the amplitude of oscillations varies with water input for a fixed L=50 km and different storage capacity V_{p} (treated as independent of N here). In each panel, the coloured curve (black in Fig. 5e) shows the minimum and maximum value of N attained for a periodic solution at the corresponding value of q_{in} (i.e. the values of N where the corresponding orbit in the phase plane crosses the N nullcline). Solid portions of these coloured curves in Fig. 5 correspond to stable periodic solutions, while dotted portions are unstable periodic orbits. The black curves in Fig. 5 generally correspond to steadystate solutions, plotting N in the steady state against q_{in} (we use black for steady states and oscillatory solutions in Fig. 5e, but these are easily distinguished by comparing the plots in Fig. 5e with those in the remaining panels). Again, solid black lines are stable steady states and dashed black lines are unstable steady states. Each panel in Fig. 5 corresponds to a different storage capacity V_{p} and therefore represents a horizontal slice through Fig. 2; V_{p} decreases from Fig. 5a to Fig. 5d.
In all cases, an unstable steady state corresponds to a limit cycle. As q_{in} crosses the lower critical value at which instability first occurs (associated with the change in the steadystate conduit from cavity to channellike behaviour), a limit cycle of small amplitude is formed, with that amplitude growing progressively as q_{in} increases, and the emerging limit cycle is an oscillation around an unstable steady state. This is consistent with a supercritical Hopf bifurcation (Wiggins, 2003).
As the upper critical value of q_{in} is approached, the stability analysis of Sect. 3.1 predicts that the steady state returns to stability. As this happens, we see that the amplitude of oscillations only shrinks continuously back to zero for the smallest value of ${\overline{V}}_{p}$ considered (in which case the upper critical value corresponds to another supercritical Hopf bifurcation). In Fig. 5a–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 increase 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 detail in the inset in Fig. 5b. (In technical terms, this is known as a saddlenode bifurcation of the Poincaré map of the dynamical system 2; see Wiggins, 2003.) Where this abrupt disappearance of the stable periodic solution occurs, it is accompanied by an unstable periodic solution near the upper critical value of q_{in}, corresponding to a subcritical Hopf bifurcation. Figure 5d (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 at the stability boundaries in Fig. 2 can be determined more quickly than by the numerical continuation method used above, using a weakly nonlinear stability analysis (see Sect. S3.3 of the Supplement). This allows us to map out in Fig. 2 where that boundary corresponds to a supercritical Hopf bifurcation, in whose vicinity a smallamplitude stable periodic solution emerges (stability boundaries shown as a solid curve), and where the Hopf bifurcation is subcritical and the smallamplitude periodic solution is unstable (dashed stability boundary). As in the few samples shown in Fig. 5, 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}.
We can also compute the period of the drainage oscillations, that is, the recurrence interval of floods. Figure 6 shows the period of the stable periodic solutions shown in Fig. 5 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 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 models are designed for. Instead of the conduit expanding through creep, as occurs in the models (1) and (2) when N<0, we expect that the glacier should instead detach from the bed, and a sheet flow of water initiates the outburst flood in this case; prototypes of models for this behaviour can be found in Schoof et al. (2012), Hewitt et al. (2012) and Tsai and Rice (2010).
It is important to know where in parameter space the outburst flood mechanism would in reality change away from conduit growth due to cavities becoming channellike 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 where the minimum effective pressure during the flood cycle is zero. Figure 7 shows that parameter regime boundary, computed numerically (Sect. S3.4 of the Supplement) and superimposed on the stability boundary plot of Fig. 2. Clearly, sheetflowinitiated 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.
3.4 Asymptotic solutions
We can also address limit cycle solutions through asymptotic methods in some parametric limits in our model. The most relevant limit for real glacierdammed 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 cavitylike 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 Sect. S4 of the Supplement. This region of parameter space lies in the upper left of the unstable region of Fig. 7, between the nearvertical 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 that the flood cycle consists of three distinct stages. During the main flood stage, the evolution of the conduit is rapid and dominated by dissipationdriven melt c_{1}qΨ and creep closure v_{c}(S,N) in Eq. (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 quasiequilibrium 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 drops. Outflow q from the 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 dissipationdriven wall melting c_{1}qΨ starts to be significant.
Qualitatively, this solution is illustrated by the limit cycle in Figs. 3 and 4 and Fig. 9a. The asymptotic solution only becomes quantitatively accurate when V_{p} is large enough to be physically unrealistic for real glacierdammed lakes (Sect. S4.4 of the Supplement); 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
where ${\stackrel{\mathrm{\u0303}}{N}}_{f}$ is a dimensionless constant, with a value of 1.44 for the parameters of α and n chosen here, in the limit of a large flow path length L (Sect. S4.1 of the Supplement). This asymptotic formula is overlaid onto the numerically computed periods in Fig. 6; it should be clear that the asymptotic formula performs poorly for the relatively moderate values of V_{p} used here. This should not be a surprise since Fig. 5 demonstrates that, for the same values of V_{p}, the amplitude of oscillations is in fact sensitive to the inflow rate q_{in}.
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 (Sect. 4.4 of the Supplement) as
where γ_{c}≈0.25 for the values of α and n chosen here. This critical inflow rate as a function of V_{p} is superimposed on Fig. 7. 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} (Sect. S5 of the Supplement). This solution is valid towards the righthand edge of the grey region in Fig. 7 and 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. As discussed above, such solutions are clearly not physically viable, and the model must be amended to take account of ice–bed separation.
The analysis above provides a comprehensive picture of the lumped model (2). Here, we consider how well that lumped model represents the behaviour of its more complete, spatially extended counterpart (1). We begin by recreating the stability boundary diagrams of Fig. 2. The method by which the latter were computed explicitly (Sect. S3.1 and S3.2 of the Supplement) cannot be applied directly to the extended model (1), since we have no closedform solution to a linearized version of the model analogous to Eq. (9). Instead, we grid the $({q}_{\mathrm{in}},{\overline{V}}_{p})$ parameter space used in Fig. 2 finely. For each $({q}_{\mathrm{in}},{\overline{V}}_{p})$ pair, we discretize the model (1), compute steady states and perform a linear stability analysis numerically as described in Schoof et al. (2014, Appendix B); this allows us to delineate regions of instability, although in a less sophisticated way. We adapt that method for the transient nonlinear calculations in Figs. 9–11 by solving Eqs. (B1)–(B2) of the same Appendix in Schoof et al. (2014) using a backward Euler step.
Results are shown in Fig. 8. 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 channellike conduit behaviour remains robust except for small domain lengths L (Fig. 8c), where the system appears to be unstable for combinations of small storage capacities ${\overline{V}}_{p}$ and inflow rates q_{in}. The lumped model typically underestimates instability 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 of 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 happens at values that can be several orders of magnitude larger than predicted by the lumped model, especially for large flow path lengths L. Furthermore, for a given storage capacity ${\overline{V}}_{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 of ${\overline{V}}_{p}$, with a narrow region of stability between (the diagonal bands of solid black diamonds in Fig. 8a and b).
To understand this discrepancy better. we have solved for the nonlinear evolution of the drainage system as described by the spatially extended model (1) for the parameter values indicated by magenta circles in Fig. 8a.
The three 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 9 shows results. While the spatially extended solution is an infinitedimensional dynamical system and strictly speaking cannot be visualized using a phase plane, we can overlay plots of conduit size S(0,t) 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 model 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 we approach the upper critical value of q_{in} at which the lumped model stabilizes.
The lower row of panels in Fig. 9 shows snapshots of N(x,t) against x for the corresponding limit cycle shown in the top panel of the same column. For the smaller two values of q_{in}, we see that pressure gradients $\partial N/\partial x$ are 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 wavelike structures, causing the approximation of negligible pressure gradients to break down, and the discrepancy between lumped and extended model grows.
The extended model remains unstable beyond the stability boundary of the lumped model, but our numerical solutions no longer support the conclusion that the system necessarily settles into a limit cycle. Figure 10 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}_{\mathrm{in}}=\mathrm{2.18}\times {\mathrm{10}}^{\mathrm{3}}$ m^{3} s^{−1} as in the magenta marker labelled “D” in Fig. 8a (the fact that this is an inflow rate of biblical proportions is not as relevant as it may at first seem, as we discuss in Sect. 5 below). Figure 10a and b show the growth of rapid oscillations in the maxima of N and S along the flow path over time, but without the growth becoming bounded. Figure 10c and d show snapshots over the last two oscillations in the computation. Note that the vertical scale in Fig. 10c 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. 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 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 Fig. 11. This instability appears to lead to bounded growth, again marked by rapid oscillations whose amplitude first grows and then saturates, 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 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 glacierdammed lake systems: the lumped model appears to be robust for the latter. This includes the prediction by the lumped system of where overpressurization of the drainage system (N<0 in the limit cycle) first occurs.
This does not render some of the more exotic instabilities shown in Figs. 9–11 irrelevant: the parameter values we have chosen here were deliberately chosen to reflect typically large reservoirs dammed by lowangled and large glaciers. Much smaller reservoirs in shorter, more steeply angled glaciers are liable to give rise to some of these exotic instabilities, even though they lie outside the reach of standard glacierdammed lakes.
5.1 Glacierdammed lakes and small reservoirs
We have shown that a drainage system that switches between cavitylike and channellike behaviour spontaneously is capable of supporting periodic outburst floods from a glacierdammed reservoir. At issue here is the initiation of the flood: as discussed by Ng (1998) and Fowler (1999), if the conduit behaves purely as a Röthlisbergertype 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 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 linkedcavitytype drainage system that spontaneously becomes channellike and undergoes runaway enlargement through dissipationdriven melting as the lake fills.
At the heart of the oscillatory behaviour of glacierdammed lakes is that runaway growth, or instability, of a channellike conduit, which prevents steady discharge of the water supplied to the reservoir (Nye, 1976; Ng, 1998; Fowler, 1999). As described in Schoof et al. (2014) and Sect. 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.
For typical glacierdammed lakes, moderate inflows q_{in} exceed the drainage capacity of a linkedcavity system, but the large storage capacity V_{p} ensures that the lake fills slowly compared with the timescale over which a basal conduit can evolve. The cycle of filling and draining then follows the characteristic sequence of a long filling period in which outflow from the lake is negligible, followed by a brief onset period in which the conduit starts to experience significant meltdriven enlargement but the lake level continues to rise, as well as 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 (Fig. 5), while significantly shortening the length of the flood cycle (Fig. 6).
As in Ng (1998) and Fowler (1999), our flood initiation mechanism can also be too slow to respond to water input and lead to our models (1) and (2) predicting negative effective pressures at the upstream end of the conduit: 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 Flowers et al. (2004), 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 approximate expression (valid in the limit of very large reservoir sizes V_{p}) for the water supply q_{in} at which the change in flood initiation mechanism occurs (switching between unstable conduit growth and partial glacier flotation and a sheet flow); this is Eq. (15). Adapting our model to account for this alternative flood initiation mechanism is the obvious next step to take. One straightforward adaptation of the lumped model to this effect, inspired by the model in Hewitt et al. (2012), is the following: if we suppose that lake level cannot rise beyond the point at which N=0, then the opening of an ice–bed gap must be such as to permit outflow balancing inflow to the lake at that point, which forces us to alter the constitutive relation for q once N=0. We can similarly build the possibility of flood termination due to the lake running dry at an upper bound on effective pressure at $N={N}_{max}={\mathit{\rho}}_{\mathrm{i}}g\left(s\right(\mathrm{0})b(\mathrm{0}\left)\right)$ into the model. Once the lake is fully empty, discharge q in the conduit will be the lesser of inflow into the conduit q_{in} (corresponding to a partially filled conduit) and the discharge the conduit would carry if it were completely filled with water. The appropriately adapted lumped model (2) is
The equivalent for the spatially extended model (1) is the channel model in Hewitt et al. (2012), with a reservoir at the upstream end. Importantly, the stability of steadystate solutions as in Sects. 3.1 and 4 is unchanged by these adaptations unless the steady state corresponds to an empty lake basin and a partially filled conduit, in which case the conduit is presumably stable since discharge in the conduit does not increase when its size is increased. The upper and lower bounds on flux introduced in the model (16) only affects the nonlinear dynamics of the system, once the amplitude gets large enough. We leave the nonlinear analysis of the adapted model to future work.
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 reservoir 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 Sect. 4).
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 outburstflood type behaviour. The stability diagram of Fig. 8 was generated for parameter values chosen to represent single reservoirs dammed by glaciers with a small (1:50) surface slope. For the more limited flow paths of 50 and 10 km length, any small reservoir like a large moulin (say, with a maximum side length of 10 m, or a base area of V_{p} of no more than 10^{−4} km^{2}) is stable for even small discharge rates of 1 m^{3} s^{−1}; only a very long flow path of 250 km combined with a single such reservoir could be unstable (Fig. 8b) with (for such a large moulin) moderate flow rates of 10 m^{3} s^{−1}.
The picture changes significantly when we look at shorter, steeper glaciers. We can reduce the size of the parameter space that we need to sample in order to understand different glacier geometries by scaling the problem (Sect. S2.2 of the Supplement). If we change the background hydraulic gradient Ψ_{0} (representing steepness) by a factor χ_{Ψ} and the cavity opening rate v_{o} (representing sliding velocity and bed roughness) by a factor χ_{v}, the effect on stability is equivalent to keeping Ψ_{0} and v_{o} fixed but instead changing q_{in} by a factor of ${\mathit{\chi}}_{\mathrm{v}}^{\mathrm{1}}{\mathit{\chi}}_{\mathrm{\Psi}}$, while simultaneously changing V_{p} by a factor of ${\mathit{\chi}}_{\mathrm{v}}^{[\mathit{\alpha}(n+\mathrm{1}\left)\right]/\left(\mathit{\alpha}n\right)}{\mathit{\chi}}_{\mathrm{\Psi}}^{[\mathrm{2}\mathit{\alpha}n+\mathrm{3}(n+\mathrm{1}\left)\right]/\left(\mathrm{2}\mathit{\alpha}n\right)}={\mathit{\chi}}_{\mathrm{v}}^{.73}{\mathit{\chi}}_{\mathrm{\Psi}}^{\mathrm{2.6}}$ and L by a factor of ${\mathit{\chi}}_{\mathrm{v}}^{(\mathit{\alpha}\mathrm{1})/\left(\mathit{\alpha}n\right)}{\mathit{\chi}}_{\mathrm{\Psi}}^{\mathrm{1}\mathrm{3}/\left(\mathrm{2}\mathit{\alpha}n\right)}={\mathit{\chi}}_{\mathrm{v}}^{.066}{\mathit{\chi}}_{\mathrm{\Psi}}^{.6}$.
The most significant impact here is that of steepening Ψ on the volume V_{p} at which stabilization or destabilization occurs. Consider the concrete example of a 10 km long glacier (Fig. 8c) and steepening its surface slope by a factor 8.5 from 1:50 to 17:100, while also reducing the cavity opening rate v_{o} by a factor of three due to slower sliding. In terms of stability, this is equivalent to keeping the 1:50 surface slope and original cavity opening rate but instead increasing q_{in} by a factor of 25.5, increasing storage capacity V_{p} by a factor 584 and increasing the glacier length L by a factor 3.88 (that is, to 38.8 km). In other words, we would go from Fig. 8c to something that resembles Fig. 8a (in which L=50 km rather than 38.8 km), but with a flux of 25.5 m^{3} s^{−1} as indicated by axis labels in the figure representing an actual flux of 1 m^{3} s^{−1} and a volume V_{p} of $\mathrm{5.84}\times {\mathrm{10}}^{\mathrm{2}}$ km^{2} actually representing a value of 10^{−3} km^{2}. At this point, a relatively modest 10 mby10 m reservoir may become unstable at reasonable water fluxes q_{in}, being located in parameter space towards the righthand end of the stability boundary marked by solid circles in Fig. 8a. As a result, some of the more exotic manifestations of instability explored in Sect. 2.1 are more likely to be within reach of typical water input rates. Figure 12 shows the actual stability diagram for this steeper glacier using the same plotting scheme as Fig. 8.
As our discussion above indicates, relatively small, isolated storage volumes could lead to instabilities either for steeper glaciers or for very long flow paths. The latter in particular are however an unlikely scenario: take the example of 10 mby10 m moulin at the head of a 250 km flow path. In reality, there would most likely be many moulins spread out along the flow path instead. With this as motivation, we discuss a second exotic instability next, in which we opt for the opposite endmember of possible storage geometries by spreading out storage capacity evenly along the flow path. Much of the necessary groundwork is already in place from the analysis in Sect. 3.1.
5.2 Distributed storage: a modification of Nye's instability
In Sect. 2, we formulated a model for a single reservoir, for which Nye's instability predicts the onset of oscillations when the conduit becomes channellike. In Schoof et al. (2014, Sect. 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 cavitylike conduits (which are generally stable for a single reservoir) for the case of such distributed water storage. The resulting oscillations (see Fig. 8 of Schoof et al., 2014; note that the figure captions for Figs. 8 and 9 are switched in the published paper) are likely to be analogous to those found in Dow et al. (2016) for the case of no localized subglacial lake. As in the case of a single reservoir, Schoof et al. (2014) 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 cavitylike.
Here we sketch how to build on Sect. 3.1 to shed further light on the results in Schoof et al. (2014), which were built on the model
Here 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 singlereservoir model (1). Below, we will denote by
the equivalent of V_{p} in Eq. (1). This model also requires boundary conditions; an obvious choice is a prescribed flux q=q_{in} at the inflow x=0, with no actual reservoir there, and again N=0 at the terminus x=L.
With a finite domain size and the proposed boundary conditions above, a steadystate solution to the model (17) will in general have spatial structure as in Schoof et al. (2014), and a linearization around that steady state, analogous to Sect. 3.1, will lead to a boundary value problem that is not amenable to a closedform solution. To avoid this, we concentrate here on shorter length scales and assume that we can use the model (17) with periodic boundary conditions at these scales. This yields a spatially uniform steadystate solution defined implicitly by Eq. (7a) combined with ${q}_{\mathrm{in}}=q(\overline{S},\overline{\mathrm{\Psi}})$ and $\overline{\mathrm{\Psi}}={\mathrm{\Psi}}_{\mathrm{0}}$, where q_{in} is the prescribed flux through the system. This is closely analogous to Eq. (7) but simpler, as the hydraulic gradient here does not contain the gradient term retained in Eq. (7c).
Linearizing as $N=\overline{N}+{N}^{\prime}\mathrm{exp}(ikx+\mathit{\lambda}t)$, $S=\overline{S}+{S}^{\prime}\mathrm{exp}(ikx+\mathit{\lambda}t)$, we find an analogue to the problem (8) as
using the same notation as in Sect. 3.1. The eigenvalue λ is again given by an equation of the form (9), but now with a_{1} and a_{2} given by
The only term that differs intrinsically from Eq. (10) is a_{2}, which now has an imaginary part. We see that the real part of a_{2} is still invariably 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}.
We can identify two mechanisms for instability. The first is essentially the same as Nye's melt–drainage feedback: the more storage capacity there is, the smaller all the terms containing k are in Eq. (20) for a given wavenumber k, and the more likely the first term $({c}_{\mathrm{1}}{q}_{S}+{v}_{\mathrm{o},S}{v}_{\mathrm{c},S})$ in the definition of a_{1} is to dominate the eigenvalue λ, leading to instability for a channellike system. To increase the size of the potentially stabilizing terms therefore requires larger wavenumbers k, and a larger range of short wavelengths is likely to be unstable.
A second instability mechanism can occur if the conduit is cavitylike with $({c}_{\mathrm{1}}{q}_{S}+{v}_{\mathrm{o},S}{v}_{\mathrm{c},S})<\mathrm{0}$ and therefore a_{1}<0. It is still possible to have an eigenvalue with a positive real part and hence instability, because a_{2} has a nonzero imaginary part. One particular case in which an instability due to this term may occur is with limited storage capacity (so ${\overline{v}}_{p}^{\mathrm{1}}$ is large) and at an intermediate wavelengths range (so k is small but not too small). Then it is possible, purely by controlling the size of the storage capacity ${\overline{v}}_{p}$, to ensure that ${a}_{\mathrm{1}}^{\mathrm{2}}\mathrm{4}{a}_{\mathrm{2}}\sim i\mathrm{4}{\overline{v}}_{p}^{\mathrm{1}}{q}_{S}{v}_{\mathrm{c},N}k$ and
Full details can be found in Sect. 6 of the Supplement. Choosing the plus (+) sign ensures an eigenvalue with a positive real part. Importantly, this instability corresponds to a growing wave 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 cavitylike. 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.
To make these considerations more concrete, Fig. 13 shows the stability diagram for the system (17) for the parameter values used in Fig. 8, but with the upstream reservoir V_{p} replaced by distributed storage v_{p} as indicated on the vertical axis. Compared with Fig. 8, it should be clear that there are larger regions of parameter space that exhibit instability for cavitylike conduits (to the left of the vertical portion of the red stability boundary). As in the case of a single reservoir, a larger domain length is destabilizing, but a larger total storage volume is in general necessary to cause instability for channelized systems compared to a localized reservoir.
For instance, for longest flow path of L=250 km, one 10 mby10 m moulin per kilometre of flow path (equivalent to ${v}_{p}={\mathrm{10}}^{\mathrm{4}}$ km) is marginally unstable at q_{in}=10 m^{3} s^{−1} in Fig. 13b, while a single such moulin at the upstream end of the conduit is unstable at even larger q_{in}. Put another way, compared with a single reservoir, distributed storage favours instability at smaller q_{in}. The possibility of oscillatory flow occurring at low q_{in} outside of the melt season due to this instability was previously explored in Schoof et al. (2014).
As demonstrated previously using a much more restricted sweep of parameter space in Schoof et al. (2014), we have shown that drainage systems capable of switching spontaneously between channel and cavitylike 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 lowtomoderate water input, where a simplified, lumped model of reservoir drainage generally reproduces the results of a more sophisticated, spatially extended 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 the flood initiation mechanism can prevent flood magnitude from increasing progressively from cycle to cycle as observed in Ng (1998) and Fowler (1999), as well as 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. The alternative mechanism here is the ability of a cavitylike drainage system that remains open during the reservoir recharge phase to switch to a channelized drainage mode.
The fact that we illustrate this by showing evolution towards a limit cycle is not at odds with the chaotic behaviour observed by Kingslake (2015) using a variant of Fowler's (1999) model: the chaotic behaviour is intrinsically the result of a timevarying 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 timevarying water input rates q_{in}(t). The point of demonstrating 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.
The lower cutoff to the drainage instability that leads to outburst floods corresponds to the drainage system switching to a cavitylike state under steady flow conditions when water input to the reservoir can be drained steadily by those cavities. The mechanism for the cutoff 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 the lumped model relatively to its spatially extended counterpart, which develops more wavelike instabilities at higher water throughput rates.
In closing, we have also investigated how wavelike instabilities can occur when the water reservoir is not localized but 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 channellike 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 that causes water to accumulate in regions where effective pressure is already low arises.
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.
This appendix provides only a brief sketch of the derivation of an asymptotic solution for a limit cycle solution in the case where the reservoir volume is large and inflow is sufficient to ensure that the reservoir cannot be drained by a cavitylike conduit but also that the conduit size evolves much faster than the timescale over which the reservoir fills. Full details are given in Sect. S4 of the Supplement.
The solution is developed from a parametric limit of the model (2) with Eq. (1g), where we assume S_{0}=∞. We scale as ${N}^{**}=N/[N{]}^{\prime}$, ${S}^{**}=S/[S{]}^{\prime}$, ${t}^{**}=t/\left[t\right]$, where the scales are defined through $[S{]}^{\prime}/[t{]}^{\prime}={c}_{\mathrm{1}}{c}_{\mathrm{3}}[S{{]}^{\prime}}^{\mathit{\alpha}}{\mathrm{\Psi}}_{\mathrm{0}}^{\mathrm{3}/\mathrm{2}}={c}_{\mathrm{2}}[S{]}^{\prime}[N{{]}^{\prime}}^{n}$ and ${V}_{p}[N{]}^{\prime}/[t{]}^{\prime}={c}_{\mathrm{1}}[S{{]}^{\prime}}^{\mathit{\alpha}}{\mathrm{\Psi}}_{\mathrm{0}}^{\mathrm{1}/\mathrm{2}}$. Substituting and omitting the asterisks immediately, we obtain
where $\mathit{\nu}=[N{]}^{\prime}/({\mathrm{\Psi}}_{\mathrm{0}}L)$, $\mathit{\delta}={u}_{\mathrm{b}}{h}_{\mathrm{r}}/\left({c}_{\mathrm{2}}\right[S{]}^{\prime}\left[N{{]}^{\prime}}^{n}\right)$ and $\mathit{\u03f5}={q}_{\mathrm{in}}/\left({c}_{\mathrm{1}}\right[S{{]}^{\prime}}^{\mathit{\alpha}}{\mathrm{\Psi}}_{\mathrm{0}}^{\mathrm{1}/\mathrm{2}})$.
We assume that the exponents n and α satisfy $n+\mathrm{1}>\mathit{\alpha}>\mathrm{1}$, in which case the asymptotic solution we develop is valid when
and ν≲1.
The main flood phase is described by omitting terms of O(δ) and O(ϵ)
where matching with the initiation phase described later requires that $(N,S)\to (\mathrm{0},\mathrm{0})$ as $t\to \mathrm{\infty}$. The transformation $P=N/S$ allows the system (A2) to be recast 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={\stackrel{\mathrm{\u0303}}{N}}_{f}$ as t→∞; this then sets the amplitude of the lake level fluctuations that give the asymptotic formula (14) for the flood cycle period.
The refilling phase is described by the rescaling $\stackrel{\mathrm{\u0303}}{N}=N$, $\stackrel{\mathrm{\u0303}}{S}={\mathit{\delta}}^{\mathrm{1}}S$ and $\stackrel{\mathrm{\u0303}}{t}=\mathit{\u03f5}(t{t}_{f})$, where t_{f} is the time of the last flood. At leading order,
conduit size is quasisteady and cavitylike, while lake level evolves purely because of inflow.
The refilling phase ends as N→0 and cavity size becomes large. The relevant rescaling becomes
and gives at leading order
where we have assumed $\mathit{\u03f5}\sim {\mathit{\delta}}^{(\mathit{\alpha}\mathrm{1})(n+\mathrm{1})/\left(\mathit{\alpha}n\right)}$; the alternative case of ϵ∼δ is described in the Supplement. The key to Eq. (A4) is that it predicts finitetime blowup in $\widehat{S}$ at some finite $\widehat{N}={\widehat{N}}_{\mathrm{c}}$ that depends purely on the value of $\mathit{\u03f5}{\mathit{\delta}}^{(\mathit{\alpha}\mathrm{1})(n+\mathrm{1})/\left(\mathit{\alpha}n\right)}$. This is the smallest value of effective pressure (rescaled, of course) that is reached during the drainage cycle. The larger $\mathit{\u03f5}\sim {\mathit{\delta}}^{(\mathit{\alpha}\mathrm{1})(n+\mathrm{1})/\left(\mathit{\alpha}n\right)}$, the smaller and eventually more negative ${\widehat{N}}_{\mathrm{c}}$ becomes; there is therefore a critical value of $\mathit{\u03f5}\sim {\mathit{\delta}}^{(\mathit{\alpha}\mathrm{1})(n+\mathrm{1})/\left(\mathit{\alpha}n\right)}$ at which ${\widehat{N}}_{\mathrm{c}}=\mathrm{0}$. When rendered in dimensional form, this value gives the formula (15).
The MATLAB code used in the computations reported is included in the Supplement, except for the code used to solve the transient calculations displayed in Figs. 9–11, which is included in the Supplement to Rada and Schoof (2018).
The supplement related to this article is available online at: https://doi.org/10.5194/tc1431752020supplement.
The author declares that there is no conflict of interest.
Discussions with Rob Vogt, Ian Hewitt and Garry Clarke are gratefully acknowledged, as are reviews by Mauro Werder and an anonymous referee. This work was supported by an NSERC Discovery Grant.
This research has been supported by the Natural Sciences and Engineering Research Council of Canada (grant no. RGPIN201804665).
This paper was edited by Daniel Farinotti and reviewed by Mauro Werder and one anonymous referee.
Anderson, S., Walder, J., Anderson, R., Kraal, E., Cunico, M., Fountain, A., and Trabant, D.: Integrated hydrologic and hydrochemical observations of Hidden Creek Lake jökulhlaups, Kennicott Glacier, Alaska, J. Geophys. Res., 108, 6003, https://doi.org/10.1029/2002JF000004, 2003. a
Bigelow, D., Flowers, G., Schoof, C., Mingo, L., Young, E., and Connal, B.: The role of englacial hydrology in the filling and drainage of an ice‐dammed lake, Kaskawulsh Glacier, Yukon, Canada, J. Geophys. Res., 125, e2019JF005110, https://doi.org/10.1029/2019JF005110, 2020. a
Björnsson, H.: Jökulhlaups in Iceland: prediction, characteristics and simulation, Ann. Glaciol., 16, 95–106, 1988. a, b
Clague, J., Huggel, C., Korup, O., and McGuire, B.: Climate change and hazardous processes in high mountains, Revista de la Asociación Geológica Argentina, 69, 328–338, 2012. a
Clarke, G.: Glacier outburst floods from “Hazard Lake”, Yukon Territory, and the problem of flood magnitude prediction, J. Glaciol., 28, 3–21, 1982. a
Clarke, G.: Hydraulics of subglacial outburst floods: new insights from the SpringHutter formulation, J. Glaciol., 49, 299–313, 2003. a, b, c
Dow, C. F., Werder, M. A., Nowicki, S., and Walker, R. T.: Modeling Antarctic subglacial lake filling and drainage cycles, The Cryosphere, 10, 1381–1393, https://doi.org/10.5194/tc1013812016, 2016. a
Fisher, D.: Subglacial leakage of Summit Lake, British Columbia, by dye determinations, IASH Pub., 19, 111–116, 1973. a
Flowers, G., Björnsson, H., Palsson, F., and Clarke, G.: A coupled sheetconduit mechanism for jökulhlaup propagation, Geophys. Res. Lett., 31, L05 401, https://doi.org/10.1029/2003GL019088, 2004. a, b
Fowler, A.: Breaking the seal at Grimsvötn, J. Glaciol., 45, 506–516, 1999. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t
Fowler, A. and Ng, F.: The rôle of sediment transport in the mechanics of jökulhlaups, Ann. Glaciol., 22, 255–259, 1996. a
Gudmundsson, M., Sigmundsson, F., and Björnsson: Icevolcano interaction of the 1996 Gjálp subglacial eruption, Vatnajökull, Iceland, Nature, 389, 954–957, 1997. a
Hewitt, I.: Seasonal changes in ice sheet motion due to melt water lubrication, Earth Planet. Sc. Lett., 371, 16–25, 2013. a
Hewitt, I. and Fowler, A.: Seasonal waves on glaciers, Hydrol. Process., 22, 3919–3930, 2008. a
Hewitt, I., Schoof, C., and Werder, M.: Flotation and free surface flow in a model for subglacial drainage. Part 2. Channel flow, J. Fluid Mech., 702, 157–188, 2012. a, b, c, d, e, f, g, h
Huss, M., Bauder, A., Werder, M., Funk, M., and Hock, R.: Glacierdammed lake outburst events of Gornersee, Switzerland, J. Glaciol., 53, 189–200, 2007. a
Kessler, M. and Anderson, R.: Testing a numerical glacial hydrological model using spring speedup events and outburst floods, Geophys. Res. Lett., 31, L18503, https://doi.org/10.1029/2004GL020622, 2004. a, b, c
Kingslake, J.: Modelling icedammed lake drainage, PhD thesis, The University of Sheffield, 2013. a, b
Kingslake, J.: Chaotic dynamics of a glaciohydraulic model, J. Glaciol., 61, 493–502, 2015. a, b
Kingslake, J. and Ng, F.: Modelling the coupling of flood discharge with glacier flow during jokulhlaups, Ann. Glaciol., 54, 25–31, 2013. a
Ng, F.: Mathematical Modelling of Subglacial Drainage and Erosion, PhD thesis, Oxford University, available at: https://ora.ox.ac.uk/objects/uuid:346ee63ca1764381aafb980057dd6f5c (last access: 20 August 2020), 1998. a, b, c, d, e, f, g, h, i, j, k, l, m
Ng, F.: Canals under sedimentbased ice sheets, Ann. Glaciol., 30, 146–152, 2000. a
Nye, J.: Water flow in glaciers: jökulhlaups, tunnels and veins, J. Glaciol., 17, 181–207, 1976. a, b, c, d, e, f, g, h
Post, A. and Mayo, L.: Glacier dammed lakes and outburst floods in Alaska, Hydrologic Atlas 455, USGS, 1971. a
Rada, C. and Schoof, C.: Channelized, distributed, and disconnected: subglacial drainage under a valley glacier in the Yukon, The Cryosphere, 12, 2609–2636, https://doi.org/10.5194/tc1226092018, 2018. a
Röthlisberger, H.: Water pressure in intra and subglacial channels, J. Glaciol., 11, 177–203, 1972. a, b, c
Schoof, C.: Icesheet acceleration driven by melt supply variability, Nature, 468, 803–806, 2010. a, b, c, d, e, f, g
Schoof, C., Hewitt, I., and Werder, M.: Flotation and free surface flow in a model for subglacial drainage. Part 1. Distributed drainage, J. Fluid Mech., 702, 126–156, 2012. a, b, c, d
Schoof, C., Rada, C. A., Wilson, N. J., Flowers, G. E., and Haseloff, M.: Oscillatory subglacial drainage in the absence of surface melt, The Cryosphere, 8, 959–976, https://doi.org/10.5194/tc89592014, 2014. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q
Schuler, T. and Fischer, U.: Modeling the dirunal variation of tracer transit velocity through a subglacial channel, J. Geophys. Res., 114, F04 017, https://doi.org/10.1029/2008JF001,238, 2009. a
Strogatz, S.: Nonlinear Dynamics and Chaos, Perseus Books, Reading, Mass., ISBN 0201543445, 1994. a
Tsai, V. and Rice, J.: A Model for Turbulent Hydraulic Fracture and Application to Crack Propagation at Glacier Beds, J. Geophys. Res., 115, F03007, https://doi.org/10.1029/2009JF001,474, 2010. a
Werder, M.: Modeling Antarctic subglacial lake filling and drainage cycles, Geophys. Res. Lett., 10, 1381–1393, 2016. a
Werder, M., Hewitt, I., Schoof, C., and Flowers, G.: Modeling channelized and distributed subglacial drainage in two dimensions, J. Geophys. Res., 118, 2140–2158, https://doi.org/10.1002/jgrf.20146, 2013. a, b
Wiggins, S.: Introduction to Applied Nonlinear Dynamical Systems, Text in Applied Mathematics, Springer, New York, 2003. a, b, c