Modeling oscillations in connected glacial lakes

Abstract

Mountain glaciers and ice sheets often host marginal and subglacial lakes that are hydraulically connected through subglacial drainage systems. These lakes exhibit complex dynamics that have been the subject of models for decades. Here we introduce and analyze a model for the evolution of glacial lakes connected by subglacial channels. Subglacial channel equations are supplied with effective pressure boundary conditions that are determined by a simple lake model. While the model can describe an arbitrary number of lakes, we solve it numerically with a finite element method for the case of two connected lakes. We examine the effect of relative lake size and spacing on the oscillations. Complex oscillations in the downstream lake are driven by discharge out of the upstream lake. These include multi-peaked and anti-phase filling–draining events. Similar filling–draining cycles have been observed on the Kennicott Glacier in Alaska and at the confluence of the Whillans and Mercer ice streams in West Antarctica. We further construct a simplified ordinary differential equation model that displays the same qualitative behavior as the full, spatially-dependent model. We analyze this model using dynamical systems theory to explain the appearance of filling–draining cycles as the meltwater supply varies.

1. Introduction

Ice-dammed and subglacial lakes have been observed to periodically fill and drain. The classical example of this behavior is the cyclic jökulhlaups from the Grímsvötn caldera lake on the Vatnajökull Ice Cap, Iceland (Nye, 1976; Fowler, 1999; Björnsson, 2003; Fowler, 2009). The Grímsvötn jökulhlaup cycle is characterized by a slow filling phase that lasts several years followed by a sudden draining of water beneath the glacier that results in intense flooding of proglacial rivers and outwash plains (Björnsson, 2003). Predicting the timing and magnitude of jökulhlaups is important for hazard mitigation in many areas (Hewitt and Liu, 2010; Carrivick, 2011; Kingslake and Ng, 2013b).

Glacial lakes are often hydraulically connected through subglacial drainage systems. Flooding from an upstream lake thereby strongly influences the behavior of the downstream lake. Synchronous filling–draining cycles have been observed on the Kennicott Glacier in the Wrangell Mountains of Alaska (Anderson and others, 2005; Bartholomaus and others, 2008, 2011; Bueler, 2014; Brinkerhoff and others, 2016). Hidden Creek Lake sits upstream of the much smaller Donoho Falls Lake. Hidden Creek drains on a yearly cycle, causing Donoho Falls to fill and drain in response to subglacial water pressure variations. To first order, Donoho Falls acts as a pressure gauge that records the passing of the flood (Anderson and others, 2005).

In Antarctica, hundreds of subglacial lakes have been discovered with satellite altimetry and ice-penetrating radar (Smith and others, 2009; Wright and Siegert, 2012). Observations of surface elevation changes within subglacial drainage catchments show that subglacial lakes are hydraulically connected (Wingham and others, 2006; Fricker and Scambos, 2009). These active subglacial lakes often exist in groups beneath Antarctic ice streams. At the confluence of the Whillans and Mercer ice streams in West Antarctica, complex oscillations have been observed where downstream lakes display multiple drainage events per filling–draining cycle of the upstream lake (Fricker and Scambos, 2009; Siegfried and others, 2016; Siegfried and Fricker, 2018). Draining events have been linked to periods of enhanced ice flow in the Mercer–Whillans system (Siegfried and others, 2016), Byrd Glacier that drains East Antarctica through the Transantarctic Mountains (Stearns and others, 2008) and Crane Glacier on the Antarctic Peninsula (Scambos and others, 2011). However, similar lake drainage events on Thwaites Glacier suggested no strong connection between drainage and ice flow (Smith and others, 2017). The coupling between lake drainage and ice flow is uncertain partly because these complex filling–draining cycles are not well-understood.

Models for ice-dammed lake drainage were originally proposed to help explain the Grímsvötn jökulhlaup cycle (Nye, 1976; Fowler, 1999). Fowler's (1999) model of single-lake oscillations supplies a modified version of the subglacial channel equations from Nye (1976) with boundary conditions that are determined by a simple lake model. Subsequent modeling studies investigated jökulhlaup discharge magnitudes (Ng and Björnsson, 2003), the effect of seasonal meltwater input (Kingslake, 2015), the coupling between flooding and basal sliding (Kingslake and Ng, 2013a), jökulhlaup predictability (Fowler, 2009; Kingslake and Ng, 2013b), ice cauldron formation (Evatt and Fowler, 2007) and applications to Antarctic subglacial lakes (Evatt and others, 2006; Carter and others, 2017). While single-lake systems have been an active research area, there have been few theoretical studies of connected lake systems (Peters and others, 2009; Carter and others, 2017).

Here, we introduce a model for connected glacial lakes that extends previous single-lake models (Fowler, 1999, 2009). The aims of this paper are (I) to explore the effect of relative lake size and spacing on the nature of filling–draining oscillations and (II) to study how meltwater supply variations lead to the appearance and disappearance of these oscillations. We use numerical experiments to address these aims and discuss physical mechanisms behind the coupled oscillations. We then derive a reduced ordinary differential equation (ODE) model that has similar qualitative behavior to the full partial differential equation (PDE) model but is orders of magnitude faster and more amenable to analysis. We use this reduced model to analyze the transition from periodic flooding to steady water discharge. We conclude by comparing numerical results with data from natural lake systems to assess the model and its implications.

2. The model

We assume that a sequence of ice-dammed lakes are hydraulically connected by subglacial channels (Fig. 1). The mathematical model consists of two components. The first component is a set of PDEs that governs the evolution of the subglacial channels. The second component is a simple evolution equation for the lakes. Following Fowler (1999, 2009), the lake evolution equations provide the boundary conditions for the subglacial channels where they are connected to the lakes. While we only explore two connected lakes here, we formulate the model for an arbitrary number of lakes.

Fig. 1. Schematic diagram of connected lake model setup. Upstream and downstream lake depths are noted by h ₁ and h ₂, respectively. The upstream and downstream channel lengths are noted by L ₁ and L ₂. The meltwater supply rate to the lakes is Q _in. Arrows in the subglacial channels note the water flow direction.

2.1 Subglacial channel model

We let S be the cross-sectional area of the semi-circular subglacial channel, Q the water discharge, Ψ the hydraulic gradient, x the distance along the channel and t time. We define the effective pressure N to be the difference between the ice overburden pressure and the water pressure in the channel. We assume that the water in the channel and the ice at the channel walls are always at the melting point and we neglect heat advection (Nye, 1976). We further assume that the ice creeps as a shear-thinning viscous fluid (Nye, 1953; Glen, 1955; Cuffey and Paterson, 2010).

