Glacier Outburst Floods From “Hazard Lake”, Yukon Territory, and the Problem of Flood Magnitude Prediction

Abstract

In August 1978 “Hazard Lake” released 19.62 × 10⁶ m³ of water through a subglacial tunnel beneath Steele Glacier, Yukon Territory, Canada. The discharge during the outburst flood was measured by recording lake level changes with time, and a peak discharge of approximately 640 m³ s^–1 was estimated from the data. We have attempted to model the 1978 flood from “Hazard Lake” using an adaptation of Nye’s (1976) theoretical model for jökulhlaups from Grimsvötn. Our aim has been to calibrate the Nye model as a first step toward using it as a peak discharge estimator for other glacier–dammed basins. The agreement between our measured and simulated hydrographs is good, and we find that creep closure, though included in our analysis, appears to play an insignificant role in limiting the discharge of “Hazard Lake”. Release of thermal energy from the relatively warm lake water is the dominant factor contributing to tunnel enlargement.

The Manning roughness of outlet channels from glacier–dammed lakes is not known a priori and must either be assumed or estimated after the fact from the flood hydrograph. For “Hazard Lake” our fit implies Manning roughness in the range n′ = 0.105 m^–1⁄3 s, consistent with Nye’s estimate of n′ = 0.12 m^l⁄3s for the 1972 Grimsvötn flood and our estimate of n′ = 0.12 m ^1/3 s for the 1967 Summit Lake flood. If the Manning roughness for flood conduits can be shown to lie within a narrow range, this would constrain one of the least certain variables of the Nye model.

By making several simplifying assumptions, we have succeeded in reducing our adapted version of Nye’s model to a simple mathematical description involving dimensionless numbers characterizing reservoir geometry and the relative magnitudes of creep closure and tunnel enlargement by melting. In this simplified form, the influence of lake temperature, reservoir geometry, and creep closure on the character of flood hydrographs can be conveniently studied.

Résumé

Les debâcles glaciaires à partir du “Hazard Lake” dans le territoire du Yukon et le problème de la prévision de l’ampleur de la crue. En août 1978 le “Hazard Lake” lança 19,62 × 10⁶ m³ d’eau à travers un chenal sous–glaciaire sous le Steele Glacier dans le territoire du Yukon au Canada. Le débit pendant la crue de débâcle a été mesuré en enregistrant les variations dans le temps du niveau du lac et on a pu estimer la pointe de crue à approximativement 640 m³ s d’après ces données. Nous avons essayé de “modéliser” la crue de 1978 du “Hazard Lake” en utilisant une adaptation du modèle théorique de Nye(1976) pour les jökulhlaups de Grimsvötn. Notre but a été de caler le modèle de Nye à litre de premier pas vers son utilisation comme un estimateur du débit pour d’autres bassins avec barrage glaciaire. La concordance entre nos hydrogrammes mesurés et simulés est bonne et nous trouvons que l’obturation par le glissement du glacier, bien qu’incluse dans notre analyse semble jouer un rôle insignifiant pour limiter le débit du “Hazard Lake”. La libération d’énergie thermique à partir de l’eau relativement chaude du lac est le facteur dominant contribuant à l’élargissement du tunnel.

Le coefficient de rugosité de Manning des canaux émissaires des lacs proglaciaires n’est pas connu a priori et doit être supposé ou estimé après coup à partir de l’hydrogramme de la crue. Pour le “Hazard Lake” notre ajustement implique un coefficient de Manning de l’ordre de n′ = 0,105 m^–1⁄3 s, cohérent avec l’estimation de Nye de n′ = 0,12 m^₋1/3 s pour la crue de Grimsvötn en 1972 et notre propre estimation de n′ = 0.12 m ^–1⁄3 s pour la crue du Summit Lake en 1962. Si on peut montrer que les coefficients de Manning pour les chenaux de vidange restent dans un faible domaine de variation cela réduirait l’incertitude sur l’une des variables les moins stables du modèle de Nye.

Grâce à plusieurs hypothèses simplificatrices, nous avons réussi à réduire notre version adaptée du modèle de Nye à une simple description mathématique comprenant des nombres sans dimensions pour caractériser la géométrie du réservoir et les ordres de grandeur relatif de l’obturation par le glissement et de l’élargissement du tunnel par la fusion. Sous cette forme simplifiée, l’influence de la température du lac de la géométrie du réservoir et de l’obturation par le glissement sur les caractéristiques des hydrogrammes de crue peut être convenablement décrite.

Zusammenfassung

Glaziale Ausbruchfluten aus dem “Hazard Lake”, Yukon Territory, und das Problem der Vorhersage von Ausbruchsmengen. Im August 1978 entsandte der “Hazard Lake” 19,62 × 10⁶ m³ Wasser durch einen subglazialen Tunnel unter dem Steele Glacier, Yukon Territory, Kanada. Der Ahfluss während der Ausbruchsflut wurde durch zeitliche Beobachtungen des Seespicgels gemessen: aus den Daten ergab sich eine Abflussspitze von ungefähr 640 m³ s ^–1. Wir versuchten, die Flut von 1978 aus dem “Hazard Lake” mit Hilfe einer Anpassung von Nye’s (1976) theoretischem Modell für die Gletscherläufe von Grimsvötn zu simulieren. Unser Ziel war, das Modell von Nye zu kalibrieren, als ein erster schritt zu seiner Benutzung als schutzmittel für den Gipfelabfiuss anderer cisgedämmter Seen. Die Übereinstimmung zwischen unseren gemessenen und simulierten Pegelstanden ist gut. und wir können feststellen, dass der Abschluss durch Kriechvorgänge, der zwar in unserer Analyse berücksichtigt war. eine nur unwesentliche Rolle bei der Beendigung des Abflusses aus dem “Hazard Lake” zu spielen scheint. Das Freiwerden thermischer Energie aus dem relativ warmen Seewasser ist der wesentliche Faktor für die Vergrösserung des Tunnels.

