Symbols
Introduction
The 1965–66 surge of Steele Glacier, Yukon Territory, Canada, displaced ice as much as 9.5 km (Reference StanleyStanley, 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 Reference Collins and ClarkeCollins 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 (Reference Collins and ClarkeCollins 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 Reference NyeNye (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 (Reference NyeNye, 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 (Reference Clague and MathewsClague and Mathews, 1973)
where Q MAX is measured in m3 s–1 and V 0 in m3. 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 × 106 m3 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 Reference ThorarinssonThorarinsson (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). Reference BjörnssonBjörnsson ([1975]) and Reference NyeNye (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 Reference Mathews(Mathews, 1965,Reference Mathews 1973), Reference FisherFisher (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 – 106 m3, about 30% more than that calculated by Reference Collins and ClarkeCollins 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(hw) is largest at the onset of drainage. The maximum calculated value of discharge is Q MAX = 511 m3 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 m3 s–1. No correction was made for water inflow from Hazard Creek during the flood; the likely magnitude of this inflow is 5 m3 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 m3 s–1 and the maximum discharge predicted by the simulation model is 547 w3 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 m3 s–1.
Discharge simulation
Reference NyeNye (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 Reference NyeNye (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:
and the empirical Gauckler–Manning formula relating tunnel cross–section and fluid potential gradient to water discharge
where
Expressions similar to the above appear in the earlier work of Reference MathewsMathews (1973), Reference RöthlisbergerRöthlisberger (1972), Reference ShreveShreve (1972) and Reference WeertmanWeertman (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, Reference NyeNye (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 Reference Spring and HutterSpring (1979). An extraordinarily complete, though computationally challenging, physical model for water flow through ice has been developed by Reference Spring and HutterSpring 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 Reference NyeNye (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
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: θ 0 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 θ 0 and θ′ vary with time and distance along the conduit; at the tunnel inlet θ 0(0, t) = θ LAKE. With these assumptions
and
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
Eliminating Q from Equations (5) and (10) gives
and substituting this expression into Equation (1) gives
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′, θ 0, K 0, 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 θ 0, p i, and p. Our point of view, and the simplifying assumptions we make, differ somewhat from those of Reference NyeNye (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 1 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 1, 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′ 2.
-
(3) Following Reference NyeNye (1976), we replace –∂θ⁄∂s by its spatial average 〈– ∂θ⁄∂s〉 = ρ w gz w/l 0; 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σn xy , and the tunnel–closure calculation of Reference NyeNye (1953) is used. Thus in equation (12) K 0 = 2 B (n+1)/2⁄nn . We use an “average” flow law with n = 3, B 0 = 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 0 exp(–E/Rθ i) = 2.16 – 10– 24 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 1 ≪ l 0 and θ (s 1, t) ≈ p w gz w(t) + p(s 1, t), where p(s 1, 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 1, t) = p w gh w(t) where h w(t) = z w(t) – z 1 and z 1, 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 θ 0, the contribution of the inlet temperature to the total water temperature, is equal to the lake temperature; thus θ 0 = θ LAKE. We take L′ = L + cw(θ LAKE – θ i) and neglect a small additional term c w θ′.
Introducing the above assumptions into Equation (12) gives
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
where (from assumptions (1), (2), and (3))
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 0, and the outlet tunnel has some small but finite cross-sectional area, S(s 1, 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
2
H = 1, and for all other shapes S⁄4πR
2
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 m3 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 m3 s–1 and the estimated peak discharge (using a power–law fit to the data) was 641 m3 s–1. Somewhat remarkably, the Clague–Mathews formula predicts a peak discharge of 551 m3 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 0 and lake volume V ★ = V(t)⁄V 0, 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 0 at the onset of drainage and the term (1 – ρ w h w⁄(t)⁄ρ i h i) n appearing in Equation (13) simplifies to (1 – h ★) n .
-
(4) The geometric function h w(V) for the reservoir is taken to have the simple form h w(t)/h 0 = (V(t)⁄V 0) 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 0⁄h 0 A(h 0). (Note that h 0 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 Reference BjörnssonBjö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 0, time l ★ = t⁄t 0, and discharge Q ★ = Q⁄Q 0. The characteristic values of tunnel area S 0, time t 0, and discharge Q 0 can be chosen arbitrarily, but detailed dimensional analysis of Equations (13) and (14) shows that
are especially convenient definitions. These characteristic values can be shown to have simple interpretations: Q 0 is the maximum value of discharge if the terms involving water temperature and creep closure are neglected; S 0 is the tunnel cross sectional area required to allow discharge at the rate Q 0; t 0 is the time required to empty a reservoir of volume V 0 at the rate Q 0. For “Hazard Lake” substitution of the input data from Table III gives S 0 = 21.8 m2, t 0 = 115 h and Q 0=47.6 m3 s–1.
The final dimensionless Quantities we define are the "tunnel closure parameter"
and the “lake temperature parameter”
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
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
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
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
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
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
as the shape of the hydrograph. Equation (28) is identical in form to equation (32) in Reference NyeNye (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
0
.
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) α = 102 is barely perceptible, and that only when α > 103 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 α ≈ 103 for Grimsvötn, so creep closure is likely to be important and may well stop the outburst floods as Reference NyeNye (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 α = 105 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 α = 105 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 α < 103.
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
0/h
0
A(h
0). (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
0 and discharge Q
0 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
0
Q*, t = t
0
t* and Q
MAX = Q
0
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
0 = 47.6 m3 s–1 gives Q
MAX = 600 m3 s−1 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
0 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
0.Thus
is the predicted peak discharge. The data of Table III give Q MAX = 47.6 m3 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
0 or
This estimator gives Q MAX = 497 m3 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
where Q MAX is measured in m3 s –1 and v 0 in m3. 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.
Acknowledgements
I owe a great debt to an anonymous referee whose Questions concerning lake temperature led to some major rethinking.
I thank S. G. Collins, D. C. Gayton, R. W. May, and D. E. Thompson for their help with field work; Collins also supervised the lake survey and prepared the base map for this paper. W. A. Wood provided unpublished information and photographs. D. Liverman informed me of the onset of the 1979 drainage, let me use his lake temperature data, and rescued some valuable equipment. B. B. Narod designed the controller for the pressure recorder. The data cassette was transcribed with the help of R. Wilson of the Department of the Environment, British Columbia, and U. Anderson of the Royal Canadian Mounted Police Software support Section. 1 had useful discussions with M. A. Church, D. Jones, W. H, Mathews, M. C. Quick, and G. J. Young.
I am grateful to Environment Canada, the National Research Council of Canada, the University of British Columbia Committee on Arctic and Alpine Research, and the Department of Indian and Northern Affairs for financial support and to the Arctic Institute of North America, University of Calgary, for logistic assistance. Parks Canada gave permission for this research to be done in Kluane National Park and their interest and co-operation is appreciated.