The magnitude of water discharge through the subglacial channel is partially determined by the hydraulic gradient

(1)$$\Psi = \psi(x) + \displaystyle{{\partial N}\over{\partial x}},$$

which includes the effective pressure gradient and a background hydraulic gradient, ψ. We choose the parameterization

(2)$$\psi (x) = \psi_0[1-2 \exp(-20x/L_{\rm c})],$$

where L _c is the length of the subglacial channel and ψ₀ is the hydraulic gradient scale. This parameterization is based on a fit to surface and bed slopes at Grímsvötn (Fowler, 2009). We use this parameterization because it allows water flow toward the lake between floods, which has previously been found to promote long-term stability of filling–draining oscillations (Fowler, 2009; Kingslake, 2015). Within ~ L _c/20 km of the lake, the hydraulic gradient Ψ is strongly influenced by ψ. We expect that other parameterizations would also allow for stable filling–draining cycles based on the success of the spatially-uniform model in ‘Simplified model’ section.

Conservation of ice mass is satisfied through a balance of melting and viscous closure rates at the channel walls. The channel size evolves according to

(3)$$\displaystyle\displaystyle{{\partial S}\over{ \partial t}} = \displaystyle{{1}\over{\rho_{\rm i} l}} Q\Psi - 2\hat{A}n^{-n} S N \vert N \vert^{n-1},$$

where ρ_i is the ice density, l is the latent heat of melting, and $\hat {A}$ and n are the parameters in Glen's flow law (Glen, 1955; Cuffey and Paterson, 2010). In Eqn (3), the melt and closure rates are proportional to QΨ and SN|N|ⁿ⁻¹, respectively.

Conservation of mass for water in the channel requires

(4)$$\displaystyle{{\partial S}\over{\partial t}} + \displaystyle{{\partial Q}\over{\partial x}} = \displaystyle{{1}\over{\rho_{\rm w} l}} Q\Psi + \mu,$$

where ρ_w is the water density and μ is a constant source term for water that enters the channel through drainage pathways such as cavities and moulins (Fowler, 1999, 2009; Kingslake, 2015). We note that μ has a stabilizing effect because it prevents the channel from closing completely. Combining Eqn (4) with Eqn (3), we obtain

(5)$$\displaystyle{{\partial Q}\over{\partial x}} = 2\hat{A}n^{-n} S N \vert N \vert^{n-1} - \displaystyle{{\rho_{\rm w}-\rho_{\rm i}}\over{\rho_{\rm w}} \rho_{\rm i} l} Q \Psi + \mu.$$

Finally, the Darcy–Weisbach law for turbulent water flow through pipes provides a relation between channel size, pressure and discharge through the equation

(6)$$Q = \left(\displaystyle{{2}\over{\pi}}\right)^{{1}/{4}}\sqrt{\displaystyle{{2+\pi}\over{\rho_{\rm w}} f_{\rm d}}} S^{\alpha}\Psi \vert\Psi\vert^{-{1}/{2}},$$

where f _d is the friction factor (Clarke, 2003). We list the numerical values of the above parameters in Table 1.

Table 1. Model parameters and physical constants used in simulations. These values are fixed unless noted otherwise.

A − in the fourth column denotes a dimensionless quantity.

2.2 Lake model

We suppose that the first upstream lake (Lake 1 in Fig. 1) is connected only to a downstream channel. We let h be the lake depth, H the ice thickness and A the lake surface area. We hold the ice thickness and lake surface area constant. Alternatively, if the geometry of the lake is known then one may prescribe the surface area as a function of the lake depth. We let Q _in be the volumetric rate of meltwater flow into the lake from runoff, englacial storage and other sources. Below, we denote the water discharged from the lake through the downstream channel by $Q_{\rm out}^{\rm ch}$. The lake depth evolves according to

(7)$$\displaystyle{{{\rm d}h}\over{{\rm d}t}} = \displaystyle{{1}\over{A}}[Q_{\rm in} - Q_{\rm out}^{\rm ch}].$$

Assuming hydrostatic balance, the effective pressure at the lake is

(8)$$ N_{\rm L} = \rho_{\rm i} g H - \rho_{\rm w} g h.$$

This relation between lake depth and effective pressure in Eqn (8) leads to the evolution equation

(9)$$\displaystyle{{{\rm d} N_{\rm L}}\over{{\rm d} t}} = \displaystyle{{\rho_{\rm w} g}\over{A}}[Q_{\rm out}^{\rm ch}-Q_{\rm in}].$$

For downstream lakes (e.g., Lake 2 in Fig. 1), there is an additional input $Q_{\rm in}^{\rm ch}$ from the upstream channel. In this case, the effective pressure at the lake evolves as

(10)$$\displaystyle{{{\rm d} N_{\rm L}}\over{{\rm d} t}} = \displaystyle{{\rho_{\rm w} g}\over{A}}[Q_{\rm out}^{\rm ch}-Q_{\rm in}-Q_{\rm in}^{\rm ch}].$$

2.3 Boundary conditions

The lake evolution Eqn (10) determines the boundary conditions for the channel equations everywhere except at the terminus. At the terminus, we impose the boundary condition N = 0. As discussed by Evatt (2015), the zero pressure condition at the terminus leads to the problem of a forever-widening channel there. Although observations suggest that subglacial channels widen at the terminus (Drews and others, 2017), unbounded growth is unreasonable. Therefore, we scale the viscous closure law by

(11)$$\displaystyle \zeta(S) = \textstyle \bigg[1 -\bigg(\displaystyle{{S}\over{S_{\rm f}}}\bigg)^{{{1}\over{n}}}\bigg]^{-n},$$

which arises from the assumption of finite ice depth (Evatt, 2015). The viscous closure rate in Eqn (3) becomes SN|N|ⁿ⁻¹ζ(S). For large S _f, ζ(S) ≈ 1 away from the terminus. Most importantly, S can never exceed S _f (Evatt, 2015). We have found that using S _f = 8000 m² in numerical experiments does not alter the solutions near the lake qualitatively and prevents the channel from growing without bound.

2.4 Scaling

Here we introduce a scaling for the subglacial channel Eqns (1–6). We introduce the dimensionless variables

