2.1. Channelized drainage

Reference NyeNye (1976) envisaged the channel cross-sectional area SR to evolve under the competition between melt enlargement and creep closure. We use Reference FowlerFowler’s (1999) version of his model, which allows the discharge QR to go negative (e.g. when the channel drains into the lake between floods), and where the equations for channel evolution, mass conservation and momentum conservation, are, respectively:

(1)

(2)

(3)

In Eqn (1), the two terms on the right represent melting and closure. The mass melted per unit channel length per unit time, m, is given by m = ( + <9NR/<9s)|QR|/L; this expression assumes that potential energy is instantaneously dissipated through friction to melt the channel walls. NR is channel effective pressure (=Pi-PR, with Pi being ice overburden pressure, PR being channel water pressure), pw and pi are water and ice densities, g is gravity, L is latent heat of fusion of water and f is hydraulic roughness (f= n’2(l2/SR)2/3 ?s0.07 m–2/3 s2 for a semicircular channel with wetted perimeter land Manning’s roughness coefficient n’=0.1 m–1/3s). The constants n = 3 and K0 = 10–24 Pa-3 s–1 pertain to temperate ice rheology (Reference FowlerFowler, 1999); K0 accounts for the assumed semicircular cross section of the channel.

The term Tin Eqn (2) represents the rate of water transfer into the channel from the cavities (volume flux per unit distance). Following Reference Hewitt and FowlerHewitt and Fowler (2008) we model it as proportional to the pressure difference between these drainage components:

(4)

where NC is the cavity effective pressure and k is a connectivity constant (we adopt Hewitt and Fowler’s value 10–9m2 s–1 Pa-1). R k is a dimensionless factor used in our experiments (Section 3.1) for varying the transfer rate Rkk that saves us from writing out its units each time.

Equation (3) derives from Manning’s friction formula and is the basic hydraulic gradient, given by

(5)

where 4>b is the glacier bed slope. We assume a with uniform thickness H (Fig. 1), so is constant.

Equations (1-3) neglect the lake’s internal energy that could enhance the channel melt rate m (Reference ClarkeClarke, 1982), so this model assumes lake and subglacial water temperatures equal to the pressure-melting point. Also, as NR and NC stay positive in our simulations, we do not model channel enlargement caused by ice uplift when NR

2.2. Lake evolution

The lake depth hL (Fig. 1) evolves with water input (Qin) and outflow (QR at the channel inlet) following the continuity equation

(6)

where AL is uniform lake area. Equation (6) expresses condition (1) of our full coupling (Section 1). According to condition (2) of our full coupling, the hydrostatic lake pressure fixes the channel-inlet effective pressure to be

(7)

at all times, and this provides a boundary condition for NR in Eqn (3).

2.3. Linked-cavity drainage and basal sliding

We follow Reference Hewitt and FowlerHewitt and Fowler (2008) in using a pseudo-steady description for the geometry of the cavities, with a conservation equation for the water flux through them and a sliding law. For sliding at velocity ub over bed obstacles of height R, Reference WalderWalder (1986) showed that the typical cavity cross-sectional area is

(8)

where γ ≈ 0.32 and k₁ ≈ 1.1 are constants, and that the water discharge through the cavity system, Q_C , depends on hydraulic gradient (just as for the channel) so that Q_C and S_C are related via

(9)

We adopt a power law to describe sliding,

(10)

(e.g. Reference BindschadlerBindschadler, 1983), where p = 4, q=1, cs is a constant dependent on bed roughness ( « 2 x 10–20 ms–1 Pa-3; Reference Hewitt and FowlerHewitt and Fowler, 2008), and Tb is the basal shear stress assumed to balance local driving stress pigHsincj)s where 4>s is glacier surface slope (Reference Cuffey and PatersonCuffey and Paterson, 2010). Following Reference WalderWalder (1986) and Reference KambKamb (1987), we assume that the sliding control on cavity size (the first of the two terms enclosed by the square brackets in Eqn (8)) dominates frictional heating (the second term). On neglecting the latter term, eliminations between Eqns (8) and (10) lead to

(11a)

and

(11b)

which give the cavity effective pressure and sliding velocity in terms of SC. For the bed obstacle height R, we use Reference WalderWalder’s (1986) value of 0.1 m.

A final equation describes water conservation in the cavities:

(12)

where MC is the background water supply to the cavities from basal- and surface-derived melt. We assume that the channel receives meltwater only via the cavities (through T) and none directly via moulins or other englacial routes.