Die Manning–Rauhigkeit von Abflusskanälen gletschergedämmter Seen ist nicht a priori bekannt und muss entweder angenommen oder nach der Flut aus den Abflusswerten geschätzt werden. Für den “Hazard Lake” ergeben unsere Auswertungen eine Manning–Rauhigkeit im Bereich von n′ = 0,105 m^–1/3 s, was mit Nye’s Abschätzung von n′= 0,12 m^–1⁄3 s für die Grimsvötn–Flut von 1972 und unserer eigenen schätzung von n′ = 0,12 m^–1⁄3 s für die Surnmit Lake–Flut von 1967 übereinstimmt. Wenn sich zeigen lässt, dass die Manning–Rauhigkeit von Abflusskanälen innerhalb enger Grenzen liegt, würde eine der unsichersten Variablen in Nye’s Modell erfasst sein.

Durch einige vereinfachende Annahmen gelang es uns, unsere angepasste Version des Nye–Modells auf eine einfache mathemalische Beschreibung zu reduzieren, in der zwei dimensionslose Zahlen die Geometrie des Beckens und das Grössenverhältnis des Kriecliabschlusses und der Tunnelerweiterung durch schmelzen kennzeichnen. In dieser vereinfachten Form kann der Hinfluss der Seetemperatur, der Beekengeometrie und des Kriechabschlusses auf den Charakter des Flutablaufes bequem untersucht werden.

Symbols

Introduction

The 1965–66 surge of Steele Glacier, Yukon Territory, Canada, displaced ice as much as 9.5 km (Stanley, 1969). One result of the surge was that the normal stream channel of Hazard Creek became ice–dammed, forming a proglacial lake now referred to as “Hazard Lake”. The lake probably began to fill in June 1966 and reached its maximum level between September 1966 and August 1967. Maximum water level is controlled by a marginal spillway over bedrock and water cannot rise above 1674 m a.s.1. Once filled, the lake stayed full until late July 1975 when it drained through a subglacial tunnel. The tunnel remained open and the lake basin empty through the summer of 1976. During the winter of 1976–77, the tunnel resealed. The lake again filled to its maximum level by July 1977, then drained subglacially around 2–5 August 1977. Since that time the lake has continued to fill and drain on an annual cycle (Table I). A description of the history and drainage of “Hazard Lake” covering events up to August 1977 is given by Collins and Clarke (1977). Several small discrepancies between that paper and this one arise because new information has become available and minor errors have been corrected. (For example, an aerial photograph by W. A. Wood was incorrectly dated and there is now no evidence for a filling between 5 August and 2 September 1977.)

Table I. Chronology for “Hazard Lake”

Research on “Hazard Lake” began in 1974 when the lake was surveyed and a bathymetric map prepared (Collins and Clarke, 1977). Our present interest in the lake arises from its 1975 transition from stable to cyclic drainage, and from its suitability as a field example against which to test a model of jökulhlaups proposed by Nye (1976). We have monitored the 1978 and 1979 discharges of the lake by measuring water level as a function of time. Knowing the geometry of the lake, it is a simple matter to convert level changes to volume changes, and thus calculate discharge rate. This makes it possible to compare the observed discharge to predictions of Nye’s model (Nye, 1976). Although his model successfully matches the discharge from Grimsvötn in Iceland, it is not known how widely applicable it is. In particular, there are variables such as the Manning roughness, which are not directly observable. For this reason the Nye model is uncalibrated.

It would be helpful if the model could be calibrated by testing what assumptions are necessary to match the discharge curves for a large number of jökulhlaups. If, for example, the Manning roughness (as empirically determined by fitting Nye’s model to a number of discharge records) proved relatively constant, the model might become useful for predictive simulations. This would be valuable for bridge and pipeline engineers since it would allow estimation of peak discharge from glacier–dammed basins. At present the only available tool is the purely empirical Clague–Mathews formula (Clague and Mathews, 1973)

where Q _MAX is measured in m³ s^–1 and V ₀ in m³. The formula is simple and appears to work reasonably well, but, if the Nye model is accepted, additional variables such as tunnel slope would be expected to affect peak discharge.

Although the floods associated with the drainage of “Hazard Lake” are not destructive, they provide an opportunity to test and attempt to calibrate the Nye model. Once calibrated, if this proves possible, the model might be used with some confidence to predict hydrographs and peak discharge from glacier–dammed lakes. An immediate practical application would be to assess the danger posed by outburst floods to planned developments along the Alaska Highway, Yukon Territory. There are more than 150 basins in the White, Donjek, and Kaskawulsh–Slims River drainage systems which are, or have been glacier–dammed (Canada, unpublished). Of these, three have estimated volumes exceeding 100 × 10⁶ m³ and could be expected seriously to threaten down–stream bridge and pipeline crossings. It would be helpful if further tools besides the Clague-Mathews formula were available to evaluate these hazards.