(12)$$\eqalign{S_{\ast} \quad & = \quad S/\tilde{S}, \quad N_{\ast} = N/\tilde{N}, \quad Q_{\ast} = Q/\tilde{Q}, \cr \psi_{\ast}\quad & = \quad \psi/\tilde{\psi}, \quad x_{\ast}= x/\tilde{x}, \quad t_{\ast} = t/\tilde{t},}$$

and choose the scales

(13)$$\eqalign{\tilde{S} & = \displaystyle{{Q_{\rm in}^{{1}/{\alpha}}}\over{\left(\displaystyle{{2}\over{\pi}}\right)^{{1}/{4\alpha}}\bigg({\psi}_0 \displaystyle{{2+\pi}\over{\rho_{\rm w}f_{\rm d}}}\bigg)^{{1}/{2\alpha}} }},\quad \tilde{N} = {\psi}_{0} L_{\rm c}, \quad \tilde{Q} = Q_{\rm in}, \cr \tilde{\psi} & = \psi_0, \quad \tilde{t} = \displaystyle{{A_1 \tilde{N}}\over{\rho_{\rm w} g Q_{\rm in}}}, \quad \tilde{x} = L_{\rm c}},$$

where A ₁ is the surface area of the first (upstream) lake. The scale $\tilde {S}$ is the channel size at the equilibrium discharge Q = Q _in in the absence of effective pressure gradients. The scale $\tilde {t}$ is the intrinsic filling timescale of the first upstream lake. Likewise, the melting timescale for the channel at the equilibrium discharge is

(14)$$t_{\rm m} = \displaystyle{{\tilde{S} \rho_{\rm i} l}\over{ Q_{\rm in}\psi_0}}.$$

The quantity

(15)$$t_{\rm c} = \displaystyle{{1}\over{2\hat{A}n^{-n} \tilde{N}^n}}$$

is the viscous closure timescale for the channel, which depends strongly on the channel length. We define the dimensionless parameters

(16)$$\displaystyle a_{\rm m} = \tilde{t}/t_{\rm m} \propto Q_{\rm in}^{-{{1}\over{\alpha}}}L_{\rm c}, \quad a_{\rm c} = \tilde{t}/t_{\rm c}\propto Q_{\rm in}^{-1}L_{\rm c}^{n+1},$$

which are the melting and closure timescales normalized by the filling timescale. The parameters a _m and a _c control the magnitudes of the melting and closure rates, respectively.

2.5 Generalized connected lake system

Here we complete the formulation of the model by generalizing Eqns (1)–(6) and the boundary condition Eqn (10). We consider a system with M lakes and channels on a glacier of length L _c. We let L _k be the length of the k-th subglacial channel, subject to

(17)$$\sum\limits_{k=1}^M L_{k} = L_{\rm c}.$$

The scaled length of each subglacial channel is

(18)$$L_k^{\ast} = L_{k}/L_{\rm c},$$

so that the total length of the scaled system is one. Supposing that the system (1)–(6) holds on each interval (0,  L _k*), the scaling described by Eqns (12–16) leads to the system of partial differential equations

(19)$$\displaystyle{{\partial S_k} \over{\partial t_{\ast}}} = a_{\rm m} Q_k \Psi_k - a_{\rm c} S_k N_k \vert N_k \vert^{n-1} \zeta (S_k),$$

(20)$$\displaystyle{{\partial Q_k}\over{\partial x_{\ast}}} = b_{\rm c} S_k N_k \vert N_k \vert^{n-1} \zeta(S_k) - b_{\rm m} Q_k \Psi_k + \mu_{\ast},$$

(21)$$Q_k = S_k^{\alpha}\Psi_k \vert\Psi_k\vert^{-\textstyle{{1}\over{2}}},$$

(22)$$\Psi_k = \psi_{\ast} + \displaystyle{{\partial N_k}\over{\partial x_{\ast}}},$$

for $k=1\ldots M$. Above, $\mu_{\ast} = \mu \tilde {x}/\tilde {Q}$, $b_{\rm m} = (({\rho_{\rm w}-\rho_{\rm i}})/{\rho_{\rm w} \rho _{\rm i} l}) \tilde {N}$, and $b_{\rm c} = 2\hat {A}n^{-n} \tilde {x} \tilde {S} \tilde {N}^n / \tilde {Q}$. The parameters b _m and b _c are small and the associated terms may be neglected, although we do not do so here. We let A _k be the surface area of each lake. To satisfy the lake evolution Eqn (10), we enforce the non-dimensionalized upstream boundary conditions

(23)$$\displaystyle{{\partial N_k}\over{\partial t_{\ast}}}(0,t_{\ast}) = \displaystyle{{1}\over{A_k^{\ast}}} \left[Q_k(0,t_{\ast}) - 1 - Q_{k-1}(L_{k-1}^{\ast},t_{\ast}) \right]$$

for k > 1, where we have defined A _k* = A _k/A ₁. We refer to A _k* as the relative size of the k-th lake. For the first lake, the non-dimensional boundary condition is

(24)$$\displaystyle{{\partial N_1}\over{\partial t_{\ast}}}(0,t_{\ast}) = Q_1(0,t_{\ast}) - 1.$$

Likewise, the downstream boundary conditions are

(25)$$\displaystyle{{\partial N_k}\over{\partial t_{\ast}}}(L_k^{\ast},t_{\ast}) = \displaystyle{{1}\over{A_{k+1}^{\ast}}}[ Q_{k+1}(0,t_{\ast}) - 1 - Q_k(L_k^{\ast},t_{\ast})]$$

for k < M, with the modification at the terminus

(26)$$N_M(L_M^{\ast},t_{\ast}) = 0.$$

Above we have assumed that the effective pressure is uniform across each lake. Moreover, we assume that $L_k^{\ast}\; \gtrsim\; 1/20$ so that the parameterization ψ_*(x _*) does not dominate the hydraulic gradient Ψ_k. Equations (19)–(26) constitute a complete mathematical model for connected glacial lakes. When discussing single-lake systems, we will omit the k subscripts on the variables. The system is solved numerically using a finite element method. We provide details on the numerical method and implementation in Appendix A.

3. Results

3.1 Single-lake model

We first study the behavior of a single-lake system (M=1) to develop physical intuition into the oscillatory nature of the model before considering the multiple-lake model. We assume that Q _in is constant in time for all experiments. In the figures and discussion below we denote the discharge out of the lake by Q _L = Q(0,  t). For clarity, we plot the lake depth h(t) rather than the effective pressure at the lake N(0,  t). These quantities are related through Eqn (8). The single-lake results outlined below may be compared to previous modeling studies (Fowler, 1999, 2009; Kingslake and Ng, 2013a; Kingslake, 2015).

