## Abstract

The Nye model of jökulhlaups is able to explain both their periodicity and even the detailed shape of the flood hydrograph. In this paper, we show how the model can be used to predict the shape of the hydrographs from the subglacial lake Grímsvötn beneath Vatnajökull in Iceland, and we comment on three particular issues that have proved contentious in the application of the Nye model: the role of water temperature; the shape of the flood channel; and the value of the bed roughness coefficient.

## 1. Introduction

Jökulhlaups are floods that occur in streams that flow beneath glaciers. They are known to occur in a number of different regions, but perhaps the best known are those that occur in Iceland and particularly those that originate below the Vatnajökull ice cap. Of these, the floods that come from the caldera of the volcano Grímsvötn have provided a vivid example of the typical characteristics of jökulhlaups.

The Grímsvötn jökulhlaups have been thoroughly described by Björnsson (1974, 1992), and also, in particular, in the same author's book (Björnsson 1988), which is a definitive and compelling study of Icelandic jökulhlaups. The floods from Grímsvötn occur periodically, with a variable period of between 5 and 10 yr, as shown in figure 1. The shape of the flood hydrographs has been rather variable, as shown in figure 2, but typically (unlike subaerial river floods) there is a slow rising limb and a more rapid descending limb.

There are several questions of interest, which thus arise in describing the dynamics of jökulhlaups. These are as follows: why are the floods periodic; what determines the shape of the flood hydrograph; and what is the mechanism of flood initiation? The seminal paper in modelling jökulhlaups is that of Nye (1976), who used a time-dependent version of Röthlisberger's (1972) theory of subglacial channel drainage. Nye could neither solve his model nor answer the three questions above, although he did describe why the model should yield the characteristic flood hydrograph, and was also able to explain the shape of the rising limb of the 1972 flood hydrograph using an approximate version of his theory. Later, Spring & Hutter (1981, 1982) and Clarke (1982) provided numerical results for the whole hydrograph. Spring & Hutter (1981), in particular, tried to fit the whole 1972 hydrograph, with little success. They also emphasized the importance of the temperature equation in the model, and showed that the hydrograph shape is sensitive to the temperature of the lake water. Fowler (1999) showed that the Nye model could indeed explain the periodicity of floods, but he did not seek to obtain a quantitative comparison. Simply, the floods occur periodically because the lake continually fills with geothermally melted water until it overflows. One might expect that this would occur when the lake water pressure reached the overburden ice pressure, thus breaking the seal, but, in fact, Grímsvötn floods typically occur when the lake level is some 60 m below the flotation level. This had perplexed Nye, but can also be explained using his theory, as was shown by Fowler (1999). Fowler did not attempt to reproduce the shape of the hydrograph, pleading the excuse that the Nye theory, assuming as it does a semicircular channel, might not be able to do so.

Despite this success, Nye's theory has been criticized. Clarke (2003) raised two issues, and we can add a third, following Johannesson (2002*a*,*b*). Firstly, according to Clarke, the value that Nye used for the Manning friction factor (*n*′=0.09 m^{−1/3} s) is higher than that would be expected for a smooth-walled ice-roofed channel. Secondly, Clarke commended Spring & Hutter's (1982) model, and, in particular, their solution of the temperature equation. The temperature equation is of course also in Nye's model, but Nye did not attempt to solve it. Spring & Hutter also include acceleration terms in the momentum, but these terms are actually small, and their inclusion simply renders the problem stiffer.

A more serious possible criticism of the Nye model lies in its assumption that the subglacial channel is semicircular. Guðmundsson *et al*. (2004) described the dynamics of the enormous flood that followed the Gjálp eruption of 1996. In this flood (when the lake level reached flotation), the outlet channel from the lake spread downstream as a laterally growing, overpressured fracture, and the outlet channel reached a width of some 25 km towards the snout of the glacier. Both the propagation as an overpressured fracture and the formation of a wide channel are not quantitatively consistent with Nye's model, although his conceptual description is still appropriate. A more suitable model of the 1996 flood was provided by Flowers *et al*. (2004). A recent thorough review of the theory and observations of jökulhlaups has been given by Roberts (2005).