Transition to Cyclic Drainage

The transition, in 1975, from stable to cyclic drainage requires comment. From the initial filling in 1966 to late July 1975, the lake level was controlled by a marginal spillway over bedrock, and outburst floods did not occur. Since the first tunnel–drainage event in July 1975, the lake has filled and drained cyclically with annual or biennial periodicity (Table I). The simplest explanation for such a transition is given by Thorarinsson (1939): subglacial drainage from the lake becomes possible when the hydrostatic pressure p _i of water from the lake exeeds ice overburden pressure p, in the region of the "seal" (Fig. 1). Björnsson ([1975]) and Nye (1976) favour such a mechanism for controlling the drainage of Grimsvötn in Iceland, and we think it applies to “Hazard Lake” as well. There is evidence that the Grimsvötn drainage starts with p = ρ _w gh _w somewhat less than p _i = ρ _i gh _i, and Nye explains how this could be possible. The discrepancy is not large, and we shall simplify discussion by assuming that equality of p and p _i is the condition for drainage onset; this amounts to assuming h _w = ρ _i h _i/ρ _w at the start of the outburst flood.

Fig 1 A schematic cross–section (roughly to scale) showing “Hazard Lake” and the Steele Glacier ice dam. The subglacial topography and actual location of the seal are unknown.

Our explanation of the transition to cyclic drainage is that before July 1975 ice pressure always exceeded water pressure, because the marginal spillway prevented water from rising to the value necessary for flotation of the dam. Since the end of the 1965–66 glacier surge, ice ablation has steadily lowered the dam height and flotation is now possible when the lake is full. These ideas are given quantitative support in the following discussion, but the question cannot be finally settled until an ice thickness survey of the glacier is performed.

An unpublished map, based on August 1966 aerial photography and 1967 ground control, shows the August 1966 elevation of the Steele Glacier surface near the dam to be 1700–1720 m a.s.l. Surveys by S. G. Collins give the maximum lake elevation as 1 674 m a.s.l. and the elevation at the base of the ice dam near the tunnel inlet as 1574 m a.s.l. The Steele Glacier surge ended by summer 1967, and since then ice flow rates near the dam have been negligible, less than 2.0 m a^–1 according to surveys by Collins in 1974 and 1975. From the same surveys, the 1974–75 ice ablation rate near the dam was 2.65 m a^–1. The typical 1975 elevation of the dam was roughly 1705 m a.s.l. If one accepts the flotation hypothesis, the seal lies 270 m below high–water level, implying ice thickness of 300 m in the seal region. Although we have no sounding data, this depth is at least plausible.

Flotation may be an acceptable triggering mechanism for “Hazard Lake” jökulhlaups, but there is no question that jökulhlaups can occur long before flotation is even approached. Floods from Summit Lake, British Columbia, start when water level is at least 55 m below that necessary for flotation (Mathews, 1965, 1973), Fisher (1973) has considered this problem and concludes that a well–connected water system normally exists within Salmon Glacier, and that floods from Summit Lake occur when the lake "captures" and then enlarges the normal drainage network. Flotation is probably a sufficient, but not a necessary, condition for triggering subglacial drainage through temperate ice. (Steele Glacier, though cold near its surface, is likely to be near the melting temperature at its bed.)

Instrumentation and field measurements

On 6 July 1978 a pressure transducer was placed in “Hazard Lake” at a depth of 34 m. The transducer used was a solid-state device, mounted in a fluid–filled metal housing, then placed inside a compliant hermetic enclosure to prevent water corrosion. The operating pressure range of the National Semiconductor LX1430AF transducer is 0–300 p.s.i. (0–2 070kPa) corresponding to a range of water depths from 0–210 m. We would have preferred to use an LX1420AF transducer with a smaller depth range and correspondingly higher sensitivity, but the manufacturer could not supply this in time for the field season. As a hermetic enclosure we used an extra–large rubber kitchen glove filled with castor oil. The glove was tied tightly around the transducer lead wires and sealed with epoxy. A collapsible kayak was used to place the sensor in the lake.

Pressure was recorded every 20 min throughout the summer, using a Datel DL2 cassette data logger. The recording station was unmanned, so it was necessary to protect the logger from the environment and grizzly bears. The recorder was recovered in September 1978 and the pressure record scanned for the drainage event.

In July 1979 we returned to the Yukon Territory with an improved lake–level recorder. The modified system used an LX1420AF transducer to measure water pressure and a simultaneous measurement of atmospheric pressure was taken using an LX1702A transducer. Atmospheric pressure fluctuations are sufficiently large to create a pressure variation roughly equivalent to ±0.1 m of water. Unfortunately the lake began draining unexpectedly early, and none of the improved equipment was in place while the lake drained. Instead, we arrived at the lake site during the drainage, and were reduced to monitoring lake level by painting the time of day on rocks as they emerged from the water. A subsequent survey of the rocks gave the water level as a function of time.

After the lake drained we carried out a levelling survey in the lake bed and a vertical photographic survey of the empty basin, and from these prepared an improved bathymetric map of “Hazard Lake” (Fig. 2). Using this map, the hypsometric function A(h _w) (contour area as a function of elevation) was calculated (Table II). The reservoir geometry expressed in this form is required for the simulation model and to convert observations of lake level as a function of time to equivalent water discharge. The total lake volume was found to be 19.62 – 10⁶ m³, about 30% more than that calculated by Collins and Clarke (1977) from an earlier and less accurate map of the lake.