The fundamental feature of the model is that the boundary condition in Eqn (24) forces oscillations in lake depth, effective pressure, discharge and channel size (Figs 2a, b). The physical mechanism for single-lake oscillation originates from the coupling of the melting and closure rates to discharge and channel size. The lake depth h increases slowly until the channel melts open enough for rapid drainage to begin (point h _max in Fig. 2b). The discharge Q _L grows during the flood, causing the channel to open rapidly due to melting. The positive feedback between discharge and channel size causes a sharp increase in melt rate. At the same time, increasing discharge causes the lake level to decrease rapidly. This melt-opening feedback persists to the peak of the flood (point Q _max in Fig. 2b). Eventually the channel begins to collapse once it becomes too large. Water flow through the channel becomes increasingly constricted as the channel closes, diminishing the melt rate and ending the flood. The lake then begins to refill because the channel is small again (Fig. 2b).

Fig. 2. Numerical simulations of a single-lake system for Q _in = 5 (a and b) and 15 m³s⁻¹ (c and d). The solutions are plotted over time (a and c) and in Q _L − h phase space (b and d). Arrows in panels (b) and (d) show the direction of the trajectories in phase space. Maximum lake depth and discharge are noted by the points labeled h _max and Q _max, respectively. The lake surface area is A = 1 km². Note that when Q _in = 15 m³s⁻¹ (panels c and d), the discharge approaches the equilibrium value Q _L = Q _in. The dimensionless parameters are (a _m,  a _c) = (66,  2900) (a and b) and (a _m,  a _c) = (28,  967) (c and d).

The oscillations in h and Q _L approach periodic functions over time for a range of meltwater input (Q _in) values (Fig. 2a). These periodic functions define a limit cycle – a closed curve in Q _L-h phase space (Fig. 2b). All solutions that start away from equilibrium will approach the limit cycle except at higher values of Q _in. The behavior at high Q _in is explored in the next paragraph and in greater detail in ‘Simplified model’ section. The oscillation frequency is controlled primarily by the meltwater input Q _in and the lake size A because the natural filling timescale for the lake $\tilde {t}$ is proportional to A/Q _in in Eqn (13). Large lakes oscillate slowly because they take longer to fill at a fixed Q _in. Conversely, a higher meltwater supply results in more frequent cycles.

Once Q _in increases past some critical threshold, the oscillations damp toward equilibrium where the lake drains water continuously at a constant rate (Fig. 2c). In phase space, the damping appears as a spiral (Fig. 2d). Damping occurs because the filling timescale $\tilde {t}$ in Eqn (13) becomes too small relative to the melting timescale t _m in Eqn (14) and viscous closure timescale t _c in Eqn (15) for oscillations to be sustainable. Likewise, the melting strength a _m and closure strength a _c parameters in Eqn (16) become small and perturbations to the equilibrium solution are not able to grow. In this scenario, the discharge out of the lake approaches the equilibrium value Q _L = Q _in (Figs 2c, d). We explore the transition from oscillation to damping further in ‘Simplified model’ section with a reduced ODE model.

3.2 Double-lake models

Connected lakes can oscillate in a variety of different ways that are distinct from the single-lake solutions discussed in the previous section. However, these double-lake oscillations follow the same physics of the single-lake system. We explore the amplitudes and frequencies of these oscillations as the downstream lake surface area, and the spacing between the lakes is varied since these are easy parameters to estimate in natural systems. We choose the values of A ₁ = 1 km² and Q _in = 5 m³s⁻¹ in order to determine appropriate values for a _m and a _c. However, these results apply to systems with other lake sizes and meltwater supply rates if the ratio A/Q _in is similar to what we have chosen. Unless stated otherwise, we use (a _m,  a _c) = (66, 2900). This parameter choice corresponds to the single-lake solution with undamped oscillations (Figs 2a, b). The scaling shows that the trends outlined below are due to the relative downstream lake size A ₂* = A ₂/A ₁ and the magnitude of the meltwater supply rate relative to the lake sizes. The behavior of the upstream lake always resembles the single-lake solutions due to the constant meltwater input. The downstream lake behavior is modified by the additional time-dependent water input from the upstream channel, Q ₁(L ₁*,  t), and the parameter variations discussed below.

Lake spacing examples

First, we explore the influence of lake spacing on the depth-discharge oscillations. The length of the whole domain is always preserved since we set L ₂* = 1 − L ₁* and vary L ₁*. In this series of examples, we keep the size of the second lake constant and equal to the first lake (A ₂* = 1). In general, we find that the lakes oscillate in phase when they are closely spaced ($L_1^{\ast} \lesssim 0.5$). When the lakes are far apart ($L_1^{\ast}\; \gtrsim\; 0.5$), we observe small-amplitude oscillations in the downstream lake that are delayed relative to the upstream cycles.

When lakes are closely spaced, both lakes damp toward equilibrium (Fig. 3a). For our parameter choices this occurs when L ₁* < 0.5, but sets in strongly around L ₁* = 0.1. Damping occurs in shorter channels because the closure rate becomes small relative to the melt rate, similar to the transition to damping with increasing meltwater supply rate. In single-lake systems, the closure rate decreases strongly as the channel length decreases because the closure magnitude parameter a _c is proportional to $L_{{\rm c}}^{n+1}$ in Eqn (16). The system approaches equilibrium because the channel cannot close quickly enough to sustain oscillations. The same effect occurs in the double-lake system, although the change is in the length of the subdomains rather than the parameter a _c. Interestingly, the single-lake solution exhibits undamped oscillations for the same choice of parameters (a _m, a _c).

Fig. 3. Simulations of equally-sized lakes (A ₂* = 1) as the scaled upstream channel length L ₁* changes. We choose (a) L ₁* = 0.1, (b) L ₁* = 0.5 and (c) L ₁* = 0.75.

When the lakes are far apart ($L_1^{\ast}\; \gtrsim\; 0.6$), draining of the downstream lake is consistently delayed and the oscillations have smaller amplitudes than the upstream lake (Fig. 3c). Delayed oscillations in distant lakes occur because pressure variations caused by draining of the upstream lake take longer to reach the downstream lake. In each cycle, the downstream lake depth reaches the flotation level h = (ρ _i/ρ _w)H and then drains rapidly. The reduced amplitude of the downstream lake oscillations is a consequence of the downstream lake draining through a shorter channel. The downstream lake continuously drains water between floods because the closure rate is small in the downstream channel, leading to a large inter-flood channel size. The downstream lake oscillations are forced by the upstream floods, in contrast to the self-sustained oscillations in the upstream lake.