It has recently become clear that subglacial lakes are a common feature below ice sheets, or at least below Antarctica (Siegert 2005), and the possibility that floods could occur from them also arises. In fact, just such a flood was reported by Goodwin (1988) from surface observations, and, more recently, Wingham *et al*. (2006) have used satellite observations to infer the occurrence of a small intra-lake flood in the Adventure Trough below Antarctica. Just such floods were predicted by Evatt *et al*. (2006), essentially using the same Nye model as before. A particular feature of these floods is their relatively small amplitude. We will comment on this later.

## 2. The Nye–Röthlisberger model

The Nye (1976) model of subglacial drainage describes flow through a semicircular conduit from an origin at *x*=0 where the lake is, to an exit portal at *x*=*l*, where *x* is measured along the axis of the channel. The variables of the model are: *S*, the cross-sectional area of the channel; *Q*, the volume flux of water through the channel; *N*=*p*_{i}−*p*_{w}, the effective pressure of the channel equal to the difference between overburden ice and channel water pressures; *θ*_{w}, the water temperature; and *m*, the interfacial melting rate. These five variables are given by the equations(2.1)In these equations, *ρ*_{i} is the ice density and *K* is a creep closure coefficient due to a nonlinear flow law , where is the strain rate and *τ* is the stress; specifically, *K*=2*A*/*n*^{n}. The term *m*/*ρ*_{w} represents the volumetric source due to side-wall-derived melt, where *ρ*_{w} is the water density; *M* is an additional source (not included by Nye) due to tributary flow and surface meltwater supply; *f* is a friction factor; *g* is the acceleration due to gravity; and *α* is the mean bedrock slope (of the channel). The left-hand side of (2.1)_{4} is the material rate of change of water temperature with time, the first term on the right-hand side is the frictional source and the second is the supply due to the enthalpy change on melting. We take *θ*_{i} to be the melting temperature; *a*_{DB} (≈0.2) is a constant; *η*_{w} is the viscosity of water; *L* is the latent heat; *c*_{w} is the specific heat of water; and *k* is its thermal conductivity.

These equations represent, respectively, ice closure (the first equation is, in fact, the kinematic surface condition for the viscous motion of the ice), conservation of water mass, conservation of water momentum (ignoring small acceleration terms), conservation of water energy and the empirical Dittus–Boelter heat transfer correlation. The value of the Dittus–Boelter coefficient depends on precisely how the Reynolds number is defined. The Dittus–Boelter correlation was derived from experiments on flow through smooth pipes at Reynolds numbers between 10^{4} and 10^{5}. It is arguable whether the same formula (or at least the same value of *a*_{DB}) is appropriate at the higher Reynolds numbers of interest here.

The equations are to be supplemented by boundary conditions and subsidiary relationships, as follows. We define the basic hydraulic gradient *Φ* by(2.2)where(2.3)where *z*=*s* is the ice cap surface and *z*=*b* is its base.