Fig 2 Map of “Hazard Lake” based on 1979 ground survey and aerial photography. Contours measure the depth below high–water mark; the contour interval is 5 m.

Table II. "Hazard lake" hypsometric data 1979

Discharge calculation

From the 1978 water–level record h _w(t) and the hypsometric function A(h _w), discharge is readily calculated using the relation Q(t)= –A(h _w)dh _w/dt. The water-level record and calculated discharge hydrograph for the August 1978 outburst flood are shown in Figure 3. Water–level measurements taken during the July 1979 flood are also included, but the data are too sparse to allow numerical differentiation. Errors in measurement of h _w(t) produce noise in the numerically–evaluated derivative dh _w/dt. This is particularly noticeable at the early stages of the outburst flood when the actual changes in h _w(t) are masked by error noise. The effect on the calculated discharge is even more pronounced since the weighting factor A(h_w) is largest at the onset of drainage. The maximum calculated value of discharge is Q _MAX = 511 m³ s^–1; this is not the true maximum but rather the value at the instant when the water level dropped below the pressure sensor. As already noted, the sensor was placed at a depth of 34 m, and maximum lake depth is 100 m. Because more than 88% of the lake volume lies between the 0–35 m depth contours, most of the outburst flood is recorded despite the fact that the sensor was not placed at the deepest part of the lake. Fitting our data to a power law gives an estimated peak discharge of 641 m³ s^–1. No correction was made for water inflow from Hazard Creek during the flood; the likely magnitude of this inflow is 5 m³ s^–1.

Fig 3 “Hazard Lake” water level record for August 1978 (solid curve descending to the right) and July 1979 (points marked by +) outburst floods. The August 1978 discharge (solid curve ascending to right) is calculated by numerical differentiation of the water–level record; this explains the saw–toothed “noise” component added to the discharge curve. The dashed curves are the predictions of the simulation model. The largest measured value of discharge was 511 m³ s^–1 and the maximum discharge predicted by the simulation model is 547 w³ s^–1. These discharge rates are for the net water flux from the lake; the discharge role for water flowing through the tunnel is obtained by adding the value of water flux into the lake, 5.0 m³ s^–1.

Discharge simulation

Nye (1976) analysed the physics of outburst floods and constructed a theoretical model which he used to calculate the discharge hydrograph for the 1972 Grimsvötn flood. The agreement between theory and measurements is impressive, and suggests the possibility that the theory might eventually be used routinely to predict hydrographs and peak discharge for other ice–dammed reservoirs. With this aim in mind, we have reworked the Nye model so that reservoir geometry is included and the contribution of plastic creep to tunnel closure is retained. Where practical, we have adopted the same notation as Nye (1976).

Basic equations of the model

Nye’s model for jökulhlaups is based on equations for tunnel geometry, continuity, energy conservation, and heat transfer; respectively:

(1)

(2)

(3)

(4)

and the empirical Gauckler–Manning formula relating tunnel cross–section and fluid potential gradient to water discharge

(5)

where

(6)

Expressions similar to the above appear in the earlier work of Mathews (1973), Röthlisberger (1972), Shreve (1972) and Weertman (1972). The derivative dθ _w⁄dt appearing in Equation (3) is the material derivative of θ _w defined as where v is the fluid velocity. In calculating the discharge curve for Grimsvötn, Nye (1976) simplified consideration of heat transfer at the conduit walls by assuming dθ _w⁄dt = 0 in Equation (3). A simulation model that includes this complication as well as a very small correction for the kinetic energy of the water flow is described by Spring (1979). An extraordinarily complete, though computationally challenging, physical model for water flow through ice has been developed by Spring and Hutter (unpublished); the Nye model is shown to be a special case of their more general one.

The assumption that for the 1972 Grimsvötn flood dθ _w⁄dt = 0 in Equation (3) is justified in Nye (1976) on the grounds that the lake temperature is likely to be near 0°C; the measured temperature at the tunnel outlet for the 1954 flood was only 0.05°C, and the calculated hydrograph agrees well with observations. Using the assumption to combine Equations (3) and (4) gives

(7)

implying that water temperature in the conduit adjusts itself so that the heat transfer to the ice occurs at a rate equal to the rate of release of potential energy Q(– ∂ϕ⁄∂s) by the flowing water.

We shall make a similar assumption in our adaptation of Nye’s model, but shall remove the assumption of 0°C lake temperature since, in some cases, thermal energy stored in the lake can make a large contribution to tunnel enlargement. It is helpful to think of θ _w in Equations (3) and (4) as consisting of two parts: θ ₀ due to the temperature at the tunnel inlet, and θ ’, the temperature elevation required to transfer the released potential energy to the ice walls. By definition θ _w =θo + θ′ where θ ₀ and θ′ vary with time and distance along the conduit; at the tunnel inlet θ ₀(0, t) = θ _LAKE. With these assumptions

(8)

and

(9)

We define L′ = L + c _w(θ _w–θ _i) as the “effective latent heat of fusion” for ice and combine Equations (3) and (9) to get the melting rate

(10)

Eliminating Q from Equations (5) and (10) gives

(11)

and substituting this expression into Equation (1) gives

(12)

governing tunnel cross–section as a function of time and distance along the conduit. The first term on the right–hand side is the rate of tunnel enlargement by release of potential energy, the second the rate of enlargement by release of stored thermal energy from the lake water, and the third the rate of closure by plastic creep. By considering only the first term and taking the coefficients of S as constant, Nye was able to integrate Equation (12) directly and obtain simple expressions for tunnel area and discharge as a function of time.