If the lakes are evenly spaced (L ₁* ≈ 0.5), they oscillate with similar amplitudes (Fig. 3b). A common feature in this scenario is the presence of multiple peaks per cycle in the downstream depth and discharge timeseries. When the lakes are the same size, as in this example, the downstream lake may have alternating double and single peaks (Fig. 3b). Multiple-peaked cycles are explored further in the following section.

Smaller downstream lake

Now we fix L ₁* = 0.5 and study the effect of decreasing the relative downstream lake size A ₂*. The frequency of downstream lake oscillations increases with decreasing lake size (Fig. 4). In the limit of a very small surface area (A ₂* ≈ 0.01), the downstream lake behaves as a manometer that records the mean effective pressure in the channel (Fig. 4f)(Anderson and others, 2005; Carter and others, 2013). The downstream lake depth is essentially constant except for when the upstream lake begins to drain. The downstream lake then fills rapidly and drains in phase with the upstream lake. Multiple peaks in downstream lake depth per upstream cycle occur when A ₂* is larger than ~0.01 (Figs 4a–e). When multiple peaks exist in the downstream lake depth timeseries, the final peak coincides with draining of the upstream lake. The overall period of the downstream lake's filling–draining cycle conforms to that of the first lake due to the periodic arrival of the additional water input Q ₁(L ₁*,  t). The upstream floods force the downstream lake to oscillate at the same overall frequency.

Fig. 4. Simulations with smaller downstream lake size. The scaled upstream channel length is L ₁* = 0.5. The relative downstream lake sizes are (a) A ₂* = 0.49, (b) 0.36, (c) 0.25, (d) 0.16, (e) 0.09 and (f) 0.01.

If the lakes were not coupled, the downstream lake would resemble the upstream one except with higher frequency. The closest example to uncoupled behavior occurs around A ₂* = 0.49 where the second lake has two large, well-formed peaks (Fig. 4a). The transition from large to small intervening oscillations occurs around A ₂* ≈ 0.36 (Fig. 4b). The small intervening oscillations are not associated with significant discharge out of the downstream lake. The amplitudes of the intervening oscillations decrease continuously as lake size decreases (Figs 4b–f).

Larger downstream lake

When the relative downstream lake size A ₂* is increased above unity, the lakes cease to oscillate in phase (Fig. 5). There is a small size range (A ₂* ≈ 2.25) where the lakes demonstrate antiphase oscillations (Fig. 5a). If the lakes were uncoupled, the downstream lake would oscillate at a lower frequency due to its large size. The additional input Q ₁(L ₁*, t) allows the downstream lake to oscillate at the same frequency if it is not too large. More precisely, the downstream lake will oscillate at a frequency similar to the upstream lake's frequency if $A_2/(Q_{\rm in} + \bar {Q}_1)$ is similar to A ₁/Q _in, where $\bar {Q}_1$ is the mean of Q ₁(L ₁*, t) over one upstream cycle. Because the additional water input from the upstream flood must be received during the downstream filling stage, the lakes must oscillate out of phase. The effect of the additional water input is reflected in the timeseries (Fig. 5a). Between floods, the downstream lake level increases slowly at first, ramps up during the upstream flood, then fills slowly again before draining. The slow-filling stage represents the natural, uncoupled filling rate of the downstream lake.

Fig. 5. Simulations of larger downstream lakes. The scaled upstream channel length is L ₁* = 0.5 for both examples. The relative downstream lake sizes are (a) A ₂* = 2.25 and (b) 16.

If the downstream lake is much larger than the upstream lake, it essentially acts as a storage tank (Fig. 5b). The lake level increases steadily at its natural rate until each upstream flood when a sharp uptick occurs. The sharp upticks superimposed on the natural filling rate lead to a staircase-like appearance in the downstream depth timeseries (Fig. 5b). After several upstream cycles, the downstream lake finally drains, releasing large volumes of water.

Meltwater variations

In contrast to the geometrically simple limit cycles in the single-lake system (Fig. 2b), the downstream lake oscillations (e.g., Fig. 6a) can exhibit complex patterns in phase space (Fig. 6b). As with the single-lake system, damping due to an increase in meltwater supply also occurs in the double-lake system (Figs 6c, d). In this scenario, both lakes approach equilibrium simultaneously. The physical mechanism for damping is the same as in the single-lake system. At high Q _in, the filling timescale $\tilde {t}$ is too fast relative to the melting t _m and closure t _c timescales to support oscillations. The transition from oscillation to stable equilibrium is a fundamental feature of the model. We explore this bifurcation further in the following section.

Fig. 6. Numerical simulations of a double-lake system for Q _in = 5 (a and b) and 15 m³s⁻¹ (c and d). The relative downstream lake size is A ₂* = 1 and the scaled upstream channel length is L ₁* = 0.5. The solutions for the downstream lake are plotted over time (a and c) and in Q _L-h phase space (b and d). The solutions for the upstream lake are not plotted because they are nearly identical to those in Fig. 2. Note that when Q _in = 15 m³s⁻¹ (panels c and d), the discharge approaches the equilibrium value Q _L = 2Q _in. The dimensionless parameters are (a _m,  a _c) = (66,  2900) (a and b) and (a _m,  a _c) = (28,  967) (c and d).

4. Simplified model

Here we explore an ODE model for the connected lake system in order to gain insight about lake behavior as the meltwater supply rate varies. Additionally, solving the ODE model numerically is significantly more computationally efficient than solving the PDE model. We assume that S, N and Q are spatially uniform near each lake. Then, the channel evolution in Eqn (19) and upstream boundary condition (23) form the dynamical system

(27)$$\displaystyle{{{\rm d}{S}_k}\over{{\rm d}t _\ast}} = a_{\rm m} Q_k \Psi_k - a_{\rm c} S_k N_k \vert N_k \vert^{n-1},$$

(28)$$\hskip-34pt \displaystyle{{{\rm d}{N}_k}\over{{\rm d}t _ \ast}} = \displaystyle{{1}\over{A_k^{\ast}}} [ Q_k - Q_{k-1} - 1 ].$$

In order to simplify the analysis, we have neglected the closure rate scale ζ(S) in Eqn (11). Including this leads to better quantitative agreement with the PDE model, but does not change the qualitative structure of the bifurcations. The first lake's upstream boundary condition in Eqn (23) becomes

