The role of sediment transport in the mechanics of jökulhlaups

Abstract

The classical theory of jökulhlaups used Röthlisberger’s earlier theory of ice-channel drainage to describe the development of the flood hydrograph. This theory has some drawbacks: the mechanism of initiation (breaking the seal) is opaque, the Manning roughness coefficient is too large and the hydrographs can reveal a sudden switching from channel opening to channel closure which is not simulated by the model. In this paper, we examine these features by exploring a more detailed model, which takes into account the physics of sediment erosion and its effect on channel morphology. We propose a theory in which channels need not be semicircular, but have shapes determined by a local balance between closure and melting, and in which erosion of the tunnel margins is taken into account; in particular, we derive theoretical predictions for sediment discharge, and we also propose a mechanism whereby the pressure seal over the caldera rim at Grímsvötn in Vatnajökull, Iceland, can be broken when the lake-level water pressure is still some 6 bar below the maximum overburden ice pressure.

1. Introduction

Jökulhlaups are outburst floods which occur underneath glaciers due to the emptying of subglacial or ice-marginal lakes dammed by the overlying ice. The former type is common in regions of significant geothermal heat, such as Iceland, where one of the best-known examples is Grímsvötn beneath Vatnajökull. Figure 1 shows a longitudinal profile of the Vatnajökull ice sheet. The subglacial lake lies within a volcanic caldera, and the high geothermal heal flux melts the ice, increasing the lake volume. At some critical level, the water pressure bursts through the seal exerted by the ice dome at the caldera rim, and the resulting flood allows peak discharges of up to 10⁴ m³ s⁻¹ over a period of 2–3 weeks. These floods occur every 5–10 years; alter the lake drawdown during the flood, there is a slow replenishment of the lake level until the next outburst occurs.

Fig. 1. Cross-sectional view of Vatnajökull and the subglacial lake Grímsvötn along, the jökulhlaup drainage route (from Björnsson, 1974).

An initial model was proposed by Nye (1976), following Röthlisberger (1972). In this description, a semicircular channel in the ice is opened rapidly by the frictional heat generated by the water flow. The theory simulated well the rising limb of the hydrograph but more detailed work by Spring and Hutter (1981) indicated that the falling limb was less well explained. Clarke (1982) has also used the model.

In this paper we address three particular problems concerning the simulation of jökulhlaups, with special reference to the 1972 flood from Grímsvötn. The first problem concerns the initiation mechanism. The most obvious reason for an outburst is if the effective pressure N decreases to zero in the seal region. However, observations by Björnsson (1992) indicate that the Grímsvötn seal is broken when the effective pressure is of the order of 6–8 bar. The second problem is that in order to tit the rising limb of the Grímsvötn hydrography, Nye was forced to use a high value of the Manning roughness coefficient. A related and more serious problem is the third one, namely the inability of the Nye model to simulate the relatively abrupt shut-off observed in the typical Grímsvötn hydrograph.

2. A Wide-Channel Mechanism for the Grímsvötn Jökulhlaups

The total estimated suspended-sediment yield from the 1972 Grímsvötn jökulhlaup is 2.4 × 10¹⁰ kg (Tómasson and others, 1980). This corresponds to an average erosion of 180 m³ per unit length of the 50 km long subglacial drainage pathway. This strongly suggests that the flow channel is wide. Wide channels have been advocated by Hooke and others (1990). If the channels are wide, then we can expect their effective roughness to be larger, which could explain the anomalous Manning coefficient values determined by Nye (1976). In addition, we will show that wide channels can explain the sudden shut-down of the flood. In the following section, we propose a mechanism for breaking the seal based on Walder and Fowler’s (1994) wide-channel drainage theory.

3. Breaking the Seal

At the end of a jökulhlaup, the channel pressure is low and the leakage may fie of Darcy type in the seal region at the ice divide near the lip of the caldera. At lower altitudes a larger base flow will cause channels to form. We calculate an expected effective pressure by assuming a region of Darcy flow near the seal, and an R channel below the seal. The Darcy–Röthlisberger transition point is taken at 8 km from the lake/ice boundary. To account for the effect of surface-derived meltwater, we prescribe a water flux along the R channel which increases linearly with downstream distance s, from approximately zero at the transition point to the typical observed value (~ 100 m³ s⁻¹; from Björnsson, 1988) at the front portal.