The full complexity of Equation (12) is only grasped when it is appreciated that S, N, ϕ, L′, θ ₀, K ₀, p _i, p, and n can all vary with both time t and distance s along the conduit. The presence of the water–temperature and creep–closure terms forces consideration of what values to assign to θ ₀, p _i, and p. Our point of view, and the simplifying assumptions we make, differ somewhat from those of Nye (1976), so we list them below:

• (1) We postulate that it is the interplay of melting and creep closure in the region of the seal that controls the evolution of the jökuthlaup. Letting the seal be at a distance s ₁ down–stream from the tunnel inlet, the assumption has the effect of replacing the partial derivative on the left–hand side of Equation (12) by an ordinary derivative evaluated at (s ₁, t).

• (2) The tunnel is assumed to have circular cross–section over its entire length and constant Manning roughness n′. Thus, N is constant and equal to, N′ = (4π)^2⁄3 ρ _w gn′ ².

• (3) Following Nye (1976), we replace –∂θ⁄∂s by its spatial average 〈– ∂θ⁄∂s〉 = ρ _w gz _w/l ₀; the lake elevation z _w(t) is allowed to change with time.

• (4) The creep closure of the tunnel is assumed to be governed by Glen’s flow law έ_xy = Bσⁿ _xy , and the tunnel–closure calculation of Nye (1953) is used. Thus in equation (12) K ₀ = 2 B ^(n+1)/2⁄nⁿ . We use an “average” flow law with n = 3, B ₀ = 8.75 – 10^–13 Pa ^–3 s^–1 E = 60.7 kJ mol^–1, θ _i = 273.16 K, and R = 8314 J mol ^–1 deg^–1, giving B = B ₀ exp(–E/Rθ _i) = 2.16 – 10– ²⁴ Pa^₋3 s^–l.

• (5) Ice overburden pressure is evaluated at the seal; thus p _i = ρ _i gh _i, where h _i is ice thickness above the seal.

• (6) The seal is assumed to lie close to the tunnel inlet so that s ₁ ≪ l ₀ and θ (s ₁, t) ≈ p _w gz _w(t) + p(s ₁, t), where p(s ₁, t) is water pressure (Fig. 1). This assumption is likely to be true for “Hazard Lake” because ice thickness is thought to be greatest near the dam and decrease down–glacier from it.

• (7) Water pressure at the seal is approximated by the hydrostatic pressure of the lake p(s ₁, t) = p _w gh _w(t) where h _w(t) = z _w(t) – z ₁ and z ₁, is the elevation of the seal. This assumption ignores a correction for the decrease in water pressure due to fluid flow.

• (8) The lake is assumed to be isothermal so that the inlet water temperature θ _LAKE is not a function of time.

• (9) The seal is assumed sufficiently close to the tunnel inlet that θ ₀, the contribution of the inlet temperature to the total water temperature, is equal to the lake temperature; thus θ ₀ = θ _LAKE. We take L′ = L + c_w(θ _LAKE – θ _i) and neglect a small additional term c _w θ′.

Introducing the above assumptions into Equation (12) gives

(13)

for the rate of change of tunnel cross-section near the seal, and, as stated in the assumptions,

A significant departure from Nye’s framework is to introduce reservoir geometry into the analysis. For Grimsvötn this was not essential because the outburst flood stopped before the reservoir drained completely. For “Hazard Lake” and many other ice–dammed reservoirs, drainage ends when the lake basin is completely empty, so it is important to keep track of lake volume. A second reason for including reservoir geometry as part of the simulation model is that the pressure head at the tunnel inlet decreases as drainage progresses. Thus the fluid potential gradient driving water along the tunnel decreases with time, as does the water pressure opposing creep or fracture closure of the tunnel.

Discharge is simply related to the change of lake volume with time

(14)

where (from assumptions (1), (2), and (3))

(15)

with

Q _IN is the rate of inflow to the lake and Q _OUT the rate of water loss from normal outflow and evaporation. Once the jókulhlaup is fully developed. Q _w and Q _OUT are usually small compared to Q. For “Hazard Lake” they could be neglected entirely, but we choose to retain them for completeness. The relationship between Q _w and Q _OUT will depend on the particular reservoir being studied. For “Hazard Lake” Q _OUT is the water escaping by overflow through the marginal spillway (evaporation loss is neglected). While the lake is filling Q _OUT = 0; when the lake is full Q _0UT = Q _IN – Q. As the jükulhlaup develops, the water flux Q following the subglacial path eventually exceeds the water inflow Q _IN and lake level drops below the spillway, overflow ceases, and Q _OUT = 0.

The coupled equations (13) and (14) are the basic equations of our simulation model. It is evident that a relationship of the form h _W = h _w(V) between water level and lake volume is required to link Equations (13) and (14); this can be generated from the data in Table II.

The discharge simulation consists of solving Equations (13) and (14), together with the geometric relation for the reservoir h _w = h _w(V). The initial conditions are that the lake is full V(0)= V ₀, and the outlet tunnel has some small but finite cross-sectional area, S(s ₁, 0) = S′. The “stopping conditions” are that either the tunnel has resealed (S = 0) or the lake has drained (V = 0). Equations (13) and (14) were simultaneously integrated using a fourth–order Runge–Kutta procedure.