1Equations (-6), (9), (11a), (11b) and (12) complete our model. In essence it describes the coupled evolution of channel size, cavity size, and lake level. SR, QR, NR, T, SC, QC, NC and ub are functions of t and s, while hL is a function of t only. This model differs from Hewitt and Fowler’s not just in its inclusion of the lake and a channel that evolves at a rate dictated by melt and ice creep, but also in using the total hydraulic gradients > +<9NR/<9s and + dNC/ds to drive channel and cavity flow (rather than alone). Four boundary conditions are needed in Eqns (2), (3), (9) and (12) for determining the discharge and pressure in these systems. For the channel we impose the condition in Eqn (7) and NR = 0 at the terminus (s = s0). For the cavities, an equation like (7) could be used to couple the cavities and the lake hydraulically. Preliminary investigations indicate that interesting dynamics result from this coupling, but a full exploration of its consequences is beyond our scope here. Instead, for simplicity, we isolate the cavities from the lake, by assuming that NC at the lake (s = 0) and at the terminus (s = s 0) are both constant, «1 x 105 Pa. Although this value may seem arbitrary, the qualitative features of our results reported below are not sensitive to its size.

Our assumption of pseudo-steady cavity geometry is motivated by the supposition that, unlike in channels, the cavity system’s effective pressure-discharge relationship will be negative whether the system is steady or not, and motivated by the lack of unsteady cavity drainage models that are strongly verified by subglacial observations. However, the potential limitations of our model need to be recognized. Since, in reality, cavity size will not react instantaneously to pressure changes, our model overestimates the role of water storage in cavities when water transfer changes rapidly and neglects the lag of peak cavity size behind peak cavity water pressure. In such situations, the model will underestimate actual NC and ub variations. Also, our model precludes investigation of a possible negative feedback between cavity size and pressure suggested by Reference Bartholomaus, Anderson and AndersonBartholomaus and others (2011). In a recent theory, Reference Hewitt, Schoof and WerderHewitt and others (2012) have tried to overcome these limitations by modelling the cavity system dynamically as a sheet with a cavity evolution equation analogous to Eqn (1).

2.4. Method of numerical solution

We solved the equations with the finite-difference method on a discretized space-time domain. The forward Euler method was used to integrate Eqns (1), (6) and (12), with upwinding for the spatial derivative in Eqn (12). When calculating results at the next time-step, we first used variables from the current time-step to update the lake depth hL and cavity area SC, from which the sliding velocity ub, cavity discharge QC and pressure NC could be found by Eqns (11) and (9). Turning to the channel, we then used Eqn (1) to update SR, and used SR in Eqns (2) and (3) (with their boundary conditions and NC at the new time) to find the discharge and pressure distributions QR and NR at the new time-step. In this boundary-value problem, where (and whether) QR switches sign along the channel is not known a priori. To overcome this we solved Eqns (2) and (3) by two nested iteration loops. In the inner loop, a guess of the lake outflow QR(s = 0) at the new time was used as the boundary condition for integrating Eqn (2) into s>0. We used dynamic relaxation for Eqn (2) and introduced a fictitious sdQR/dt to its left-hand side, and iterated QR to steady state on the fast (e) timescale. This inner loop thus yields the QR and NR distributions simultaneously, including the channel effective pressure at its inlet, NR(s= 0). Generally, QR(s = 0 ) and NR(s=0) at the new timestep from this loop mismatch the lake condition in Eqn (7). In the outer loop, we use the Newton-Raphson method to correct the guess of QR(s = 0) iteratively (letting the inner loop converge each time) to obtain the match.

3.1. Flood cycles

Nye’s model formulated with only time dependence and no spatial dependence can simulate recurring floods, but their size grows indefinitely with time through a feedback between decreasingly small channel size and increasing highstand and peak discharge (Reference NgNg, 1998; Reference Ng and BjörnssonNg and Björnsson, 2003). In reality many Jökulhlaups terminate through complete emptying of the lake, which would interrupt such unbounded growth, but Reference Ng and BjörnssonNg and Björnsson (2003) also found this behaviour when simulating floods from Grímsvotn, Iceland, which does not empty (Björnsson, 2002). Consequently, unbounded flood growth is seen as an intrinsic weakness of some Jökulhlaup models. Reference FowlerFowler (1999) showed that model flood cycles remain bounded if spatial dependence is included, the channel receives a uniform meltwater supply and the basic hydraulic gradient near the lake is negative (e.g. due to retrograde ice-dam surface slope towards the lake - a ‘topographic seal’). Whether uniform channel supply and a topographic seal are necessary conditions for bounded flood cycles has remained unexplored. As our model treats how the channel captures subglacial water in more detail than Fowler’s model and uses a positive basic hydraulic gradient everywhere, we address this by asking whether we can simulate stable flood cycles.