(29)$${{\rm d}{N}_1\over {\rm d}t _ \ast} = Q_1 - 1.$$

We introduce a regularized Darcy–Weisbach law

(30)$$Q_k = (S_k+\varepsilon_{\ast})^{{\alpha}}\Psi_k \vert\Psi_k\vert^{-{{1}\over{2}}},$$

where $\varepsilon_{\ast} = \varepsilon /\tilde {S}$, for some small parameter ε. Equation (30) ensures that the channel size is bounded below by some positive value. Regularization is necessary to produce the limit cycles and bifurcations that exist in the PDE model. An unregularized version of the simplified model presented here with ε = 0 and constant Ψ was studied previously by Fowler (1999) and Ng and Björnsson (2003). Local stability analysis and numerical experiments were used to suggest that limit cycles do not exist in the unregularized model (Fowler, 1999; Ng and Björnsson, 2003). Using an energy argument, we extend this result by showing that all non-equilibrium solutions to the unregularized model are globally unstable and blow up if α > 1 (Appendix B). However, for α = 1 all non-equilibrium solutions are periodic orbits that resemble the limit cycles in the single-lake system. In the case of α = 1, the level curves of the energy function provide an approximate algebraic relation between S and N during single-lake flood cycles (Appendix B).

The governing equations may be stabilized in several other ways. For example, an equivalent regularization is to modify the closure term to (S − $\varepsilon_*$)N|N|ⁿ⁻¹ so that the channel ceases to close at some small size, while keeping the usual Darcy–Weisbach law. The solutions associated with this regularization will only differ from ours by the small additive constant $\varepsilon_*$. A third possible regularization is to add a small parameter to Eqn (27), which may be interpreted as an opening rate due to basal sliding (Schoof, 2010). These ‘regularizations’ have the same stabilizing effect as the channel source term μ in the PDE model Eqn (4) because they prevent the channel from closing completely.

To complete the simplified model, we introduce a forward difference approximation of the dimensional hydraulic gradient $\Psi _k^{\rm d} = {\psi }_0 - N_k/L_{k}$. We neglect downstream pressure contributions to the upstream hydraulic gradient under the assumption that upstream pressure variations are dominated by upstream lake activity. With this parameterization, the upstream lake is completely decoupled from the downstream lake. This may be inappropriate for closely spaced lakes because the downstream lake may strongly influence the upstream channel. This approximation corresponds to the scaled hydraulic gradient

(31)$$\Psi_k = 1 - \displaystyle{{1}\over{L_k^{\ast}}} N_k.$$

The simplified model described by Eqns (27–31) produces results that are qualitatively consistent with those from the PDE model. All oscillatory behaviors discussed for the double-lake system exist in the simplified model for similar parameter choices (Fig. 7). Moreover, the oscillations damp as Q _in increases past some critical threshold, as in the PDE model (Appendix C). We provide a quantitative comparison of the ODE and PDE models in Appendix D. The advantage of Eqns (27–31) is that they may be studied using dynamical systems theory (Chicone, 2006; Kuznetsov, 2004; Strogatz, 2018).

Fig. 7. Numerical results from the simplified model in Eqns (27 and 28) for double-lake systems compared with the PDE results from 'Results' section. For the ODE results we choose the hydraulic gradient scale ψ₀ = 75 Pa m⁻¹, ε = 0.1 m² and glacier length L _c = 40 km. The scaled upstream channel length is L ₁* = 0.5 except in panels (g) and (h). The parameters we choose are (a) A ₂* = 0.01, (b) A ₂* = 0.01, (c) A ₂* = 0.16, (d) A ₂* = 0.16, (e) A ₂* = 0.81, (f) A ₂* = 1, (g) A ₂* = 0.81 and L ₁* = 0.75, (h) A ₂* = 1.0 and L ₁* = 0.75, (i) A ₂* = 2.25, (j) A ₂* = 2.25, (k) A ₂* = 16 and (l) A ₂* = 16.

4.1 Bautin bifurcation

Single-lake and multiple-lake systems damp toward equilibrium once a critical value of the meltwater source parameter is reached (Figs 2c, d and 6c, d). Based on the numerical experiments, Fowler (1999) suggested that this transition occurs due to a Hopf bifurcation. Numerical results show that the transition occurs because Eqns (27) and (28) undergo a Bautin (generalized Hopf) bifurcation as Q _in varies (Figs 8 and 9). The Bautin bifurcation is a combination of a supercritical Hopf bifurcation, subcritical Hopf bifurcation and saddle-node bifurcation of cycles (Kuznetsov, 2004).

Fig. 8. Phase portraits of the simplified model in Eqns (27 and 29) for a single-lake system with varying meltwater supply. The meltwater supply rates for each panel are (a) Q _in = 5, (b) 8.5, (c) 9 and (d) 15 m³s⁻¹. Limit cycles are plotted for t > 315 years. We choose the parameters A = 1 km², ε = 0.05 m², ψ₀ = 45 Pa m⁻¹ and L _c = 40 km. The dimensionless parameters are (a) (a _m,  a _c) = (10,  108), (b) (a _m,  a _c) = (6.2,  63.4), (c) (a _m,  a _c) = (5.9,  59.8) and (d) (a _m,  a _c) = (4,  36).

Fig. 9. Bifurcation diagram for the ODE model of a single-lake system in Eqns (27 and 29) for various values of ε. The bifurcation parameter is Q _in. Limit cycle amplitudes are defined as max (S(t)) for t > 315 years. Unstable limit cycle amplitudes are computed by solving the time-reversed system of equations. We choose A = 1 km², ψ₀ = 45 Pa m⁻¹ and L _c = 40 km. Labeled points correspond to the (i) supercritical Hopf bifurcation, (ii) subcritical Hopf bifurcation and (iii) saddle-node bifurcation of cycles.

The individual bifurcations occur sequentially as Q _in increases. The supercritical Hopf bifurcation is characterized by the appearance of a stable limit cycle at some small Q _in ≪ 1 (Fig. 9 point i). As Q _in increases, the amplitude of the stable cycle grows (Fig. 8a) until the subcritical bifurcation occurs. At the subcritical Hopf bifurcation, an unstable limit cycle appears inside of the stable cycle (Figs 8b and 9 point ii) and grows rapidly (Fig. 8c). The unstable limit cycle soon annihilates the stable cycle in a saddle-node bifurcation of cycles, leaving behind only a stable spiral node (Figs 8d and 9 point iii). The values of Q _in where these bifurcations occur are noted in Figure 9, but depend on ε and the other parameters that determine a _m and a _c in Eqn (16).