Figure 2 shows the calculated effective-pressure profiles at pre-flood and flood-initiation lake levels. In the seal region, the effective pressure N has a maximum; as the lake refills, this maximum is decreased to a value of approximately 15 bar when the lake reaches its flood-triggering level at about 1430 m a.s.l. Notice that this differs from the value 6 bar which is based on hydrostatic equilibrium. The much higher effective pressure is due to the assumption of Darcy leakage.

Fig. 2. Computed raines of effective pressure N via Darcy–Röthlisberger flow at (a) the pre-flood, low lake level 1374.5 m, and (b) the typical flood-initiation lake level 1430 m.

The distinction between the Darcy and Röthlisberger regions is not only one of water flux. The ice-surface slope over the lower reaches is about 0.02, while at the seal it is smaller. The drainage theory of Walder and Fowler (l994) indicates that the style of drainage for water flow over deformable sediments depends strongly on the ice surface slope α. If α is very low, then Darcy flow is unstable for high enough water flux Q, and canals exist with N < , where is critical effective pressure: R channels are not feasible. At higher α (e.g. 0.02), both channels and canals are viable, and then channels are preferred.

Our mechanism for breaking the seal is as follows: in between floods, the Darcy flow in the seal has maximum N > , and no transition to a canal flow is possible; lower down, α is higher, so that channelized drainage can occur. However, as the lake refills, the maximum N in the seal region drops until N = , a canal is viable and we envisage a transition to a canal system in the seal. This transition effectively acts like opening a tap, and the flood proceeds; the details of the process await further study.

A Wide-Channel Model

Our wide-channel formulation is illustrated by the schematic diagram shown in Figure 3. It can be shown that the downstream variation of Q will normally be negligible, and hence that drainage during jökulhlaups can be ascertained by analyzing the opening and closure of the channel independently of downstream distance s. Let x be a cross-stream coordinate, and let h(x, t) be the channel depth. We suppose the channel is in −l < x < l, and in general l = l(t). In the Röthlisberger–Nye model, we effectively (and arbitrarily) assume l = h, and the relaxation of tins assumption requires an extra relation to describe the width l(t).

Fig. 3. Schematic diagram of the wide channel showing the mathematical symbols used in our formulation. The four thermomechanical channel processes of interest here are the vertical ice-melting and closure, and the lateral erosion and creep in flow of till sediments.

Our analysis is based on the limit h ≪ l; if h ~ l, then we regain the Nye model. If h ≪ l, we may approximate the ice-roof and bedrock boundaries as being parallel at every cross-stream position, then a local (in x) balance leads to the dimensional equations

(4.1)

where is the local melting rate, is the closure rate, û is the local flow velocity, τ is the local wall shear stress (given via the friction factor f*), and ϕ′ is the hydraulic gradient. Note that a friction parameterization which differs from the Manning formula in Equation(4.1) has been used in describing the local momentum balance; as a result, our formulation bears similarity to that of Walder and Fowler (1994); more specifically, since τ represents the total shear stress at the roof and bed boundaries (as opposed to only at the bed). f* can be related to the more familiar friction factor f in river mechanics by f* ≈ 2f.

In Equation(4.1) there are five equations for h, , , ϕ′, τ, p _w and û, so that we require two further equations. One is conservation of mass, which takes the global integrated form

. (4.2)

where L is the interval (−l, l), and M is a surface-derived meltwater supply rate. The section S and flux Q are given by

. (4.3)

The main difference between this and the Nye model lies in the closure relation between and N = p _i , – p _w, which is our second equation. We can adopt the methods of linear elasticity as used for crack problems by England (1971), for example, to show that, for constant-viscosity ice, the closure rate is determined by inverting

, (4.4)

where ϯ denotes a principal value integral, and μ _i, is ice viscosity; and if N is independent of x, then

. (4.5)

By analogs to closure of circular ice channels, we now extend this result and allow a non-linear flow law of Glen type by putting μ _i ⁻¹ = N ² μ ₀ ⁻¹. From Equation(4.1), we obtain

, (4.6)

and this has to be coupled with a lake inlet boundary condition relating N and Q. However, we see that l(t) is undetermined in the above model, and we have to prescribe l(t) independently to complete the model.

A Parameterized Model

We adopt a lumped parameter version of Equation (4.6) by writing

, (4.7)

putting

, (4.8)

where we omit the x dependence of the closure rate and the channel depth, so that variables are now functions of time only (i.e. ḣ ≡ dh/dt). It can be shown that the equation describing lake-water mass is

, (4.9)

where m _L is the lake-refilling rate, and, A _L. its effective area.