Figure 2 shows the simulated time series from a set of model runs, each spanning 15 model years and conducted with an initial lake depth of 30 m (hL at t=0), a constant lake-water input Qin=10m3s–1, a constant background water supply to cavities M C = 1 x10–3m2s–1 , and several values of the parameter Rk between 0 and 1. Rk = 0 and Rk=1 refer to total hydraulic isolation and the hydraulic coupling used by Reference Hewitt and FowlerHewitt and Fowler (2008), respectively. The coupled lake-channel-cavity system oscillates in time, yielding asymmetric filling and draining cycles of lake depth (Fig. 2a, c, e and g) and corresponding floods in the channel discharge (Fig. 2b, d, f and h). In all runs, initially, the volume and peak size of successive floods increase with time in a transient response to initial conditions. When Rk<0.3, the lake’s post-flood ‘lowstand’ level decreases until the lake empties completely, interrupting the simulation (Fig. 2a-d). In contrast, when R k > 0 . 6 , the cycles become periodic and stable without emptying the lake (Fig. 2e-h). Runs using 0.3 <R k < 0.6 (not shown) show that increasing Rk delays the emptying of the lake as peak discharge is suppressed by increased cavity-to-channel water transfer.

Fig. 2. Modelled lake level, hL(t) (left), and channel discharge at the lake, QR(0, t) (right), for Rk= 0, 0.3, 0.6 and 1. Cavity water supply, MC =1x10–3 m2 s–1.

These results go beyond Fowler’s to show that stable limit cycles occur when high T values allow the channel to capture water efficiently from the cavities and a water divide (or seal) to form between floods. The stabilizing effect of this transfer is evidenced also by the impact of Rk on the cycles: the higher is Rk, the shorter is their period and the smaller the flood volumes. Furthermore, we show that the channel can drain water into the lake at low lake levels and stable flood cycles can occur even if the basic hydraulic gradient is positive everywhere: a topographic seal is not a necessary condition for water divide formation.

3.2. Spatio-temporal evolution of sliding and drainage

The periodicity of stable cycles provides us with a common basis for studying individual floods and how their characteristics depend on model parameters, without the obstacles associated with initial conditions. Here we analyse the spatio-temporal behaviour of subglacial drainage and ice flow in the («1.5year long) flood cycle highlighted by the boxes in Figure 2g and h.

Figure 3 shows the simulated lake depth hL(t) and lake outflow QR(s = 0, t) (Fig. 3a), basal sliding velocity ub(s,t) (Fig. 3b) and rate of cavity-to-channel water transfer T(s, t) (Fig. 3c) through the cycle. Figure 3c is a map of Tover time and distance; a similar map is also used to plot ub(s, t) in Figure 4c, studied in Section 3.3. The hydrographic sequence here has been discussed by Reference FowlerFowler (1999). In Figure 3b, the time before point A (and after point E of the last cycle) is the lake-refilling phase between successive floods. At point A, high lake level moves the subglacial ‘seal’ (water divide in the channel) back to the lake, initiating the flood and allowing Nye’s positive feedback to enlarge the channel. B marks the time of the lake-level highstand when the flood’s growing discharge instantaneously matches the lake-water input. C and D mark the flood peak and the post-flood lake lowstand, respectively.

Fig. 3. Evolution of (a) lake level, hL(t), and channel discharge at the lake, QR(0, t); (b) sliding velocity, ub(s, t); and (c) cavity-channel water transfer rate, T(s, t) in the limit cycle indicated by the boxes in Figure 2g and h. MC=1 x 10–3 m2s–1, Rk = 1.

Fig. 4. Sliding velocity ub(s, t) (filled contour maps) and channel discharge at the lake QR(0, t) (white dashed lines and right-hand vertical axis) in one limit cycle for (a) MC=2 x 10–4m2s–1, (b) MC=4 x 10–4m2s–1, (c) MC=1 x 10–3 m2s–1 and (d) MC=2 x 10–3m2s–1.

Regarding ice motion, Figure 3b shows similar responses in ub at different positions on the glacier. ub rises gradually as the lake fills, but accelerates rapidly after flood initiation (A) to peak soon after the time of the lake highstand (B). ub drops during the main flood (B to D) to a much lower velocity than during lake refilling, then increases again after flood termination (D to E). The durations of elevated and depressed ub are comparable. Spatially, ub rises monotonically with distance down-glacier, but its time variations at different positions are not precisely synchronous. Figure 3b shows that the peak and trough in ub propagate down-glacier at ?s250md–1 and «450md–1, respectively (see the dashed lines in Fig. 3c).