The friction factor *f* in (2.1)_{3} is related to the Manning roughness *n*′ by(2.4)where *R*_{H} is the hydraulic radius (=*S*/*l*_{p}, where *l*_{p} is the wetted perimeter). For a semicircular channel , so that if *n*′=0.09 m^{−1/3} s, then *f*∼0.05 m^{−2/3} s^{2}. We follow Nye (1976) in using Manning's empirical formula, which was derived to describe open river channel flows. As such, there is no real justification to use it to describe closed channel flows, moreover consisting of two distinct parts: ice roof and sediment floor. In a sense, its justification lies in the veracity of the model results, as we shall see.

### (a) Boundary conditions

The equations in (2.1) require initial conditions for *S* and *θ*, two boundary conditions for *Q* or *N*, and one for *θ*_{w}. The boundary condition for *θ*_{w} is taken to be (at least when *Q*>0 at the lake)(2.5)At the outlet, it seems we should prescribe(2.6)i.e. the water pressure becomes atmospheric (as also does the ice pressure). At the inlet to the channel, conservation of mass requires that(2.7)where *V* is the lake volume and *m*_{L} represents its rate of increase, partly due to the supply of hydrothermal fluids from below, but mostly due to melting of the overlying ice shelf by the excess geothermal heat in the caldera. Suppose the lake level is at *z*=*h* and open to the atmosphere. We assume *V*=*V*(*h*), and, in fact, *V*′(*h*)=*A*_{L} is the lake surface area (which may depend on *h*). Now, the water pressure at the inlet is *ρ*_{w}*g*(*h*−*b*), where *b* is the grounding line elevation of the bed, and this is equal to *p*_{i}−*N*. Therefore (if *b* and *p*_{i} do not vary),(2.8)and thus the boundary condition at the lake inlet is(2.9)where *x*=0 is taken to be the position of the lake margin.

### (b) Non-dimensionalization

We choose scales for the variables by writing(2.10)where(2.11)in which *h*_{0} is the lake elevation above the outlet, and the scales are chosen as follows in order to balance appropriate terms in the equations, in terms of the yet to be specified *Q*_{0}:(2.12)

The dimensionless equations (2.1) then become (dropping the asterisks on the variables)(2.13)where the parameters *ϵ*, *r*, *Ω*, *δ* and *γ* are defined, after some algebra, by(2.14)

Finally, we choose the volume flux scale *Q*_{0} by balancing terms in the lake-refilling condition (2.9); again after some algebra, we determine(2.15)and the refilling condition (2.9) can be written in the dimensionless form(2.16)where(2.17)

Using values of the constants as set out in tables 1 and 2, we can determine typical values of the scales, as set forth in table 3, and the corresponding dimensionless parameters, as shown in table 4.

### (c) A reduced model

It must be emphasized that the derived parameter values in table 4 depend quite sensitively on the particular choices made in table 2. For example, if we change *n* to the value 3.5, then *Q*_{0} drops to 0.28×10^{4} m^{3} s^{−1}, and the value of *δ*, in particular, drops to 0.027.

We derive a reduced version of (2.13) by assuming *ϵ*≪1, *δ*≪1 and *γ*≪1. The choice *γ*≪1 is unrealistic and made for convenience; it has the precise effect of enabling the deduction that *θ*≪1, so that the terms in *θ* in (2.13)_{4} can be ignored. We discuss this further in §4. The other two assumptions are reasonable, and enable a significant reduction of the model. The reduced model, as presented by Fowler (1999), uses a redefined *x-*coordinate near the lake defined by(2.18)and explicitly assumes (as is appropriate for Grímsvötn) that the dimensionless basic hydraulic gradient *Φ* varies on the *X* scale; based on a fit to the observed surface and bed profiles, we take the dimensional basic hydraulic gradient to be of the form(2.19)where the subscript d denotes dimensional values. The corresponding dimensionless function is then(2.20)where(2.21)This choice of the basic hydraulic gradient was made as a simple fit to measured data from Grímsvötn by Fowler (1999). If, as is commonly done (e.g. Ng *et al*. 2007), one takes *Φ*=const. (i.e. *b*=0), then unbounded (and unrealistic) oscillatory floods occur (Ng 1998). The main point about *Φ* is that it must be negative (i.e. *a*<*b*) at *X*=0 in order that there be a lake at all; and it is a consequence of *this* that requires the resolution of the model near the lake, since we must include the terms ∂*N*/∂*x* in (2.13) when *Φ* is negative.

We can now write the reduced model in the form(2.22)where we define(2.23)The boundary conditions that are appropriate are(2.24)The first of these is the refilling condition, and the second is a matching condition to the far-field (outer) solution, which is obtained by neglecting *δ* in (2.13). The boundary condition *N*=0 at the outlet *x*=1 is lost, being relegated to a boundary layer near the snout.

In solving (2.22), we retain the small coefficients *ω* and *ν*; the first because it is only the sign of *Q* between floods that is important (in determining where flow reversal takes place), and the second because the term in *ν* is dominant in the refilling equation between floods, when *Q* is much smaller. The numerical method we use has been described by Fowler (1999). We solve equations (2.22) in the finite domain 0<*X*<*M*, where typically we take *M*=10. From the current time step, we step forward *S* and *N*|_{x=0} explicitly, and then use the condition *N*_{X}=0 at *X*=*M* to determine the divide position *X*_{0}(*t*) in the water flux equation *Q*=*ω*(*X*−*X*_{0}). We then solve for *N* by quadrature from *X*=0. The time stepping is then repeated using the improved Euler method to give second-order accuracy.

## 3. Fitting the hydrographs

### (a) The 1972 hydrograph

We now wish to ask whether this reduced model has the ability to reproduce the observations shown in figures 1 and 2. The basic answer to this is, surprisingly, yes. Figure 3 shows the variation of *N* at the lake (essentially figure 1 inverted) and the corresponding discharge at the outlet, as functions of time. Figure 4 shows the detailed shape of the hydrograph, with the data from the 1972 flood hydrograph. The agreement is astonishing, in view of the apparent arbitrary assumption of cross-sectional shape, and also the neglect of the effects of lake temperature, which has a significant effect on the hydrograph shape (Spring & Hutter 1981).

### (b) Parameter variations

In obtaining the fit in figure 4, we had to rescale the discharge and time scales. This raises the issue of how the model can be used to fit the observational data. To address this question, we consider that the choices of the scales and parameters in tables 3 and 4 are only indicators, and that we may improve the observational fit by judicious choice of some of the dimensional inputs to the model. We now discuss these choices sequentially. We are able to do this because the effects of the parameters are largely independent. First, we note that the variables of the model, such as discharge, can be written dimensionally as(3.1)where the notation indicates that the dimensionless flux *Q*^{*} depends only on the parameters *ω* and *ν*, as well as the choice of the dimensionless basic hydraulic gradient function *Φ*.

#### (i) Shape

The principal success of this model is in mimicking the hydrograph shape, at least for the 1972 hydrograph, which was the one that Nye (1976) focused on. We discuss other hydrographs below. The shape of the hydrograph is determined by the dimensionless function *Q*^{*}(*t*^{*}), and is independent of the scales *Q*_{0} and *t*_{0}. In fact, because *ω* and *ν* are small, they can be neglected during the flood, so that the hydrograph shape can only depend on the hydraulic gradient *Φ*. For fixed *Φ*, therefore, the flood hydrographs should all have the same shape. This is borne out in figure 5, where the hydrographs for differing values of *ω* and *ν* vary in shape only through the relative sizes of peak discharge compared with baseflow.

#### (ii) Maximum discharge and flood duration

Given that the hydrographs have a similar shape, any particular flood (with such a shape) can be fitted by choosing *Q*_{0} and *t*_{0} appropriately. There are a number of quantities in the definitions of the scales that are subject to uncertainty. These are, perhaps, the lake area *A*_{L} (which is, as a matter of fact, variable), the friction factor *f*, the Glen exponent *n* and the closure rate coefficient *K*. There are other constituent parts of the model, which could be adapted, for example the assumption of a semicircular channel, the assumption of Manning's law of friction and so on. To keep things simple, we take the two adjustable quantities to be *A*_{L} and *f*. With *n*=3, we have(3.2)Thus, in order to obtain the values of *Q*_{0}=1.57×10^{5} m^{3} s^{−1} and *t*_{0}=0.123×10^{−2} yr, i.e. multiplication by 2.7 and 0.493, respectively, we need to multiply the default values of *f* and *A*_{L} by factors Δ*f* and Δ*A*_{L}, where(3.3)from which we obtain(3.4)corresponding to adjusted values(3.5)both of which are acceptable, and, indeed, perhaps preferable, estimates for these quantities.

#### (iii) Period

The default period between floods is approximately 4 yr in figure 3, but, of course, the adjustment of time scale by the factor 0.493 also changes the period, to approximately 1.9 yr. Now between floods, *Q* at the channel inlet is very small, and the linear trend in *N*, which can be seen in figure 3, suggests that the period will be proportional to 1/*ν*. This is an oversimplification, since flood initiation depends on when the seal at *X*=*X*^{*}, where *Q*=0, reaches the inlet, and this depends on solving the model to determine *N*.

Asymptotic analysis of the inter-flood dynamics is accommodated by rescaling the variables as(3.6)and the equations become(3.7)subject to(3.8)where(3.9)

We can make the quasi-static approximation that *ϵ*′=0, and then we can eliminate *s*, and the inter-flood problem becomes that of determining the seal position *ξ*^{*} from the solution of(3.10)subject to (3.8) and (3.7)_{2}. Given *ξ*^{*}, there is a unique value of *Π*=*Π*^{*} there, such that (3.8)_{2} is satisfied. Solving backwards in *ξ*<*ξ*^{*} then provides *Π*=*Π*_{0} at *ξ*=0 as a function of *ξ*^{*}, so that (3.8)_{1} gives the evolution of *ξ*^{*}. In practice, the evolution of *Π*_{0} appears to be fairly linear (figure 3), which suggests that *λξ*^{*} is small, and the period should be largely independent of *λ*. In any event, we can adjust the period by changing *ν* while keeping *λ* fixed, or, approximately, keeping *ω*/*ν* fixed. So to change the period from 1.9 yr to the apparent 5.7 yr in figure 1, we can simply adjust the values of *ω* and *ν* to be *ω*=0.12×10^{−3} and *ν*=0.09×10^{−3}. Consulting (2.12), (2.14) and (2.17), we see that this can be done without compromising the values of the scales or the other parameters by reducing baseflow *Ml* and refilling rate *m*_{L} by a factor of 3. Indeed, one can use a variable refilling rate as a way of generating the variable period observed in the data. Figure 6 shows the result of doing this. The period is in fact 6.2 yr rather than 5.7 yr, at least partly because in order to compute the solution, we have had to lower the resolution to Δ*x*=10^{−3} and Δ*t*=10^{−4}; in addition, the value of ‘infinity’ where we impose the condition (∂*N*/∂*X*)=0 has been increased from *X*_{∞}=10 to 15.

#### (iv) Onset below flotation

One of the features of the normal Grímsvötn floods is that they begin when the lake level is below flotation (i.e. the inlet effective pressure *N*_{L} is positive). Consulting figure 1, we can see, assuming a flotation level of 1505 m, that initiation occurs when *N*_{L}≈6–7 bars, and flood termination occurs when *N*_{L}≈16–18 bars (1 bar=10^{5} Pa). This is in fairly good agreement with figure 3. The onset of floods relies on the migration of the seal, and the principal controlling effect on this appears to be the shape of the basic hydraulic gradient. In particular, a preferable fit to the measured data has *a*_{d}=290 Pa m^{−1}, *b*_{d}=1200 Pa m^{−1} and *c*_{d}=0.5 km^{−1} in (2.19), but with these values, the seal is strong and floods are initiated when *N*_{L}<0. To obtain the results in figure 3, we reduced *b*_{d} to 600 Pa m^{−1}. Such laxity here is motivated by the fact that, firstly, the exponential choice in (2.19) is a rather rough fit to the measured data, and, secondly, the data themselves are hardly likely to give an accurate representation of the hydraulic gradient along the flow path in the vicinity of the lake. For these reasons, we consider it reasonable to allow flexibility in the choice of hydraulic parameters.

If the minimum effective pressure is set by the hydraulic gradient, then the maximum is set largely by the dynamics of the flood. To see this, we might ignore the hydraulic gradient altogether by putting *Φ*=1, whence *S*=*Q*^{3/4}, and *N* and *Q* during a flood are functions of time only, and satisfy the approximate equations(3.11)Starting from a small value of *Q* at *N*=*N*_{−}, *Q* first increases then decreases until *Q*=0 at *N*=*N*_{+}, say, and the terminal *N*_{+} is determined in terms of the starting value *N*_{−}. Given a result such as that in figure 3, we could fine tune the maximum and minimum effective pressures by changing *b*_{d} (to fix *N*_{−}) and then vary the effective pressure scale *N*_{0}, for example by adjusting the closure rate coefficient *K*, in order to fix *N*_{+}. We do not pursue this further here, however.

### (c) Different hydrographs

We have shown how parameters can be methodically chosen to fit most of the characteristics of the jökulhlaups from Grímsvötn, except the shape, but that the 1972 shape is also well fitted by the model. The question arises whether other Grímsvötn jökulhlaups have the same shape. Figure 7 shows six of the hydrographs from figure 2, redrawn to have the same peak and duration. From this, we see that the rising limbs are generally similar, but that there is a marked disparity in the falling limbs, associated with the collapse of the channel.

Figure 8 compares these normalized hydrographs with a model output, and shows that the model does well. The worst agreement is with the 1972 slow decay and the 1934 fast decay. But we have already seen in figure 4 that the 1972 hydrograph can be fitted well, and the enormous magnitude of the 1934 flood suggests an eruptive origin as in 1996 (Guðmundsson *et al*. 2004); in this latter flood, the channel propagated as an overpressured lateral fracture, and the assumption of a semicircular channel is wildly inappropriate. In particular, the closure rate of such a wide channel would be expected to be much more rapid than for a semicircular channel.

## 4. Modelling issues

### (a) Friction coefficient

In his discussion of the Nye model, Clarke (2003) drew attention to the high value of the Manning roughness used by Nye (1976), and inferred that the model was flawed for this reason. Espousing the Spring–Hutter version of the model that includes the acceleration terms, Clarke preferred a lower value of *n*′, approximately 0.03. In fact, we have shown that the Nye model does not require a high value, and that the choice *n*′=0.04 m^{−1/3} s enables a quantitative fit to the 1972 hydrograph. Nor are the acceleration terms important. If we add the terms *ρ*_{w}(*u*_{t}+*uu*_{x}), where *u*=*Q*/*S*, to the right-hand side of the momentum equation (2.1)_{3}, then the corresponding dimensionless terms on the right-hand side of (2.13)_{3} are *Fr*^{2}(*ϵu*_{t}+*uu*_{x}), where we scale *u*∼*u*_{0}=*Q*_{0}/*S*_{0}, and *Fr* is an effective Froude number defined by(4.1)Note that this is based on the depth of ice, not the depth of the channel. From table 3, we take *u*_{0}≈13 m s^{−1}, whence we find *Fr*^{2}≈0.01.

### (b) Temperature variation

A serious criticism of the quantitative application of the Nye model is the neglect of the water temperature equation, in the sense that we assume that *θ*=0. Recent work (Ng *et al*. 2007) has shown that lake temperature of a marginal glacial lake can have a quantitative effect on flood characteristics. From (2.13), the energy equation can be approximated as(4.2)Evidently, the issue of temperature should only be of potential interest during a flood, and for the values of *θ*_{0}=3.5 K and *γ*=5.37, we would expect a significant effect. Clarke (2003) rightly highlighted the importance of temperature, as was also done by Spring & Hutter (1981) and Björnsson (1992), and yet measurements indicate that the exit temperature of the water remains close to the freezing point. To reconcile this with the theory, Clarke suggested that other heat transfer mechanisms might be appropriate, for example involving ice particles within the flow.

If we take the observation that *θ*≈0 at the outlet at face value, then it is only consistent with (4.2) if *γ*≪1, and the advection term is only negligible over the inlet region if *γ*≪*δ*. The value of *γ* relies on the use of the Dittus–Boelter heat transfer relationship, which was determined by measurements of heat transfer in fluids at Reynolds numbers of *O*(10^{4}–10^{5}). However, the Reynolds number appropriate here is , and it seems entirely feasible that the Dittus–Boelter relationship is inappropriate at such Reynolds numbers.

Until the necessity for the inclusion of heat transfer in the model can be demonstrated by the measurement of positive exit temperatures, it seems at least arguable that we can take *γ* to be sufficiently small that the heat advection term is indeed negligible. Additionally, it is worth pointing out that solving (4.2) between floods is a hazardous adventure, since *Q* reverses sign in the domain, and the direction of integration must reverse at the seal; in particular, no boundary condition can be applied between floods.

### (c) Cross-sectional shape

The assumption of semicircular cross section is perhaps the most blatant simplification that has no justification. The surprising thing is then that, despite this, the Nye model can do such a good job of simulating the flood hydrograph. The reason presumably is that the channel is indeed approximately semicircular, and in the Nye-type jökulhlaup, this seems reasonable, since the channel is inflated by frictional heating, which applies uniformly around the perimeter. It is equally clearly inappropriate for post-eruption floods such as that in 1996, when the overpressured lake water causes the channel to open as a propagating elastodynamic fracture, and the resultant channel is much wider. The consequent closure rate scales with the width not the depth, and this suggests an explanation for the rapid shutdown seen following such floods.

### (d) Sub-Antarctic floods

Recent observations by Wingham *et al*. (2006), Bell *et al*. (2007) and Fricker *et al*. (2007) have shown that subglacial lakes are not only common beneath Antarctica, but that they undergo rapid deflation, with typical surface elevation changes of the order of metres over times of the order of years. Such deflations (and inflations elsewhere) are interpreted as occurring through drainage of the lakes in subglacial floods, and the question arises whether such long duration floods of such small magnitude can be explained with the Nye model. The answer to this is positive. It can be shown that a year-long flood lowering the ice surface (and thus lake level) by metres can be simulated in the Nye model by increasing the ice closure rate by four orders of magnitude. The inference would be that the closure is due to subglacial till, and thus that such floods are excavated through high-pressure subglacial canals (Walder & Fowler 1994), rather than through Röthlisberger channels.

## 5. Conclusions

Nye's (1976) model of jökulhlaups has provided a theme around which several variations have been woven, and some of these variations have raised issues concerning Nye's model. In this paper, we have shown that Nye's model is able to provide an exceptionally good quantitative and qualitative explanation for the jökulhlaups that originate in the lake Grímsvötn below Vatnajökull in Iceland. Nye's model can explain the periodicity, the peak discharge and the duration of the floods, and it provides an excellent fit to the shape of the flood hydrograph.

There are two issues that are unresolved in this study. The first is that different Grímsvötn jökulhlaups have hydrographs that differ largely in the rapidity of the falling limb. It seems likely that this difference is a genuine one (as opposed to simply being due to poor estimated fluxes), and this suggests one way in which the model fails to represent the data. One possibility, suggested by the abrupt collapse following a post-eruptive flood, is that channel shape may play a role, since it is known that, following the 1996 eruption, the subsequent flood propagated through a very wide channel (Johannesson 2002*a*,*b*). However, it is very difficult to model the evolution of wide-channelled floods in an objective way.

The other issue is the role of temperature. Largely for pragmatic reasons, we have neglected heat advection in the energy equation (4.2). It seems simple enough to include it, but this is not true between floods, when *Q* reverses sign at the seal position *X*=*X*^{*}, and the consequently degenerate equation must then be solved in both directions away from the seal, where the degeneracy forces *θ*=0. While this is not insurmountable, the observation that exit temperatures are close to the freezing point suggests that, in reality, the value of *γ* should be much smaller than our estimate in table 4, and this may be the case if the heat transfer at Reynolds numbers in excess of 10^{8} is much more efficient than suggested by the Dittus–Boelter relationship. Alternatively, as suggested by Clarke (2003), other buffering mechanisms may keep the temperature near the freezing point. This point is worthy of further investigation.

Other objections, such as the arbitrary assumption of a semicircular channel, are *a posteriori* justified by the agreement of the theory with observations, except for post-eruptive type floods, where the initial channel growth is promoted by fracture propagation at the bed rather than melting. Clarke's (2003) comment on the high value of *n*′ used by Nye might be argued against on the basis of the extrapolation of the Manning roughness to enormous values of the Reynolds number, but, in fact, we have shown that a best fit to the 1972 hydrograph actually involves a lower value of *n*′=0.04 m^{−1/3} s in any case.

It remains to be seen to what extent the Nye model can be applied in a quantitative way to other jökulhlaups.

## Acknowledgments

I acknowledge the support of the Mathematics Applications Consortium for Science and Industry (www.macsi.ul.ie) funded by the Science Foundation Ireland mathematics initiative grant 06/MI/005.

## Footnotes

- Received November 26, 2008.
- Accepted February 12, 2009.

- © 2009 The Royal Society