It remains to prescribe the half-width l(t). Here we suppose that bank closure occurs via inward creep of the till, while erosion outwards occurs via bank failure and sediment uptake by the water flow. We suppose bank erosion depends directly on stream power (the product τu), while inward creep of till may be inferred front related work by Fowler and Walder (1993). Our model for the evolution of l is then

(4.10)

where A _T, a, b are the Boulton and Hindmarsh (1987) flow-law parameters, and N _T is the pre-flood effective pressure in the till. The parameterization of the depth h also facilitates the calculation of a bulk channel roughness for comparison with the Nye model. We use the formula

, (4.11)

where for the present case R ≈ h/2, and remembering f ≈ f*/2, the effective Manning roughness for the wide channel is then

. (4.12)

We non-dimensionalize Equation(4.7), (4.9) and (4.10) by writing the variables h, N, l, t, Q, and |ϕ′| (which we call Φ) in terms of appropriately chosen scales h ₀, N ₀, l ₀, t ₀, Q ₀ , Φ₀, and we then have the dimensionless model

(4.13)

where

. (4.14)

Typical estimates for Grímsvötn lead to

, (4.15)

and then typical values of λ and μ are

. (4.16)

In Figure 4 we show the results of computing a flood hydrograph, where we have chosen parameter values to provide good fits to the rising limb and the maximum discharge of the Grímsvötn 1972 flood. These values are f* = 0.1, K ₁ = 1.6 × 10⁻⁴ kg m⁻⁵ s², Φ = 196 kg m⁻⁵ s⁻². We see that the model is able to simulate the hydrograph quite well, and in particular the sharp cut-off is mimicked. This is caused by the accelerated closure due to the rapidly increased value of l. The Manning roughness calculated from Equation (4.12) takes an average value of about 0.03 m^−1/3 s, more normal than the value inferred from Nye’s (1976) work. Finally, the simulation indicates a maximum channel width of about 90 m. If a layer of sediment 2 m thick is eroded in the flood, then we can also simulate the estimated figure of 180 m³ sediment eroded per unit length of channel.

Fig. 4. A simulated hydrograph (solid line) together with hydrographic data points (crosses; from Rist, 1973) of the 1972 Grímsvötn jökulhlaup: (a) linear discharge scale; (b) logarithmic discharge scale. Dashed lines are the simulated recession solution of Equation (4.19).

It will be seen that we have improved the estimate of maximum discharge at the expense of producing a very sharp fall-off in water flux. This can be rectified by incorporating the spatial dependence of the falling discharge. The mass conservation law corresponding to Equation (4.2) can be written dimensionlessly as

, (4.17)

with the base flow melt run-off being estimated as M ≈ 0.01 (corresponding to 100 m³ s⁻¹ outflow at the terminus). The point here is that, even if the inlet flow shuts off instantaneously, an expansion or rarefaction wave will propagate down the channel, so that h (and hence water flux) will shut off more gradually. Indeed we can find an explicit solution based on the method of characteristics for this decay (Carrier and Pearson, 1976).

At flood termination, l changes little, so the shutdown in h is described (ignoring the melting term) by

, (4.18)

where = M/l 10⁻² (since l ~ O(1)). Suppose that h = 1 for s > 0 at t = 0, and then h = 0 for s = 0 at t > 0 (this corresponds essentially to the results obtained so far). The solution; u the outlet s = 1 is then h = 1 + ( t/ε) for 0 < t < 2ε/3, when the expansion wave reaches the outlet, after which

, (4.19)

until t. = ε/ ^1/3, when the base flow with h = ^2/3 resumes.

We have fitted the descending limb of the hydrograph using Equation (4.19) with values of = 2.5 × 10⁻³ and ε = 2, as shown in Figure 4. Although we successfully fit the hydrograph, the value of ε is abnormally high (an appropriate value is ε = 0.03). We consider that this may be due to a distributed stream pattern during floods, which effectively may increase the tortuosity of the channels by the requisite amount.

5. Conclusions

If the restriction of semicircular channels in flood models is removed, then we need separate equations for channel depth h and half-width l. We propose an ice-closure equation for the depth, and suppose that width is determined primarily by erosion and inward creep of sediments. As a result, we are able to provide an improved simulation hydrograph of the 1972 Grímsvötn jökulhlaup with a plausible value of Manning roughness.

In order to explain the flood-trigger mechanism when the lake pressure is some 6–8 bar below the maximum ice-overburden pressure, we invoke the fledgling theory of basal drainage through subglacial sediments, and suggest that Darcy flow in the seal region becomes unstable when the effective pressure becomes sufficiently low, leading to canalised flow and subsequent flood conditions. A more precise description awaits further work.