Simulation results Jar “Hazard Lake” outburst flood

The approach we take in modelling the August 1978 “Hazard Lake” outburst flood is to regard the Manning roughness coefficient n′ as the only unknown parameter of the model. The actual cross–sectional shape of the tunnel cannot be observed over much of its length, nor can we determine whether our replacement of – ∂ϕ⁄∂s by its spatial average is a good approximation. If these assumptions are unsatisfactory, the value of n′ which gives the best fit to the measured hydrograph will be contaminated by other uncertainties in the model, and may differ substantially from the true Manning roughness of the outlet path. If, for example, the tunnel has non–circular cross section where R _H is the hydraulic radius (cross–section divided by wetted perimeter of the conduit). For circular cross–section S/4πR ² _H = 1, and for all other shapes S⁄4πR ² _H > 1. If the “Hazard Lake” outlet tunnel has non–circular cross–section, the effect would be indistinguishable from that of increasing n′. Uncertainty in – ∂ϕ⁄∂s could increase or decrease the apparent value of n′.

The input parameters for the “Hazard Lake” jökulhlaup model are summarized in Table III. Lake temperature was not measured in 1978, but measurements taken by D. Liverman in July 1979 gave surprisingly warm temperatures, in the range 5.0–6.5°C. The agreement between the 1978 and 1979 water₋level records (Fig. 3) suggests we can safely transpose the 1979 temperature data to the 1978 flood model. As stated, only the Manning roughness coefficient was adjusted in an effort to fit the simulated curve to observations. (To be perfectly precise, the onset time for the drainage was also adjusted so that the simulation curve and the measured discharge reached peak values at roughly the same instant.) Figure 3 shows the agreement between the observed and calculated water–level records and discharge hydrographs. Agreement between measurements and model predictions is reasonably good, and this inspires confidence in the soundness of Nye’s physical model.

Table III. "Hazard Lake" discharge model

We used a value of n′ = 0.105 m^−1⁄3 s in our simulation. A better fit to the data could have been obtained by adjusting both the Manning roughness and lake temperature. Our aim, however, was not to get the best possible fit, but to test whether an acceptable fit could be found (Table III)using simple assumptions. A good fit can also be obtained by taking the lake temperature as 0°C and the Manning roughness as n′ = 0.009 m ^1⁄3 s; this is an extremely low value for Manning roughness and corresponds to a highly polished circular tunnel. The latter model is unacceptable because it ignores the large amount of thermal energy stored in the lake, but it does illustrate that strikingly different models can produce visually acceptable hydrographs.

For the “Hazard Lake” model, the simulation gives Q _MAX = 547 m³ s^–1 as the maximum net discharge from the lake. (Net discharge dV⁄dt, differs slightly from total discharge; the latter is obtained by adding the input water flux Q _IN to the net discharge.) This of course is not a true prediction of the model because n′ was chosen to fit the data. The maximum measured net discharge for the August 1978 flood was 511 m³ s^–1 and the estimated peak discharge (using a power–law fit to the data) was 641 m³ s^–1. Somewhat remarkably, the Clague–Mathews formula predicts a peak discharge of 551 m³ s^–1 an amazingly good prediction considering the formula depends only on reservoir volume and requires no prior knowledge of tunnel geometry and slope, Manning roughness, or the other variables that appear in the Nye model.

Our value of n′ = 0.105 m^−1⁄3 s for apparent Manning roughness is consistent with Nye’s estimate of n′ = 0.12 m ^1⁄3 s for the 1972 Grimsvötn flood. We have also studied the September 1967 flood from Summit Lake and find n′ = 0.12 m^–1⁄3 s for that event. This narrow range for Manning roughness encourages the hope that the Nye model can be calibrated and used to predict flood magnitude.

Estimation of peak discharge

In the Foregoing section we demonstrated how, by matching measured and simulated discharge records, the apparent Manning roughness of the outlet path could be estimated. If instead, peak discharge is to be estimated, the point of view must be reversed: Manning roughness (as well as the other input variables in Table III) must be known, or guessed, and the model used to predict discharge. Uncertainty in the input data probably contributes as much as the idealized nature of the model to uncertainties in model predictions. We take this as justification for further simplifying the model, so that rough predictions of peak discharge can be obtained without having to undertake computer simulation.

Dimensionless formulation

We begin by defining dimensionless variables describing water depth h ^★ = h _w(t)⁄h ₀ and lake volume V ^★ = V(t)⁄V ₀, and add the following new assumptions:

• (1) The terms Q _IN and Q _OUT in Equation (14) are assumed to be negligibly small compared to the flood discharge Q.

• (2) The decrease with time of the fluid potential gradient ∂ϕ⁄∂s is neglected, so that 〈–∂ϕ⁄∂s〉 is taken as constant in Equations (13) and (14).

• (3) Flotation at the seal is assumed to trigger the outburst flood. Thus ρ _i h _i = ρ _w h ₀ at the onset of drainage and the term (1 – ρ _w h _w⁄(t)⁄ρ _i h _i) ⁿ appearing in Equation (13) simplifies to (1 – h ^★) ⁿ .

