, where 1. Introduction
The release of water from glaciers in catastrophic floods, known as jökulhlaups, has a long history of study prompted by their environmental consequences and potential threat to human activity. The challenge clearly is to predict their timing and magnitude. But we are still short of this goal even though a physical understanding of flood mechanics is emerging.
In this paper, we discuss a long-standing issue in the problem of flood prediction. By analyzing flood data from 10 marginal ice-dammed lakes, Reference Clague and MathewsClague and Mathews (1973) identified an empirical relation between the peak discharge, Q max, and the water volume drained during the flood, Vt , in the form
The constants are K = 75 and b = 0.67 (with Q max measured in m3 s−1 and Vt measured in 106 m3) and the statistical fit was found to be remarkably good. A careful and extended study by Reference Walder and CostaWalder and Costa (1996) revealed greater scatter when obtaining a fit to 26 lakes, but the exponent, b ≈ 2/3, seems robust for lake systems that produce jökulhlaups subglacially (as opposed to those that are drained by water flowing along the ice margin, bounded by valley walls). It is thus natural to suppose that the power law in Equation (1), the so-called “Clague–Mathews relation”, has an origin in the physics governing flood evolution, and this is what we set out to examine here.
The starting point of our enquiry is the modern thermomechanical description of jökulhlaups. Its constituents have been investigated by several authors, notably Reference RöthlisbergerRöthlisberger (1972), Reference MathewsMathews (1973) and Reference Glazyrin and SokolovGlazyrin and Sokolov (1976), but a complete theory first appeared in the celebrated paper by Reference NyeNye (1976). In this description the ice-dammed lake drains via a subglacial water channel linking it to the snout. During the rising flood stage, the channel expands by melting of its ice walls due to heat dissipation by the water being conveyed. A positive feedback between the rate of melting and channel enlargement leads to an escalating discharge (the rapid growth observed in the flood hydrograph), but this does not continue indefinitely, as later, at lower lake levels and reduced water pressure, ice deformation closes down the channel and the flood subsides. Reference NyeNye (1976) formulated the governing equations to model this behaviour. Although this model is not without shortcomings (Reference Fowler and NgFowler and Ng, 1996), it has been used to simulate hydrographic data, with a fair degree of success, from jökulhlaups from both marginal and subglacial lakes.
Explaining the trend behind Equation (1) is motivated by the potential for predicting flood magnitude given a known lake volume. Moreover, the trend persists over a vast range in Q max and Vt , and poses a crucial test for our current understanding of flood mechanics. Previous work has focused on simulating discharge variations over time, especially the flood hydrograph (e.g. Reference Spring and HutterSpring and Hutter, 1981; Reference ClarkeClarke, 1982; Reference BjörnssonBjörnsson, 1992). Instead, we use the Nye model here to illuminate the nature of existing Q max−Vt data, and show that insight can be gained into the mechanics of the floods through their collective behaviour.
2. Background
2.1. Discharge–volume data for jökulhlaups
Pertinent data compilations have been presented by other investigators. Reference Walder and CostaWalder and Costa (1996, fig. 2 and table 1) studied a dataset of “subaerial” ice-dammed lakes in a variety of settings (e.g. at the glacier margin or at the confluence of glacier tributaries), including a revised version of the original data of Reference Clague and MathewsClague and Mathews (1973, fig. 2 and table 1) as subset. Reference BjörnssonBjörnsson (1992, fig.16) worked on a dataset based on jökulhlaups originating from subglacial lake Grímsvötn in Iceland, one of the most extensively studied systems in the world. These two datasets are probably the most comprehensive to date and are shown in Figure 1. Only floods that are thought to have drained subglacially are included. The power laws derived by earlier authors are straight lines in the figure, and the corresponding K and b values are listed in Table 1. Full reference to the data sources is provided by Reference Walder and CostaWalder and Costa (1996).
Fig. 1 Peak flood discharge plotted against volume drained for subaerial lakes (crosses) (taken from the dataset of Reference Walder and CostaWalder and Costa, 1996) and for Grímsvötn (circles) (from Reference BjörnssonBjörnsson, 1992; Reference Gudmundsson, Björnsson and PálssonGudmundsson and others, 1995). Envelopes group together multiple events from the same lake if there are three or more. Numbering of regression power laws followsTable 1.
Table 1. Constants in the regression power laws of Clague–Mathews form derived by earlier authors for flood discharge–volume data. r 2 is the correlation coefficient. Power law 3 is based on the Grimsvotn dataset only
Although Figure 1 displays a coherent trend, several complications must be taken into account when interpreting it. First and foremost, the scatter of the data is significant. Where data on multiple flood events are available for a given lake, they commonly occur in clusters spanning up to an order-of-magnitude range in Q max and Vt . If one were to construct a Clague–Mathews relation for each lake, the resulting best-fit slope (the value of b) would generally differ significantly from the overall slope derived for all lake systems, which is about 2/3.
Another feature is that many of the outbursts are known to have emptied the lake (“complete draining”), while others have terminated with the lake still partially filled (“incomplete draining”). (Reference Clague and MathewsClague and Mathews (1973) referred to Vt as the water volume stored in the lake, or the lake capacity, which is appropriate because draining was thought to be complete in their examples. In this paper, to avoid ambiguity, we define Vt to be the volume drained.) The Grímsvötn jökulhlaups belong consistently to the second category. A Clague–Mathews fit for them also seems to be exceptionally good compared to other systems (Fig. 1), despite a drastically different value of the exponent; Table 1 shows that b calculated by Reference BjörnssonBjörnsson (1992) for Grímsvötn is almost three times the other b values. These observations may not be coincidental, and we discuss them later.
A plausible strategy for us to follow is to study how a physics-based relation between Q max and Vt can arise in a model. This notion is not new. Equation (1) is suggestive of a scaling law in the description of the floods, for instance, and this aspect has been investigated by Reference ClarkeClarke (1982). But the complications noted above reveal important subtleties. The general trend pertaining to b ≈ 2/3 seems to be one that exists between the different data clusters, not within individual clusters, and so the factors governing the discharge–volume characteristics of a number of lakes and those of a single lake, while they may be closely related, are probably different. In addition, the possibility of the lake running out of water creates a (discontinuous) threshold in the flood system. Because of this, the data recorded in Figure 1 may not be carrying the same information about the floods’ evolutionary mechanics, in the sense that some of them have been influenced by this threshold while others have not.
2.2. The Nye model
Figure 2 shows a hypothetical jökulhlaup system, consisting of a source lake drained by a subglacial channel with crosssectional area S(x; t), where x denotes downstream distance from the lake and t denotes time. Reference NyeNye’s (1976) theory describes the evolution of this channel as a result of ice melting and ice deformation (creep closure) along its length, using the equation
where m is the melt rate of ice (mass per unit length per time), ρ i is ice density, and the last term of this equation relates the closure rate of the channel to the difference between the overburden ice pressure, p i, and the channel water pressure, p w. It is convenient to define an effective pressure, N = p i − p w, for subsequent calculations. n (≈ 3) is the exponent in Glen’s flow law. The classical theory assumes a cylindrical channel, and the constant K 0 has been derived by solving the corresponding creep closure problem for the ice (Reference NyeNye, 1953).
Fig. 2 Schematic diagram of a jökulhlaup system and the definitions used in our mathematical model. Inset (top right) illustrates the depth correction arising from a floating ice shelf, described in the text.
The key idea in this model is that a flood arises as an imbalance between melting and closure when the channel is coupled to the lake, the draining of which controls N (through water pressure). m is deduced from hydraulic considerations. If Q(x; t) is the water flux in the channel, then the equations for mass, momentum and energy conservation, respectively, are
where ρ w is the density of water and L is latent heat. (An equivalent version of Nye’s equations written in terms of the effective pressure N has been given by Reference FowlerFowler (1999).) Equation (4) is a Manning parameterization for turbulent flow, in which the constant F 1 takes into account the wall friction of the channel and its cross-sectional geometry. Ψ + ∂N/∂x is the hydraulic gradient driving the flow, and Ψ is that component due to topography alone (a function of x only), defined by
αb is the slope of the glacier bed, and g (= 9.8 m s−2) is gravitational acceleration.
We need not be concerned with the full model here. Indeed, Equation (5), which relates the rate of heat dissipation to the melt rate, is already simplified, as it neglects the heat that raises the water temperature (slightly) above the melting point. Reference NyeNye (1976) and Reference ClarkeClarke (1982) considered other simplifications, the most useful of which allows the model’s x dependence to be removed. By a scaling analysis, Reference NgNg (1998) has shown that such an approximation is valid when pressure conditions at the snout do not influence the lake-end of the system strongly.Footnote 1 Under these conditions, N refers to the effective pressure at the seal (close to the lake), and Equations (3–5) become
On eliminating m and S from Equations (2), (8) and (9), we find
which is a compact description of channel evolution arising from melting and closure (the first and second term on the right, respectively). In our applications, we use a spatially averaged value of Ψ calculated from the ice thickness at the seal and the bed elevation drop from there to the glacier snout. We also assume F 1 = 530 kg m−8/3 based on a Manning roughness coefficient 0.1 m−1/3 s found by Reference NyeNye (1976) and Reference ClarkeClarke (1982), and the constants ρ i = 900 kg m−3, ρ w = 1000 kg m−3, L = 3.34×105 J kg−1.
Lake dynamics
If h i is the ice thickness at the seal and h w is the lake surface elevation relative to the lowest point of the lake (Fig. 2), then the effective pressure is
and we may describe the time variation of N through surface lowering of the lake as its volume V is reduced. This formulation is similar to the one introduced by Reference ClarkeClarke (1982), but not identical since we allow a vertical offset, h
s, between the lake bottom and the seal. Given the lake hypsometry A = A(h
w), where A is lake area, we have 
. The mass-continuity equation for the lake is simply dV/dt = −Q + Q
IN. An alternative way of writing this is
This equation couples with Equation (10). Following Reference ClarkeClarke (1982), we describe the lake hypsometry using a power-law function and let
in which h 0 is the flotation lake level (at which point N = 0), given by
A 0 is the corresponding lake area (Fig. 2).
When Equation (12) is applied, a correction is necessary if a floating ice shelf covers the lake, as encountered at Grímsvötn. Then it is the basal area of the ice shelf that determines the lake volume change for a known adjustment of water level. However, because of the way terms representing water pressure have been defined in Equation (11), h w in this case refers to an equivalent lake water depth taking into account the ice shelf, and we evaluate A in Equation (12) at the elevation h w − ρ i d s = ρ w, where d s is the average ice-shelf thickness, and not at the elevation h w (see inset in Fig. 2).
Non-dimensionalization
To facilitate analysis of Equations (10–13), we let [t] be the time-scale and assign scales for the other variables [h w], [Q], [N], [A], [V], writing t ⋆ = t/[t], h ⋆ = h w/[h w], and Q ⋆ = Q/[Q], etc. (The stars label dimensionless variables.) If we use the natural scales [h w] = h 0, [A] = A 0, [V] = h 0 A 0/(β + 1) and balance the scaled size of all the terms in Equation (10), and also the left of Equation (12) with [Q] arising from its righthand side, then the dimensionless model is (dropping the stars⋆)
and these are supplemented by the equations
The definitions of [Q], [N] and [t] are given in the Appendix. We have introduced the new parameters
as dimensionless measures of the lake depth, meltwater input and ice-shelf thickness, respectively. Equations (15–18) provide a framework in which jökulhlaups from different lake systems can be analyzed together once their topographic settings have been defined via the constants h 0 (flotation water depth), A 0, β (hypsometry), Q IN (influx to the lake) and Ψ (topographical part of the hydraulic gradient). In this dimensionless model, γ controls the size of N as the lake drains, and thus it plays a part in determining the relative importance of channel closure and melting in Equation (15). The water level varies in the range 1 ≥ h ≥ ξ, so the dimensionless (maximum) lake depth is given by γ(1 − ξ). We have ξ = 0 for a subaerial lake (no ice shelf), and typically ν ≪ 1.
3. Theoretical Development
3.1. Fundamentals
Properties of the model may be studied on the Q−N phase plane (Fig. 3a), which shows the solution trajectory as time proceeds. This construction is instrumental for understanding the Clague–Mathews data of section 2.1. It is useful to consider the state of the flood system with reference to three special curves on the phase plane, known as nullclines (a mathematical term from the study of dynamical systems).
-
(i) Q = N 4n , on which the rates of melting and closure balance, and the time derivative dQ = dt = 0. This is a dimensionless version of Reference RöthlisbergerRöthlisberger’s (1972) condition for steady channel flow.
-
(ii) Q = ν (≪ 1), on which lake inflow and outflow equalize, and dN/dt = 0.
-
(iii) N = γ(1 − ξ), on which the lake is empty, and dN/dt “blows up” (except where it meets nullcline (ii)).
Fig. 3 (a) A sketch showing the phase-plane solutions of the coupled model of Equations (15–17) and nullclines (i–iii) described in the text. Arrows point in the direction of time. (b) Definition of Q max , N1 and N2 for a flood trajectory. (c) Peak-discharge–effective-pressure functions N1(Q max ) and N2(Q max ) derived from (a) and (b).
These nullclines may also be deduced from the non-linear ordinary differential equation
which follows from dividing Equation (15) by Equation (16) and substituting for A with Equation (17). The phase trajectories have zero slope where they meet nullcline (i) or (iii), and infinite slope where they meet nullcline (ii). Because the formulation of creep closure breaks down for N < 0, the model cannot describe what happens when the lake is above the flotation level. We consider solutions in the range 0 ≤ N ≤ γ(1 − ξ) and Q ≥ 0.
We highlight two important problems when the model is used for flood simulation. Given suitable initial values for Q (16) and N that resemble conditions at the beginning of a jökulhlaup (typically, low Q and low N), Equations (15–17) may be solved, giving a clockwise arc on the phase plane, with the rise and fall of Q depicting growth and recession of a flood. An example is shown in Figure 3a (thick curve). The (17) first problem is that the physical basis of the Nye model fails to describe the shut-down of channelized subglacial drainage at flood termination (we suppose this happens in reality). The solution therefore becomes invalid after a flood, even though, mathematically, it provides a logical (but bogus) prediction of a continued, spiral trajectory on the phase plane (Fig. 3a), mimicking a succession of flood events separated by stages of slow refilling. This shortcoming is confirmed by observations (e.g. at Grímsvötn), where neither the magnitude nor the triggering lake level of successive floods increases indefinitely.
With the model’s applicability restricted to one flood cycle, there is a second problem, relating to prediction. The initial conditions for Q and N determine the precise trajectory (out of an infinite number of possible ones on the phase plane) followed by a flood. However, it is not clear how to specify them because the model lacks a mechanism of flood initiation, so the peak discharge, Q max, cannot be predicted. Furthermore, although the model may be tuned to match an observed hydrograph, the result is extremely sensitive to the chosen initial conditions. To investigate how the dam-break mechanism might determine these conditions, Reference FowlerFowler (1999) has modelled the spatial development of subglacial drainage at the seal. His theory indicates that during flood initiation the lake level begins to drop at a critical N-value (N crit) that depends on the basal melting rate (at the seal) and also the refilling rate ν. He thereby proposed an explanation for the regular lake levels (below flotation) at which the normal Grímsvötn jökulhlaups tend to be triggered, corresponding to N crit ≈ 6 bar (Reference BjörnssonBjörnsson, 1988, Reference Björnsson1992; Reference Gudmundsson, Björnsson and PálssonGudmundsson and others, 1995). Tentatively, this would suggest taking (N, Q) = (N crit, ν) as initial conditions on the phase plane. But while Fowler’s theory is insightful and may lay the path for flood prediction, there are serious obstacles in applying (indeed, testing) it, owing to the difficulty of measuring either the basal melting rate or ν accurately in the field.
3.2. Peak-discharge–effective-pressure relation
Despite uncertainty as to which trajectory is selected, we may determine from the model how the change in effective pressure in each (and every possible) flood is related to its peak discharge. Figure 3 illustrates this idea. Each trajectory has associated with it a value of Q max where it crosses nullcline (i), but also unique values of the effective pressure: N = N 2 at the end of the flood and N = N 1 at the beginning (see Fig. 3b, where we define flood initiation and termination by the interception points of the trajectory with nullcline (ii)). Moreoever, the relation between the three quantities may be studied by plotting Q max vs. N 1 and Q max vs N 2 for different trajectories (Fig. 3c).
We think of N 1 and N 2 as functions of Q max here, even though the implied dependences are implicit (and do not involve time). This treatment may not be intuitive but it provides a conceptual basis for understanding discharge–volume data. During a flood, the change from N 1 to N 2 corresponds to a certain lake volume reduction. Therefore, the functions in Figure 3c can be used to calculate the volume drained, Vt , directly, knowing Q max. Notably, Figure 3c indicates that a more destructive jökulhlaup (with a large value of Q max) begins at a lower N value (higher lake level) and terminates at a higher N value (lower lake level) than a less destructive jökulhlaup (with small Q max) would do. Thus Vt increases with Q max for a given lake system.
Crucial to the preceding argument, Figure 3 succeeds in demonstrating that in the Nye model: (a) a flood culminates when the rates of melting and closure balance, so Q max determines the state of the system completely — the value of N then, and hence the entire history of the flood, especially the values of N 1 and N 2; and (b) different phase trajectories do not cross, in such a way that N 2 increases (and N 1 decreases) with Q max. It is much more difficult to arrive at the dependence of Vt on Q max (or vice versa) by an alternative argument based directly on Equations (15) and (16), where rates of change with time are involved.
Notice we impose lower and upper cut-offs to N 1 and N 2, respectively, at N = 0 and N = γ(1 −ξ), for those trajectories high on the phase plane (with sufficiently large Q max) that do not meet nullcline (ii) (see Fig. 3a). In this connection, we define a critical value Q crit in Figure 3c, where, if Q max ≥ Q crit, the corresponding flood empties the lake completely (its trajectory meets nullcline (iii), and N 2 = γ(1 − ξ)), whereas below this critical value (Q max < Q crit) the lake is only partially drained by the flood. This distinction between the two “drainage styles” (complete vs incomplete draining or emptying) becomes relevant in section 4, where we analyze flood data on both the Q max– N and Q max– Vt domains.
3.3. Results for a single lake
Analytical solution of Equation (19) is not straightforward, so we evaluate N
1(Q
max) and N
2(Q
max) numerically. We treat Equations (15) and (16) as parametric equations and integrate them forward and backward in time from different points on nullcline (i), where 
, with Q
max
> ν, to obtain N
1 and N
2 subject to the cut-offs in effective pressure. In back integrations N
1 (≥ 0) is given by the crossing with nullcline (ii), at the time of flood initiation hypothesized in our model. Typically the observed value, N
1o, differs from this owing to the actual initiation process. N
1o may be calculated from the initial lake level, which has been measured for some flood events.
Knowing the change from N 1 to N 2 in a flood, the lake volume drained, derived from Equations (17) 2 and (17) 3, is
Since N 1 and N 2 are functions of Q max, a model-based Clague–Mathews relation for floods originating from a given lake can now be established. In dimensionless form, we write
or
An example of computed results is shown in Figure 4, where we plot Q max against N 1 and N 2 (solid lines), taking n = 3 and the parameters ξ = 0, β = 1, γ = 3 and ν = 0.001. The corresponding function f calculated from Equation (22) is also indicated (dashed line).
Fig. 4 Double logarithmic plot of the discharge–pressure functions N1(Q max ) and N2(Q max ) (solid lines) and the discharge–volume function Q max = f(V t ) (dashed line), computed from the dimensionless model with n = 3, ξ = 0, β = 1, γ = 3 and ν = 0.001.
3.4. Comparing results for different lakes
The fact that f is an increasing function appears consistent with the discharge–volume data. However, comparing f directly with the overall trend in Figure 1 would be highly misleading, because Equation (21) (or (22)) is dimensionless and its application requires the Clague–Mathews data to be rescaled first. Since the scales [Q], [N], [V] depend on the values of A 0, h 0, Ψ and β , the data for different lakes are rendered dimensionless in different scales; the effect of this is to translate the data clusters on the log–log plot. When the size and geometry of the systems are different, it clearly is not possible to recover f by fitting a power law to the dimensional data.
Alternatively, consider the reverse of putting Equation (21) associated with two different flood systems, 
and 
, respectively, into dimensional form. (We reintroduce the⋆ here to highlight the dimensionless variables.) Then Q
max/[Q]
i
= fi
(Vt
/[V]
i
) (i = 1 or 2), and by defining the functions g
1 and g
2 via 10
gi
(
x
) = fi
(10
x
), we can write
These equations show that even for two similar systems — similar in the sense that they share identical β , γ, ν and ξ, such that g 1 ≡ g 2 — the corresponding Clague–Mathews relations would appear displaced from each other in Figure 1 by the amounts log10([V]2/[V]1) horizontally and log10 ([Q]2 = [Q]1) vertically. Generally, g 1 ≠ g 2 means that the relations will differ in shape also.
Our theory reveals other factors influencing the behaviour of each lake system in Figure 1 that complicate its interpretation. For multiple, incomplete-draining flood events, Q max−Vt data for f within a limited range may be approximated by a power law with a well-defined Clague–Mathews exponent b. (For instance, b ≈1–2 in Figure 4.) In practice though, this property will be destroyed if: (a) the flood-initiation lake levels differ markedly from those corresponding to N 1 defined in our model (then Vt departs from that predicted by theory, “blurred” by the actual values of N 1o), (b) the flood system geometry has changed significantly between events, or (c) some of the floods involve complete draining of the lake. The last factor is particularly problematic because the corresponding data will not reflect the curved portion of f at all. Armed with the prerequisites, we are ready to examine existing data more closely.
4. Model–Data Comparison
4.1. Grímsvötn
Dimensionless Qmax–N domain
We begin by comparing N 1(Q max) and N 2(Q max) derived from theory with N 1o and N 2o calculated from the lake levels. (Subscripto denotes the observed values.) The Grímsvötn data are well suited for this purpose because the initial and final values of h w have been documented. Moreover, draining of the lake is incomplete (Reference BjörnssonBjörnsson, 1988), which makes possible a comparison between N 2 and N 2o below their upper cut-off.
We use the constants listed in Table 2 and analyze the data for 10 jökulhlaups of the last century compiled by Reference Gudmundsson, Björnsson and PálssonGudmundsson and others (1995), reproduced in Table 3. (The volcanogenic floods of 1934, 1938 and 1996 are omitted because they probably involved a high lake-water temperature several °C above the melting point, not described by our simplified Nye model.) As can be seen from the definitions in the Appendix, the scales for non-dimensionalizing the data depend on the values of K 0 and n specified. Good agreement between model and data is achieved (in Fig. 5, discussed below) with K 0 = 5 × 10−24 Pa−3 s−1 when n = 3. Then [Q] = 4.27 × 105 m3 s−1 and [N] = 15.3 bar, and the model parameters are
The corresponding Q max, N 1o and N 2o values (dimensionless) appear in Table 3. Figure 5a illustrates these discharge–pressure data, which are divided into two groups by nullcline (i) (Q = N 12 in this case).
Fig. 5 (a) Results of matching the model-derived functions N1(Q max ) and N2(Q max ) with the observed values N1o and N2o for 10 floods from Grímsvötn. Model constants include n = 3, K0 = 5 × 10−24 Pa−3 s−1 and the ones listed in Tables 2 and 3. Circles refer to the N2 values calculated with the ice-shelf thickness specific to each flood. (b) Results of matching the model-derived Clague–Mathews function Qmax = f(V t ) with discharge–volume data. Model parameters are the same as in (a). All figure axes are dimensionless.
Table 2 Model constants for the Grímsvötn system
Next we calculate N 1 and N 2 from the theory, assuming the parameters in Equation (24). Because the ice shelf at Grímsvötn has been thickening, we are dealing with a system whose geometry has changed between the floods, and this is reflected in the different values of ξ (Table 3). Our first comparison between data and theory assumes the mean value ξ = 0.403 to define a fixed system. The (continuous) functions N 1 and N 2 have been computed and plotted over the data in Figure 5a. N 2 agrees quite well with the observed N 2o. The discrepancy between N 1 and N 1o is expected because the former is based on a hypothetical definition of flood initiation, but it is relevant to the volume prediction below (Fig. 5b). In the second comparison, we correct for the variability of ice-shelf thickness and calculate N 2(Q max) for each flood by assuming the value of ξ for that year. This gives a much better agreement than before (circles in Fig. 5a), because the larger (earlier) floods involved an ice shelf that was generally thinner and had a larger basal area, leading to lake-level drops that were less significant.
Table 3 Characteristics of 10 recent jökulhlaups from Grímsvötn. Values in the last four columns are dimensionless and based on the scales defined in section 4.1
Dimensionless Qmax–Vt domain
Using Equation (22), we proceed to convert N 1o (replacing N 1) and N 2 into Vt , the volume drained. The model results are shown in Figure 5b, together with data derived from the measured Vt in Table 3. We give the theoretical curve (f) based on the mean value of the ice-shelf thickness and the mean value of N 1o, and also isolated predictions (circles) based on the year-specific values of these parameters. Knowing the value of N 1o for each flood is crucial to the agreement, since a small change in the initial lake level introduces large errors in Vt . The predicted volumes may be seen to be as much as 50% in error with the observed, but this is probably similar to the uncertainty in the field measurements (for both Q max and Vt ). Since the data span an order-of-magnitude range, the agreement is really rather impressive.
We have thus demonstrated that the discharge–volume characteristics of the Grímsvötn jökulhlaups (in dimensional form (Fig. 1)) can be explained satisfactorily by the model. Given the agreement in both stages of our comparison (Fig. 5a and b) and the well-constrained geometry and data from Grímsvötn, this result is convincing. While the operation of subglacial water drainage during a flood has not been observed directly, the idea of competing conduit “expansion” and “constriction” seems capable of capturing its dynamics, lending credibility to the mechanisms envisaged by Reference NyeNye (1976).
4.2. Other lake systems
The previous analysis needs to be applied to other systems, for which, similarly, knowledge of their size and geometry is required to scale the observations. Unfortunately, a search through the literature shows that usable data are scarce, as few studies provide reliable estimates of the ice-dam and lake geometry (h i, h s, β), not to mention lake levels. This inevitably limits what more can be deduced, but the potential value of the model is not exhausted.
Implications of complete-draining flood events
Complete emptying of the lake seems to be a common occurrence (e.g. Table 4). In such cases, our model shows that the Q max−N data contain little information on flood mechanics because N 2 lie on nullcline (iii), which is a threshold in the system. The drained volume is then equal to the initial lake volume (in dimensionless terms, (1 − ξ − γ −1 N 1o) β +1), so that variations in the flood-triggering lake levels are directly reflected in a scatter of the data in the Vt axis. This effect may play a leading role in shaping some of the data clusters in Figure1. In contrast, the scatter is subdued in the Grímsvötn data, where it appears superimposed on a dominant rising trend of the function f associated with partial draining of the lake. (In cases both of complete and of incomplete draining, variations in the conditions at the time of flood initiation are responsible for different values of Q max in the data (section 3).) For a given lake system, a power-law fit to discharge–volume data involving (any) complete-draining events will therefore misrepresent its dynamics. Data interpretation is further complicated if the drainage style of some of the floods is not known, or if the ice-dam configuration changed drastically between floods, altering the maximum lake capacity (e.g. the right and left clusters for Grænalón in Figure 1, corresponding to the pre- and post-1950 periods, may be due to this).
Table 4 Parameters for 11 ice-dammed lake systems in order of decreasing lake volume
This raises interesting questions: What sets Grímsvötn apart from the other lakes in terms of drainage style? Is this related to its subglacial, rather than marginal, location? Moreover, if Figure 1 is plagued by data referring to complete drainage, and f is irrecoverable, why is there an overall trend at all, and what does it indicate?
The question of lake drainage style
According to Figure 3c, complete draining occurs when (dimensionlessly) Q max exceeds the value Q crit (section 3). For every lake in Figure 1, one could in principle calculate Q crit from the model, and test whether or not the measured value of Q max in each flood is consistent with the observed drainage style. But since Q crit and the discharge scale involved in non-dimensionalizing Q max both depend on the closure rate constant K 0 that is assumed, the results will not be conclusive. Moreover, we have sufficient data only for a handful of lakes (Table 4) for such analysis.
To explore why some lakes tend to drain partially, or completely, or do both, we use a more general approach based on Figure 3. The dimensionless lake depth γ(1 − ξ), which controls the position of nullcline (iii), measures the importance of channel closure in the model, and Q crit is an increasing function of it. Suppose the dam-break conditions are “invariant” — we assume the flood is initiated at a fixed place on the phase plane — then a sufficiently “shallow” lake (γ small, and Q crit small) will empty in the flood, whereas a “deep” lake (γ large, Q crit large) will drain only partially, other model parameters being equal. To look for such pattern, note that γ (= ρ w gh 0/[N], where [N] is defined in the Appendix) is proportional to
and in Figure 6 we plot this ratio (taking n = 3) against the maximum volume drained for the lakes listed in Table 4. There is some evidence for the pattern mentioned above, but this result is not unequivocal. Although γ(1 − ξ) for Grímsvötn and γ for Grænalón (after 1950) lie near the top of the range, they are not the highest values among the lakes. (Grænalón emptied in the floods before 1940 and the corresponding γ is highest.) Neither are they much greater than the γ values for other lakes that drain completely. We have verified that this result holds for 2 ≤ n ≤ 4.
Fig. 6 Plot showing the ratio in Expression (25) (proportional to the dimensionless lake depth, γ) against lake volume for 12 jökulhlaup systems. The Grímsvötn value has been corrected by the multiplicative factor 1 − ξ where ξ = 0.403, due to the presence of an ice shelf. Open/filled circles refer to complete/incomplete lake-emptying. Lake hypsometry exponents β , where known, are shown next to the data points.
The preceding comparison is exact for fixed β , because then the drainage style is controlled by γ alone. However, β does vary between systems (Fig. 6), and several of them have rather high values, for instance, indicating that they involve “horn-shaped” reservoirs rather than “bowl-shaped” ones (using Reference ClarkeClarke’s (1982) terminology) in the side-valleys. The effect of this is to accelerate lake area reduction at low water levels, promoting complete drainage. A small β has the opposite effect, but comparison with the pre-1940 Grænalón value (β ≈ 1.79) in Figure 6 suggests that the Grímsvötn value (β ≈ 1.94) is not small enough to guarantee flood termination before the lake empties. Partial draining of Grímsvötn and post-1950 Grænalón probably reflects a tendency of their flood initiation mechanism to launch trajectories lower on the phase plane than in the case of the other lakes, which is against the supposition of invariant initial conditions. An understanding of flood initiation thus seems to be necessary when comparing data for different lakes.
Multiple systems: rising trend of the Clague–Mathews relation
Our results so far indicate that model and data are compatible. For a single lake, the discharge–volume data should:
Record the function f predicted by theory (Equations (21) and (22)) and follow a coherent trend (with a log–log slope of ∼1 to 2 (e.g. Fig. 4)) if they refer to flood events involving only partial draining of the lake. The Grímsvötn jökulhlaups are the definitive (and only) example (section 4.1).
Show little or no coherent trend when complete draining is involved. This may be responsible for the scatter of data within clusters in Figure 1.
That voluminous floods should proceed with high peak flows might seem intuitive, and happens in the first scenario above. What remains perplexing is the overall rising trend among different lake systems (Fig. 1), given that many of the lakes drain completely. Explaining this has important implications for predicting the maximum possible Q max. As we argue below, part of the answer lies in the issue of “scaling”. The trend reflects different dynamics of the flood systems arising from their geometry and topographical setting (which are conditioned by other processes). But these are not the only factors, and the precise reason for the exponent b being close to 2/3 seems to be inextricably linked to the nature of flood initiation.
Two lake systems are “similar” to each other if, after non-dimensionalization (section 2.2), their corresponding model parameters are identical. Then the phase-plane solutions governing their dynamics are the same, and the natural scales of the systems, [Q] and [V], determine the locations of the (lake-specific) Clague–Mathews relations in Figure 1.Footnote 2 It is not possible to predict where the next flood event will occur along each Clague–Mathews relation, because of uncertainty in the conditions at the time of flood initiation. If, however, the same conditions of flood initiation (initial values of dimensionless N and Q on the phase plane) apply to the lakes, then a Clague–Mathews relation describing the flood data for both lakes would have the slope of the line that passes through (log10 [V]1 , log10[Q]1) and (log10[V]2, log10[Q]2). (Here subscripts label the lakes, and we refer to the line as the scaling line.) In actual observations, departure from the theoretical scaling slope brings out the dissimilarity between the systems, in terms of dynamics, conditions at flood initiation or both.
This idea may be generalized for several lakes. In Figure 7, we show the values of [Q] and [V] derived for the lakes of Table 4 on a logarithmic plot (assuming n = 3 and K
0 = 5 × 10−24 Pa−3 s−1), and construct an approximate scaling line through them in order to make a comparison with actual flood data from Figure 1. [Q] and [V] are scales, so our aim is not to produce a match by tuning K
0, even though increasing it lowers the scaling line in the figure towards the data. In Figure 7 the axes are dimensional, and flood data are recorded at positions (log10
Vt
, log10
Q
max), where Vt
≈ [V] (due to the common occurrence of complete-draining events) and
. The dimensionless peak discharge
depends on the systems’ parameters and the trajectories picked up at flood initiation. The vertical offset between scaling data and flood data indicates that
is not the same for all lakes; typically
.
Fig. 7 Plot showing discharge–volume scales [Q]−[V] (plus signs) and the corresponding “scaling line” for Grímsvötn and other lakes listed in Table 4, alongside Q max −V t data of Figure 1 (circles or envelopes). Roman letters label some of the lakes appearing in Figure 1. c1 = Grænalón <1940; c2 = Grænalón >1950; d1 = Vatnsdalslón 1898; d2 = Vatnsdalslón 1974–78. The scales are based on n = 3 and K0 = 5×10−24 Pa−3 s−1.
According to Figure 7, scale dependence alone would cause a rising trend in the flood data, with a slope b close to 1. This arises because 
for n = 3 (Appendix), and the dataset happens to obey A
0 ∼ [V]2/3 (roughly indicative of geometric scaling) while Ψ has no significant statistical dependence on [V]. (In the Appendix, we also show that the result b ≈ 1 is insensitive to the value of n.) Reference ClarkeClarke (1982) was the first to propose a scaling argument. His mathematical derivation is a special case (of ours) where the same dimensionless closure rate is implied for all systems. This is equivalent to using a fixed value of γ in our model, which leads to [Q] ∝ Ψ11/6[V]4/3 (Clarke’s equation (29)) and a theoretical slope of 4/3. Our revised scaling removes the restriction on γ and is necessary for the cross-system analysis.
Because the slope predicted from scaling exceeds b ≈ 2/3 for the flood data, large systems (with large Vt
) release less violent floods (with lower Q
max) than is expected from simple scale extrapolation of what is observed for small systems. (We are referring to relative magnitudes.) This result points to systematic variation of other factors among the different lakes, specifically with lake volume. These factors include: (1) the thermomechanical parameters governing flood evolution (e.g. constants K
0 and n), and (2) the conditions at flood initiation, as measured on the dimensionless phase plane. (1) is motivated by the fact that the scaling slope will be altered if constants in the model vary systematically with [V] of the different lakes. (2) invokes variation of
(the ratio of Q
max to [Q]) with [V], without a change in the scaling slope.
In connection with (1), K 0 and n for a given lake system may be estimated by the method of section 4.1. Alternatively, they can be derived from known lake levels at the peak of flood events, because the corresponding Q−N pairs lie on nullcline (i) (the “Röthlisberger line”) and a simple regression would reveal K 0 and n needed to scale them correctly. These methods, however, may be used only if the peak occurs before the lake empties, which is not generally the case for floods involving complete draining of the lake. Despite this difficulty, it seems very unlikely that the different lake volumes should exert any control on the creep closure constants of subglacial drainage. Other parameters to be considered within the first factor include β and γ, but their statistical dependence on Vt is weak or inconsistent (Fig. 6), unable to account for much (and certainly not all) of the discrepancy between the predicted and observed values of b over a several orders-of-magnitude range in Vt .
We suggest that the second factor above — variation in the conditions of flood initiation among different lakes — is a probable candidate in the explanation of b ≈ 2/3. With the scaling slope unchanged, Figure 7 implies that
decreases with the flood volume (while approximately [Q] ∝ [V]), and thus large lakes (with large Vt
) are drained in the floods via lower trajectories on the phase plane when compared with what happens for small lakes (with small Vt
). If this is true, then the largest lakes may also exhibit incomplete-draining behaviour. This prediction is consistent with what we deduced from Figure 6 earlier for Grímsvötn and post-1950 Grænalón.
As yet, we do not fully understand the mechanism of flood initiation, but there is evidence that the dam-break process may be unique to each lake. The data clusters have limited range, and a triggering threshold in the lake level seems to operate in certain systems. A critical effective pressure for Grímsvötn has been mentioned before. Other examples include Grænalón and Vatnsdalslón, with N crit of about 4–5 and 2 bar respectively (Reference BjörnssonBjörnsson, 1988). Although Reference FowlerFowler’s (1999) theory can potentially connect these values to the physics of flood initiation, it needs to be tested against independent data or experiments, and its mathematical formulation adapted to the present context. Given the need to predict the timing of the floods as well as their magnitude, the initiation process is clearly the most important unknown to be resolved in future work.
5. Conclusions
The Clague–Mathews relation for jökulhlaups,
, is one of the more curious results in glacier hydrology. Its physical origin has remained obscure despite its apparent success as a peak discharge estimator to within an order of magnitude. For the first time, we give in-depth treatment of the data on which this relation is based, and we make a serious attempt at unravelling it within the framework of the Reference NyeNye (1976) model.
Phase-plane analysis of the model reduced to a pair of coupled equations clarifies how flood discharge–volume data (Fig. 1) ought to be interpreted. We propose a theory to show that for a given jökulhlaup system, there is in fact a monotonically increasing function, f, relating Q max to Vt , as long as draining of the lake is not complete. In this case, data from multiple flood events contain information on the floods’ evolutionary mechanics. The information may be extracted with known system geometry and lake-level history. Grímsvötn provides a unique example where we are able to demonstrate this in a direct comparison between model and data, and impressive agreement is found (Fig. 5). Thus, besides its success in simulating the flood hydrograph, Nye’s description of the floods in terms of the growth and collapse of channelized subglacial drainage is immensely robust in capturing their collective characteristics.
Existing data (Fig. 1) refer to outburst floods from a large number of lake systems, and the Clague–Mathews exponent b ≈2/3 is one that describes their rising trend among the lakes. The trend is the combined result of: (1) scale-dependent effects, due to the different system size and geometry and their modification through time, (2) the precise conditions at flood initiation, which may differ between systems and between floods, as well as (3) the fundamental flood evolutionary mechanics arising from the coupling of subglacial drainage to an ice-dammed lake. Although complete draining of the lake means that information on (3) is lost, making it more difficult to understand the characteristics of individual lake systems, this is not an obstacle when examining the other two contributions in a cross-system data analysis. We show that scaling alone leads to b ≈ 1 (Fig. 7), and conjecture that the observed b ≈ 2/3 at least partially reflects the conditions in which the floods are initiated in different systems.
And what are the implications for flood prediction? The Clague–Mathews relation has lost some but not all of its mystery, as we do not have an entirely physics-based theory yet at our disposal, but some of the components of this theory and the unknowns are becoming clear. Our study stresses the importance of continuously monitoring lake levels during flood events, and of making reliable field measurements, especially on ice-dam thickness (via radar-echo or seismic sounding techniques) and lake geometry. Flood initiation is the obvious missing piece in the puzzle. Combining Reference FowlerFowler’s (1999) theory (which provides a mechanistic description of seal-breaking) with the framework advocated in this paper (which connects model to observations in the large scale) seems to be a promising way forward. Such an approach calls for an extended assimilation of data from ice-dammed lake systems and a close exchange between jökulhlaup scientists worldwide.
Acknowledgements
We thank G. K. C. Clarke for directing us to his published Summit Lake bathymetry data, and J. Desloges for providing field data for Ape Lake. We appreciate the very helpful comments of J. Walder (Scientific Editor), R. LeB. Hooke and an anonymous reviewer on an earlier version of the paper. F. Ng acknowledges the support of a Junior Research Fellowship at St John’s College, University of Oxford.
Appendix
Scale definitions
If we define the constants
then
The scaling law and its sensitivity on n
To derive the scaling law discussed at the end of section 4.2, notice that the first definition in Equation (A2) implies
With the data in Table 4, it is possible to show that Ψ does not depend statistically on the lake volume scale [V] to have any significant effect here, while approximately A 0 ∼ [V]2/3. Therefore, we write [Q] ∝ [V] k where k = 8n/3(3n − 1). For n = 2, 3, 4 we obtain k = 16/15, 1, 32/33, respectively, hence the scaling slope is close to unity and varies little with n.