The stability transitions are related to the sign of tr(J), the trace of the Jacobian of Eqns (27) and (28) evaluated at the equilibrium. The trace may be written as

(32)$${\rm tr}(J) = \displaystyle{{\alpha a_{\rm m}}\over{(S_{\rm e}+\varepsilon_{\ast})^{1+2\alpha}}} - \displaystyle{{a_{\rm m}}\over{S_{\rm e}(S_{\rm e}+\varepsilon_{\ast})^{2\alpha}}} - \displaystyle{{1}\over{2}}(S_{\rm e}+\varepsilon_{\ast})^{2\alpha},$$

where S _e is the equilibrium channel size (Appendix C). For fixed ε, we have that $\varepsilon_{\ast} \propto Q_{\rm in}^{-\textstyle{{1}/{\alpha}}} \to \infty$ as Q _in → 0. Equation (32) then implies that tr(J) < 0 for sufficiently small Q _in, meaning that the equilibrium becomes stable. This stability transition is a necessary condition for the supercritical Hopf bifurcation (Fig. 9 point i).

At the subcritical Hopf bifurcation (Figs 8b and 9 point ii), the equilibrium becomes stable because the meltwater supply increases past some critical value, as observed in the numerical results. The sign of tr(J) is equal to the sign of the function

(33)$${\cal T} = 2a_{\rm m} (1-N_{\rm e})\left[\alpha(1-N_{\rm e})^{1-({1}/{2\alpha})} - \displaystyle{{a_{\rm c}}\over{a_{\rm m}}} N_{\rm e}^n\right] - 1,$$

where N _e is the equilibrium effective pressure. For a constant ratio a _c/a _m and fixed ε_*, N _e and the quantity in square brackets in Eqn (29) are also fixed (Appendix C). Thus, varying a _c or a _m at a fixed ratio changes the sign of ${\cal T}$, resulting in stability transitions. In particular, increasing Q _in leads to decreasing a _m and, eventually, ${\cal T} \lt 0$. This condition implies that the equilibrium is locally stable, and, therefore, explains the damping with increasing meltwater supply. Likewise, ${\cal T} \gt 0$ implies that the equilibrium is locally unstable. We plot the sign of ${\cal T}$ in the a _m-a _c plane in Figure 10a and in the L _c-Q _in plane in Figure 10b. Physically, the equilibrium becomes stable because the filling timescale $\tilde {t}$ in Eqn (13) becomes small relative to the melting t _m and closure t _c timescales. Since the channel evolves too slowly relative to the meltwater supply rate, perturbations to steady state return to equilibrium. These stability transitions generalize to systems with arbitrary numbers of lakes (Appendix C).

Fig. 10. Sign of tr(J) evaluated at equilibrium for a range of a _m and a _c for the single-lake ODE model. The positive region (red) corresponds to an unstable equilibrium and the negative region (blue) corresponds to a stable equilibrium. We hold ε_* = 0.05 constant. For panel (b), we choose the parameters ψ₀ = 45 Pa m⁻¹ and A = 1 km² to compute ${\cal T}(a_{\rm m}(Q_{\rm in},\,L_{\rm c}),\,a_{\rm c}(Q_{\rm in},\,L_{\rm c}))$ in Eqn (33), where the equilibrium pressure N _e depends on a _m and a _c through solution of the non-linear Eqn (C5). Note that this stability diagram is for fixed ε_* so that the stability transition at small Q _in with fixed ε (Fig. 9 point i) is not observed.

For a range of Q _in near the subcritical Hopf bifurcation, there exists an unstable limit cycle around the stable fixed point (Figs 8b, c). In this region of parameter space, the asymptotic behavior of the solution depends on the initial condition. In particular, channels that are far from equilibrium approach the stable limit cycle. Intuitively, this difference in asymptotic behavior occurs because trajectories starting far from equilibrium pass through larger values of S and Q that admit higher melting and closure rates. Sustained oscillations are then possible since the channel evolution keeps pace with the filling rate. However, the melting strength a _m and closure strength a _c parameters are small enough to also admit solutions that converge to equilibrium. The equilibrium and stable limit cycle are very weak attractors in this parameter range. Convergence to the limit set can take more than one hundred years.

5. Discussion

Several natural lake systems exhibit behavior that is qualitatively similar to our numerical results. The Hidden Creek–Donoho Falls lake system in Alaska provides an example of a smaller downstream lake that is in phase with the upstream filling–draining cycle (Bartholomaus and others, 2008, 2011). Hidden Creek lake drains on a yearly cycle and the downstream lake (Donoho Falls) drains in phase with this, but also oscillates intermittently. Numerical results in Figures 4b–f show similar behaviors where the downstream lake has several small-amplitude oscillations between the main draining events that are in phase with the upstream lake. In this case, smaller downstream lakes may act as pressure gauges that are sensitive to upstream lake cycles (Anderson and others, 2005). Observations suggest that ‘through-flowing’ subglacial lakes may behave in a similar manner (Carter and others, 2013; Siegfried and Fricker, 2018).

The multiple-peaked downstream cycles explored here (Fig. 4) are similar in nature to the water pressure doublets studied by Schoof and others (2014), despite the difference in timescales and water storage mechanisms. The pressure doublets result from a distributed storage model forced with time-varying water input (Schoof and others, 2014). The similarity to the multiple-peaked lake cycles likely arises because the downstream lake in our model is forced by time-varying water input from the upstream channel. If the storage capacity is small relative to the filling rate, higher oscillation frequencies should be expected in both models.

Observations of lake systems in Antarctica also support our model. Delayed filling–draining events have been observed in East Antarctica (Wingham and others, 2006), the lower MacAyeal Ice Stream, and Slessor Glacier (Siegfried and Fricker, 2018). The out-of-phase oscillations depicted in Figure 5 are most similar to these observations. Out-of-phase oscillation requires the downstream lake to be larger than the upstream lake in the model. In the case of the MacAyeal Ice Stream, the downstream lake Mac1 is indeed slightly larger than the upstream lake Mac3 (Siegfried and Fricker, 2018). Similarly, filling of the larger lake on Slessor Glacier (Slessor₃) during 2005–2015 was coincident with the draining of the smaller adjacent lake (Slessor₂) (Siegfried and Fricker, 2018). However, in the case that these are in fact a single lake (Slessor₂₃), as suggested by Siegfried and Fricker (2018), this analysis is not applicable. An improved theory of how connected subglacial lakes behave could help clarify whether or not these are distinct lakes.