• (4) The geometric function h _w(V) for the reservoir is taken to have the simple form h _w(t)/h ₀ = (V(t)⁄V ₀) ^M so that in dimensionless variables h ^★ = V ^★M describes reservoir geometry. We have found that h ^★ – V ^★M gives acceptable fits to all lake basins modelled. The geometric parameter M can be found by least–squares fitting to the hypsometric function for the reservoir, or simply by taking M = V ₀⁄h ₀ A(h ₀). (Note that h ₀ is initial lake surface elevation above the seal, not lake depth.) For a reservoir with vertical walls M = 1; for a paraboloid of rotation M = 1⁄2, and for a cone M = ⅓. For bowl–shaped reservoirs M lies in the range ⅓ < M < 1, and for “horn–shaped” reservoirs 0 < M < ⅓. For “Hazard Lake” M is approximately 0.06; small values such as this seem typical of natural reservoirs. Reservoirs resembling an inverted bowl, such as the subglacial “water cupola” shapes envisaged by Björnsson ([1975]), cannot be modelled using a function of the form h ^★ = V ^★M .

Next we define dimensionless variables for tunnel cross–sectional area S ^★ = S⁄S ₀, time l ^★ = t⁄t ₀, and discharge Q ^★ = Q⁄Q ₀. The characteristic values of tunnel area S ₀, time t ₀, and discharge Q ₀ can be chosen arbitrarily, but detailed dimensional analysis of Equations (13) and (14) shows that

(16)

(17)

(18)

are especially convenient definitions. These characteristic values can be shown to have simple interpretations: Q ₀ is the maximum value of discharge if the terms involving water temperature and creep closure are neglected; S ₀ is the tunnel cross sectional area required to allow discharge at the rate Q ₀; t ₀ is the time required to empty a reservoir of volume V ₀ at the rate Q ₀. For “Hazard Lake” substitution of the input data from Table III gives S ₀ = 21.8 m², t ₀ = 115 h and Q ₀=47.6 m³ s^–1.

The final dimensionless Quantities we define are the "tunnel closure parameter"

(19)

and the “lake temperature parameter”

(20)

The number α characterizes the relative importance of creep closure in controlling the rate of enlargement or closure of the tunnel. For α = 0 there is no creep closure of the tunnel. The number β characterizes the relative importance of lake temperature in controlling the rate of enlargement of the tunnel. For β = 0 the lake and ice temperature are identical and there is no stored thermal energy available to melt ice.

From the above definitions and assumptions, the key equations (13) and (14) can be written in a simple dimensionless form

(21)

(22)

The mathematical description has been simplified to one involving three dimensionless numbers: M describing reservoir geometry and α and β describing the importance of creep closure and lake temperature. This allows the influence of creep closure, lake temperature, and reservoir shape on the magnitude of jökulhlaups to be explored in a simple and fairly general way.

Simple analytic solutions to Equations (21) and (22) can be found by setting α = 0 (i.e. neglecting creep closure). As initial conditions we take V ^★(0)= 1 and S ^★(0) = 0 (the lake is full, the tunnel vanishingly small) and find

(23)

(24)

Discharge is related to tunnel cross–section through the Gauckler–Manning formula, which in dimensionless variables is simply Q ^★ = S ^★4/3; thus the predicted form of the hydrograph is

(25)

Maximum discharge occurs at the instant preceding complete drainage of the reservoir (when V ^★(t ^★) = 0), and the time of maximum discharge can be shown to satisfy the transcendental equation

(26)

Figure 4 shows how dimensionless peak discharge varies with β when α is assumed to vanish; note that the influence of lake temperature on discharge becomes important when β > 0.1. For β ≫ 1 the lake-temperature term dominates and maximum discharge approaches

(27)