This sliding response can be explained by considering our model’s physics and is the direct result of the full coupling between the lake and the drainage system. Equation (11b) shows that ub varies directly with cavity size, SC. In turn, SC depends on how much water the cavities gain from the background supply MC and lose to the channel via water transfer, whose rate T increases with channel effective pressure, NR (Eqns (12) and (4)). The sliding response thus stems ultimately from changes in lake level, which, through its effect on the channel pressure, governs how fast the channel draws water from the cavities. The timing associations in Figure 3 between hL, T and ub confirm this explanation. When the lake fills between floods (E to A), reducing NR at the channel inlet, its hydrostatic pressure is not transmitted past the seal (located in Fig. 3c by the black line) to modulate Tand ub directly far down-glacier (further simulations show that this effect would be diminished if the cavities and lake were coupled hydraulically); there is a slight increase in ub only because the lake level weakly affects the distributions of NR and Tin the channel through its control on the migrating seal position. However, after the flood starts at t« 11.4 years, the seal’s absence means that the pressure transmission is unimpeded and able to reduce NR and T markedly down most of the channel (see descent into the dark trough in Fig. 3c); this causes the cavities to expand and accelerates sliding (A to B; 11.4< t< 11.7years). Later, as flood discharge grows, peaks, and recedes under Nye’s mechanisms (B to D, 11.7 < t < 12 years), the lake level drops rapidly, raising NR along the channel so that enhanced water transfer out of the cavities shrinks their size and decelerates sliding. After the seal re-forms at flood termination (E, t«12.2 years) we enter a new cycle. During the flood when the lake strongly modulates the cavity-to-channel water transfer, the nonlinear wave properties of Eqns (1-3) determine the propagation of changes in Tand ub down the channel.

3.3. Sensitivity to the cavity background supply, MC

Finally, we study how the ice-flow and flood evolution depends on the supply MC, which is the other control on cavity drainage beside T. For temperate glaciers, MC likely encompasses water reaching the bed from the surface, so it can be viewed as a proxy for weather conditions.

Figure 4a-d plot the maps of ub(s, t ) and time series of lake outflow from four model runs that assumed different values of MC between 2 x10–4 and 2x10–3m2s–1. Each plot’s duration is one flood cycle. The third run (Fig. 4c), identical to that discussed in Section 3.2, serves as the control. An enduring pattern in these runs is that the ice flow speeds up after flood initiation and begins to decelerate before the flood peaks. In the run with the lowest MC (Fig. 4a), the speed-up is very sudden and a strong compression wave (dub/ds<0) propagates down-glacier some 70 days prior to the flood peak (around 12.5 < t< 12.7 years).

Physically, we expect that decreasing MC reduces cavity size and discharge via mass conservation (Eqn (12)), and hence reduces the sliding velocity ub (Eqn (11a)). We see this trend in Figure 4a-d (note their different scales) as we move from right to left; however, the percentage variation in ub in the cycles increases. Figure 5 quantifies both trends by plotting the mean, maximum and minimum values of ub at 4 km from the lake in each cycle, including results from extra runs where MC is varied within the same range. As MC is raised, the mean value of ub increases but its maximum-to-mean ratio decreases.

Fig. 5. Mean sliding velocity at 4 km from the lake, ub (circles and left axis), and peak lake discharge, QPK (crosses and right axis), for MC =2x10– 4 m2 s–1 to 2 x 10–3 m2 s–1. The top and bottom ends of the vertical lines indicate maximum and minimum ub at 4 km from the lake.

We expect two outcomes of decreased MC on cavity-channel interaction. One is increased dominance of the transfer term in the cavity mass conservation equation (Eqn (12)). This explains the variable maximum-to-mean ub ratio; when MC is large, the remaining, time-varying terms in Eqn (12), dQC/ds and T, are comparatively small and SC and ub vary less than when MC is small; then, SC and ub are modulated strongly because Tdominates the mass balance. In addition, when MC<4 x 10–4m2 s–1, the cavities collapse after the lake reaches its highstand due to high effective pressure in the channel, preventing significant reduction of ub below its mean value (Fig. 5); then, most of the background water supply to the cavities reaches the channel, and the cavity area is controlled by MC=k(NR-NC) and Eqn (11a). In the case where MC = 2 x 10–4m2 s–1, the cavity system does not recover until after the divide reaches the lake (Fig. 4a; t« 12.3 years). These model results imply that under cold weather conditions when subglacial water is less abundant, glacier flow velocities would be lower during quiescence but respond more strongly in absolute and relative magnitudes to subglacial floods.

The other expected outcome of decreasing MC is that the cavity effective pressure NC increases (Eqn (11a)), reducing water transfer Tto the channel (Eqn (4)). Analysis of channel equations (1-3) shows that the seal then migrates more slowly towards the lake as it refills; this delays flood initiation, so the lake reaches a higher highstand - a higher flood initiation threshold, to cause a higher flood peak discharge, QPK. Figures 4 and 5 show that the negative dependence of QPK on MC is nonlinear and becomes sensitive only for MC< 6 x 10–4 m2 s–1.