For similar-sized lakes, we have shown that double-peaked drainage events are possible in the downstream lake. Double-peaked events have been observed in a lake system at the confluence of the Whillans and Mercer ice streams in West Antarctica (Fricker and Scambos, 2009; Siegfried and others, 2016; Siegfried and Fricker, 2018). Between 2004 and 2008, two downstream lakes, Mercer Subglacial Lake (SLM) and L7, exhibited double-peaked filling–draining events that were roughly in phase with the upstream Subglacial Lake Conway (SLC). Following the first peaks in the downstream lakes, SLC drained and the downstream lakes drained again soon after. The filling–draining cycle in the Whillans–Mercer system repeated between 2010 and 2015 (Siegfried and others, 2016; Siegfried and Fricker, 2018). During this cycle, both SLM and L7 displayed double peaks again, although SLM's second peak was much smaller than in the first cycle. The SLM cycle is qualitatively similar to the alternating single–double peak cycles in Figures 6 and 7.

The alternative to oscillation is steady lake drainage. Our ODE model results suggest that oscillations may not be possible in cases where the lake filling rate is either much slower or faster than the channel evolution timescales (Fig. 9). An example of this behavior may occur at Lake Vostok, which has been inactive for at least a decade (Richter and others, 2014). However, Evatt and others (2006) simulated filling–draining cycles in Vostok-sized lakes with a model similar to the one used here. The oscillation period of Vostok was estimated to be at least 5000 years (Evatt and others, 2006). It is possible that Lake Vostok is so large that its filling rate would be too low for us to detect during quiescent, inter-flood periods.

While we see similarities between our numerical results and data from Antarctica, a complete model for subglacial lake drainage would include ice deformation over the lake. Inclusion of this would lead to a different functional relationship between the effective pressure and the lake depth, or different boundary conditions at the lakes (Evatt and others, 2006; Evatt and Fowler, 2007). This should result in oscillations that are closer in magnitude to observations of Antarctic subglacial lake activity. Ice flow also influences the evolution of the hydraulic gradient around the lake because filling and draining events lead to changes in surface slopes. Likewise, changes in ice thickness may also influence subglacial lake drainage events. We leave exploration of such models for future work but expect that the basic trends outlined here will hold.

Subglacial flooding in Antarctica has led to enhanced ice flow on several glaciers (Stearns and others, 2008; Scambos and others, 2011; Fricker and others, 2016; Siegfried and others, 2016). Ice deformation is usually related to effective pressure at the base through a sliding law. Understanding basal sliding in a drainage system during floods requires determining the evolution of effective pressure in all of the drainage elements. Realistic drainage basins include interacting channels, cavities, water sheets, canals, moulins and other drainage pathways (Creyts and Schoof, 2009; Schoof, 2010; Hewitt, 2011; Werder and others, 2013). Our single-lake results are similar to those obtained by Kingslake and Ng (2013a) who included a one-dimensional cavity model. Even so, it is unclear how a two-dimensional network of drainage elements would respond to oscillatory lake drainage. In the same vein, alternative hydrology models like the canal model of Walder and Fowler (1994) may lead to better quantitative agreement with observations of Antarctic filling–draining events (Carter and others, 2017).

An important question is whether infrequent large floods or frequent small floods lead to longer periods and higher magnitudes of enhanced sliding. Large lakes may temporarily shield the downstream drainage system from frequent upstream floods (Fig. 5). Counterintuitively, a larger upstream lake may also protect the downstream glacier from enhanced sliding by forcing downstream oscillations to occur less frequently (Fig. 4). However, a more comprehensive model would be required to test these ideas.

In our analysis, we developed an ODE model in Eqns (27–31) that agrees with the PDE model in Eqns (19–26) in the sense that it displays the same oscillation patterns as the meltwater supply rate, lake sizes and channel lengths are varied (Fig. 7). The simplifying assumptions and parameterizations used to develop the ODE model provide insight into the nature of the physical system. We showed that the ε-regularized Darcy–Weisbach law in Eqn (30) allows for the existence of limit cycles. In some sense, this regularization substitutes for the stabilizing effect of the background hydraulic gradient parameterization ψ(x) in Eqn (2) and channel source term μ in Eqn (4) that allow for filling–draining cycles in the PDE model. The existence of limit cycles in the ODE model suggests that oscillations are not uniquely produced by the particular background hydraulic gradient chosen here and elsewhere (Fowler, 1999; Kingslake, 2015). Furthermore, the agreement of the double-lake ODE results obtained with the ‘uncoupled’ hydraulic gradient approximation in Eqn (31) with the PDE results (Fig. 7) implies that the downstream pressure has little influence on the upstream system for the set of simulations considered here. We leave a more thorough comparison of these two models and alternative regularizations for future work.

6. Conclusions

Here, we introduced a generalized model for hydraulically connected ice-dammed lakes that extends previous work on single-lake systems (Fowler, 1999; Ng and Björnsson, 2003; Fowler, 2009; Kingslake and Ng, 2013a; Kingslake, 2015). Using this model, we explored how lake size and spacing influence the coupled depth-discharge oscillations in double-lake systems. We found that the downstream lake

• • oscillates in-phase with the upstream lake and exhibits multi-peaked cycles if it is moderately smaller ($A_2^{\ast} \lesssim 1$),

• • acts like a pressure gauge for the upstream lake if it is much smaller ($A_2^{\ast} \ll 1$),

• • oscillates out-of-phase with the upstream lake if it is moderately larger ($A_2^{\ast}\; \gtrsim\; 1$),

• • fills in a stepwise manner if it is much larger ($A_2^{\ast} \gg 1$), and

• • exhibits delayed and low-amplitude oscillations if it is far from the upstream lake ($L_1^{\ast} \gt 0.5$).

These oscillations are distinct from the single-lake oscillations and have features that are qualitatively similar to several natural lake systems. While we assumed that the meltwater supply rate to each lake was the same, other trends are likely to arise if the meltwater supply rates are different. We derived a simplified ODE model that is a good approximation of the PDE system near the lakes in the sense that it reproduces all of the same qualitative behavior. Using the simplified model, we explained how limit cycles appear and vanish as the meltwater supply varies. The simplified model is a useful tool for future work aiming to understand the complex coupling between subglacial lakes, drainage systems and ice-sheet flow.