Fig 4 Effect of lake temperature parameter (S on dimensionless peak discharge – Large values of β correspond to tunnel enlargement by release of stored thermal energy from the lake water rather than the release of potential energy by the flow. β = 0 corresponds to a lake at 0°C in which no stored thermal energy is available to melt ice and enlarge the tunnel. For this graph α = 0 so that tunnel closure is neglected.

For α = β=0, the effects of lake temperature and creep closure are both neglected and a similar analysis gives

(28)

as the shape of the hydrograph. Equation (28) is identical in form to equation (32) in Nye (1976); the singularity is placed at t*= 0 by an appropriate choice of the integration constant. In Equation (28) maximum discharge occurs at t* = – 3 and = 1. The above results can be used for practical flood magnitude prediction by restoring dimensions to obtain Q _MAX = Q ₀

.

For α ≠ 0, creep closure complicates the analysis and Equations (21) and (22) must be solved numerically. We have explored solutions of these equations for a wide range of values of α, β, and M and present some results in Figure 5. As expected, creep closure is found to be less influential for small values of M than for larger values. This is because horn–shaped reservoirs store most of their volume at the highest level, thus the pressure head in the tunnel remains high for most of the drainage and only begins to drop rapidly at the final stages of the flood. In Figure 5a, the influence of lake temperature has been suppressed by setting β = 0 and M has been taken as 0.05, approximately the value for “Hazard Lake”; in Figure 5b, β = 0 and M has been taken as 0.30 corresponding, we believe, to a rather high value for natural reservoirs. From the data of Table III, we note that the calculated value of a for “Hazard Lake” is α = 1.22 indicating what one might have already suspected, that creep closure is relatively insignificant for outburst floods from this reservoir. Note that for M = 0.05, the difference between the curves for α = 0 and (fig 5) α = 10² is barely perceptible, and that only when α > 10³ does creep closure actually terminate the flood before the reservoir has emptied. For M = 0.30, creep closure becomes important at much smaller values of α. A rough calculation gives α ≈ 10³ for Grimsvötn, so creep closure is likely to be important and may well stop the outburst floods as Nye (1976) suggests.

Fig 5 Effect of tunnel closure parameter α and lake temperature parameter β on the character of discharge hydrographs for outburst floods. Large values of the closure parameter (α ≫ 1) correspond to creep closure dominating the effects of tunnel melting; α = 0 corresponds to no creep closure. Large values of lake temperature parameter p correspond to release of stored thermal energy from the lake water as the dominant mechanism for tunnel enlargement; β = 0 corresponds to a lake at 0°C. Discharge and lime are given in dimensionless variables.

(a) For this figure the reservoir geometric parameter is taken as M = 0.05. This value is similar to that for “Hazard Lake” and corresponds to a “horn shaped” reservoir which stores most of its volume at the highest levels. The lake temperature parameter β is set to zero. The curces for α = 10⁵ and larger values show the flood prematurely terminated by creep closure of the tunnel.

(b) For this figure the reservoir geometric parameter is taken as M = 0.30 and the lake temperature parameter is set to zero. The water pressure in the outlet tunnel falls more rapidly than for case (a); thus the outlet tunnel is more susceptible to creep closure and peak discharge for a given value of a is less than that for case (a). The curves for α = 10⁵ and larger values show the flood prematurely terminated by creep closure of the outlet tunnel.

(c) For this figure the reservoir geometric parameter is taken as M = 0.05. as for “Hazard Lake”, and the lake temperature parameter is taken as β = 10, similar to that estimated for the 1978 outburst flood. Note the different scale for this figure. Comparison with (a) shows that elevated lake temperature affects the magnitude and shape of the discharge curve for a given value of creep closure parameter α.

Figure 5c illustrates how the competing influences of lake temperature and creep closure affect the shape of the discharge curve. For this example M = 0.05, as for “Hazard Lake”, β = 10 and α varies. The calculated value of β for the 1978 outburst flood is β = 11.3 (cf. Table III), so the curve for α = 0 most closely corresponds to the predicted hydrograph for this event. The large value of β indicates that it is the release of the thermal rather than the gravitational potential energy of “Hazard Lake” that makes the greater contribution to tunnel enlargement. Note that creep closure has little influence on the hydrograph for α < 10³.

In concluding, we show how results from our simplified dimensionless model can be applied to the practical problem of predicting the magnitude of outburst floods. First one must calculate α, β, and M using known or guessed information about the glacier and lake. Expressions for α and β arc given in Equations (19) and (20). The geometric parameter M can be found by least–squares fitting using hypsometric data as in Table II or by simply taking M= V ₀/h ₀ A(h ₀). (For “Hazard Lake” least squares fitting gives M = 0.055 5 and the simpler approach gives M = 0.057 0 so the discrepancy is not significant.) Next the characteristic time t ₀ and discharge Q ₀ are evaluated from glacier and lake data using Equations (17) and (18). Finally Equations (21) and (22) are integrated using the appropriate values of α, β, and M to find Q* (t*) and Dimensions are then restored using the relations Q = Q ₀ Q*, t = t ₀ t* and Q _MAX = Q ₀

to obtain the predicted hydrograph and peak discharge.

If the calculated value of α is small (say 0 ≤ α ≤ 100 for a reservoir with M < 0.1), a simpler prescription can be followed. Figure 4 shows how dimensionless peak discharge varies with β. For the “Hazard Lake” flood, β = 11.3, and Figure 4 gives the corresponding value of maximum discharge as roughly = 12.6. Multiplication of this value by the calculated value of characteristic discharge Q ₀ = 47.6 m³ s^–1 gives Q _MAX = 600 m³ s⁻¹ as the estimated maximum discharge. This estimate is in acceptable agreement with the results of the full simulation using Equations (13) and (14). (The example is somewhat artificial since the value of n′ from Table III used to calculate β and Q ₀ was originally chosen to give a good fit to the observed hydrograph.)

In the limiting cases of small and large β (with α small), further simplifications are possible. For β≪ 1 the contribution of lake temperature to tunnel enlargement is negligible so that = 1 and Q _MAX = Q ₀.Thus

(29)

is the predicted peak discharge. The data of Table III give Q _MAX = 47.6 m³ s^–1 in Equation (29). That this estimator badly underestimates the peak discharge for the “Hazard Lake” flood is not surprising because the calculated value of β is large for that event; Equation (29) predicts the discharge magnitude if the lake were at 0°C.

For β ≫ 1 (and α small), thermal energy from the lake is the dominant contribution to melting and = (5β⁄3)^4⁄5. Thus = (5β⁄3)^4⁄5 Q ₀ or

(30)

This estimator gives Q _MAX = 497 m³ s^–1 for the “Hazard Lake” outburst flood, a value that compares well with the predictions of the full simulation based on Equations (13) and (14). The peak discharge is underestimated because tunnel enlargement by release of gravitational energy is ignored.

One aim of this paper has been to put new tools, albeit crude ones, in the hands of engineers and hydrologists who need to evaluate flood potential. Expressions (29) and (30) are intended to supplement the empirical Clague-Mathews formula

(31)

where Q _MAX is measured in m³ s ^–1 and v ₀ in m³. It is interesting that the exponent of the volume term varies from 2/3 to 4/3 in the various prediction formulas. The relationship, if any, between the Clague–Mathews formula and Nye’s model remains mysterious.