## Abstract

Subglacial floods (jökulhlaups) are well documented as occurring beneath present day glaciers and ice caps. In addition, it is known that massive floods have occurred from ice-dammed lakes proximal to the Laurentide ice sheet during the last ice age, and it has been suggested that at least one such flood below the waning ice sheet was responsible for a dramatic cooling event some 8000 years ago. We propose that drainage of lakes from beneath ice sheets will generally occur in a time-periodic fashion, and that such floods can be of severe magnitude. Such hydraulic eruptions are likely to have caused severe climatic disturbances in the past, and may well do so in the future.

## 1. Introduction

Jökulhlaups, or subglacial floods, occur when an ice-dammed lake drains, and the resulting discharge can be very large. Sometimes such floods occur from ice-marginal lakes, and sometimes from subglacial lakes. Some of the best-known examples of jökulhlaups occur from beneath the Vatnajökull ice cap in Iceland which, sitting as it does astride the mid-Atlantic ridge, is subject to enhanced geothermal heat flow and occasional volcanic eruptions (Björnsson 1974, 1988). The outbursts from the subglacial caldera lake Grímsvötn are particularly well documented.

Floods have also been associated with ice sheets; in particular, the Channelled Scablands of Eastern Washington State in the USA are thought to have been caused by a succession (some 40 in all) of floods from the ice-dammed pro-glacial lake Missoula (Bretz 1923, 1969). These floods are thought to have been some of the largest on Earth, causing the gouged terrain that remains in place today. Similar terrain has been identified on Mars, and the apparent massive outflow channels there have led some to suggest that Mars used to have great ice sheets which discharged huge jökulhlaups (Baker & Milton 1974; Baker 2001; Chapman *et al*. 2003). There is also evidence from Antarctica of similar massive floods in the past (Sugden & Denton 2004).

### (a) Dansgaard–Oeschger events

One may ask what the climatic consequences of such floods would be. On the Earth, it is known that during the last ice age, there was a sequence of rapid climatic change events, in which the air temperature over central Greenland changed by 5–10 °C over a matter of decades (Johnsen *et al*. 1992; Dansgaard *et al*. 1993; Taylor *et al*. 1993). These events have been termed Dansgaard–Oeschger events, and it has been suggested that their cause was a comparable rapid change in North Atlantic Ocean circulation, itself driven by changes to the freshwater input to the North Atlantic (Stocker & Wright 1991; Ganopolski & Rahmstorf 2001). Essentially, the idea is that the North Atlantic circulation has different operating states: a strong, interglacial circulation and a weak, glacial circulation, each of them driven, ultimately, by the equator to pole differential solar heating. If excess fresh water is supplied to the North Atlantic, the resulting buoyancy may be sufficient to switch off a weak, glacial circulation entirely. The resulting static ocean is then so unstable that it will restart with a bang, causing the sudden warming which is observed. In all the model studies which support this view, the origin of this freshwater pulse is left unexplained, but the magnitude of the freshwater flux which is required is generally of the order of 0.1 Sverdrup (Sv), i.e. 10^{5} m^{3} s^{−1}. More specifically, model studies show that steady perturbations of freshwater flux to the North Atlantic of this order can cause hysteretic switches between different modes of circulation, and a sinusoidal variation of freshwater flux with amplitude of only 0.03 Sv over a period of 1500 years can cause model Dansgaard–Oeschger events to occur (Ganopolski & Rahmstorf 2001).

A further suggestion is that there is a rough periodicity to the Dansgaard–Oeschger events, with the interval between them being about 1500 years (Bond *et al*. 1997). These events are clustered in cycles, and terminate in Heinrich events (Heinrich 1988), which are thought to occur because of massive discharge of icebergs from the Laurentide ice sheet every 10 000 years or so (Bond *et al*. 1992; MacAyeal 1993).

If this story is true, it remains to establish the cause of the anomalous freshwater fluxes, and the reason for their apparent 1500 year periodicity. The most obvious difference between the glacial Dansgaard–Oeschger events and the present apparent lack of them (although Bond *et al*. (1999) suggest that they continue in muted form, for example, in the Little Ice Age) is the presence of massive land ice, particularly in North America and Europe. The most obvious source of water is from subglacial or pro-glacial lakes, and the most obvious source of periodicity is that of a jökulhlaup. Hence, we adduce to the theory of ocean circulation fluctuations the idea that there were massive outburst floods from beneath one or more of the ice sheets on the Earth at that time.

### (b) Jökulhlaups

In the absence of direct corroborative evidence for this idea, we need a viable theory which suggests that it is plausible. This theory exists, and has been frequently applied to the study of jökulhlaups. The basic theory is due to Nye (1976), who used the hydraulic theory of Röthlisberger (1972). This theory has subsequently been elaborated by a number of authors, including Spring & Hutter (1981, 1982), Clarke (1982, 2003), Björnsson (1992) and Fowler (1999). A recent thorough review is by Roberts (2005).

The Nye theory is incredibly successful in predicting both the shape of the flood hydrograph and the flood periodicity. It has not previously been applied to the possibility of floods from beneath ice sheets, and it is the central purpose of this paper to provide a coherent application of this theory. Since our idea is that such floods should provide a robust explanation of Dansgaard–Oeschger events, we pose as goals the identification of a millennial time scale and a peak discharge of order 0.1 Sv. It turns out that the time scale can be easily identified, but the peak discharge requires a genuine solution of the Nye model. We also examine possible floods from Lake Vostok. We show that such floods are possible and have slightly larger peak discharges than those from Grímsvötn, but the periods are much longer, being of the order of 10^{4}–10^{5} years.

The idea of floods from below ice sheets is not a new one. In particular, Shaw and his co-workers (Shaw 1983; Shaw *et al*. 1989) have been keen on positing the idea of massive subglacial floods in order to explain the origin of drumlins. The idea of floods causing subglacial landforms is also not new, and dates back to Hall (1815)—before the glacial theory! In the spirit of the times, Hall's theory was placed in a scriptural culture.

Shaw's pursuit of his thesis has been controversial: a recent defrocking is that by Clarke *et al*. (2005). Despite this, it has to be said that many recent discoveries are consistent with Shaw's viewpoint: the idea that the 8200 b.p. cooling was caused by a massive flood from Lake Agassiz (Alley *et al*. 1997; Clarke *et al*. 2004), the wide channel propagation of the 1996 Gjálp eruption and flood from Vatnajökull (Guðmundsson *et al*. 1997, 2004), the discovery of many subglacial lakes below Antarctica (Siegert *et al*. 1996, 2001; Siegert 2005). While not wishing to eschew imaginative possibility, in the present paper, we seek to establish possibilities that are securely based on reliable physics.

In the remainder of the paper, we provide a version of Nye's hydraulic theory which is appropriate to lakes below ice sheets, in particular, where the lake is not connected to the ice surface and the ice over the lake responds dynamically to lake level fluctuations. The resulting model is solved numerically for the case of a single lake draining to the sea, and we show that floods generally occur. For a lake beneath the Laurentide ice sheet with a large catchment area, their period can be of the order of millennia, and their peak discharge, obtained over several years, can be of the order of 0.1 Sv. It is on the basis of these findings that we suggest that sub-ice sheet floods are a viable explanation for Dansgaard–Oeschger events. We then apply the results to the possible catatrophic drainage of Lake Vostok. The catchment for this lake is much smaller, and consequently periodic floods, if they occur, should have a period in the range 10^{4}–10^{5} years.

## 2. The Nye theory of subglacial floods

### (a) Channel hydraulics

The mathematical theory of jökulhlaups is given in a seminal paper by Nye (1976), based on the earlier theory of subglacial hydrology of Röthlisberger (1972). In this theory, a flood occurs by drainage from the source region, usually a lake, along a single channel. Fundamental to the concept of the theory is that the channel water pressure is lower than the overburden ice pressure , which allows the channel to maintain its integrity and not flood the bed. This comes at a price, because the resulting positive effective pressure, defined as(2.1)will cause the channel to close through viscous creep of the ice. Röthlisberger's concept was that this closure could be balanced by an opening rate caused by meltback of the channel walls. The heat to cause this meltback is supplied by the viscous dissipation of the turbulent water flow in the channel.

This radical theory has been very successful in describing steady drainage, and it has proven to be equally successful in describing jökulhlaups. In Nye's formulation, the channel cross-section *S* is thus determined by the combination of creep closure and meltback, which leads to the equation(2.2)where *K* is proportional to the rate coefficient *A* in Glen's law(2.3)in which is the second strain rate invariant, *τ* is the second stress invariant and *n* is the Glen exponent. Specifically, . The other quantities in equation (2.2) are the melt rate *m* and the ice density . We discuss an appropriate modification of this equation below.

Equation (2.2) is derived as the kinematic condition of the ice–water interface, assuming a locally two-dimensional closure rate and, in order to obtain a closed form solution, a cylindrical channel in an infinite expanse of ice. One usually thinks of the channel as semi-circular, when the same equation will apply, provided the glacier bed is essentially stress-free. Some thought has been given to the task of modelling non-circular channels (Ng 1998), but with limited success. It is possible that the semi-circular channel description is accurate during major floods, when the rapid channel opening during flood initiation is essentially isotropic.

Equations for the hydraulics of the channel are taken to be(2.4)These represent conservation of mass and momentum of the water flow. In them, *x* is the distance along channel, *Q* is the volumetric water flux, is the density of water, *M* is a source term to represent the inflow of tributary water from elsewhere on the bed; *g* represents gravity, *α* is the local angle of inclination (downhill) of the bed and *f* is a friction factor. Equation (2.4)_{2} represents the empirical Manning law used by Nye; *f* is related to the Manning roughness (designated here) by the relation(2.5)where is the hydraulic radius (, where *l* is the wetted perimeter). For a semi-circular channel , so that if, for example, (a rough channel), then . Clarke (2003) has queried the use of the Manning law, and in particular the high roughness coefficient. His argument stems from the idea that the high values of used by Nye (1976) were required because of an incorrect assumption that thermal advection in the channel flow was small. This assumption is discussed further below. A consequence of the retention of thermal advection found by Clarke (2003) is that exit water temperatures are higher than suggested by observations. To accommodate this, Clarke suggests the possibility (see also Björnsson (1992)) that heat transfer may be more effective at the very high Reynolds numbers of jökulhlaups. As we indicate below, a consequence of this suggestion is the validation of Nye's original assumption.

The momentum equation (2.4)_{2} also neglects the inertia terms. Spring & Hutter (1981, 1982) introduced these, as did Clarke (2003). Despite this, it is straightforward to show that these terms are small, and hence their inclusion suffices only to cause numerical headaches, as the problem becomes stiffer.

The energy balance for the water flow is given by(2.6)Here, is the water temperature, is the ice temperature (which we suppose to be the pressure melting point), is the specific heat of water, and *L* is the latent heat. The first term on the right-hand side is the frictional heat generation of water flow and the second is the latent and sensible heat necessary to melt the ice wall and raise its temperature to that of the water. A final equation relates this extra heat to that supplied by turbulent heat transfer from the water, taken by Nye to be the empirical Dittus–Boelter relation (Szekely *et al*. 1976),(2.7)In this last equation, is a numerical constant, is the viscosity of water, and *k* is its thermal conductivity. Heat transfer has been something of an issue in modelling jökulhlaups (Björnsson 1992; Clarke 2003). Use of equation (2.7) tends to give exit water temperatures which are significantly higher than the freezing point, while measurements indicate emerging floodwater is essentially at the freezing point. One simple reason for this may be that the Dittus–Boelter relation was determined by measurements at Reynolds number approximately 10^{5}, whereas the Reynolds number is much higher in jökulhlaups, perhaps of order 10^{8} (Nye 1976), and it is not obvious that the Dittus–Boelter relation can be extrapolated to these higher flow rates (Björnsson 1992). A simple way of representing this possibility is to retain equation (2.7), but use a higher value of . The five equations in (2.2), (2.4), (2.6) and (2.7) provide a sufficient model to determine the five variables *S*, *Q*, *N*, *m* and , providing the other quantities in the equations are prescribed.

#### (i) Boundary and initial conditions

The schematic geometry of the situation we wish to consider is indicated in figure 1. We suppose an ice sheet flows over a depression in the bedrock, so that a subglacial lake exists there through entrapment of basal meltwater if the basal ice reaches the melting point, and there is net melting of water. Examples in Antarctica include Lake Vostok (Siegert 2005) and candidate lakes under the Laurentide ice sheet would be the Great Lakes, as well as Great Slave Lake and Great Bear Lake.

We consider only flow in two dimensions, with *x* being horizontal and *z* vertically upwards. We denote the upstream and downstream margin positions of the lake by and . With water flowing from left to right, there will be a volume flux at into the lake and a volume flux at out of the lake. (If then water at flows into the lake.)

The conditions we appear to require are initial conditions for *S* and , an upstream boundary condition for and two boundary conditions, one upstream and one downstream, for the water flux *Q* and the effective pressure *N*. We require these conditions on each of the segments , , where we locate the ice summit at *x*=0 and the ice margin at *x*=*l*. In figure 1, the ice sheet margin is taken to be at the grounding line.

On the upstream segment, we have(2.8)Because of this, the temperature equation is degenerate, and no boundary condition for can be provided: the water temperature at the summit selects itself.

We suppose the lake is sealed from the ice surface, so that the weight of ice over the lake is supported by the water pressure. In fact, flexure of the ice will depend on the lake pressure, and the dynamics of this will determine the effective pressures at the margin. For the moment, we specify(2.9)and will comment further on the choice of (and particularly ) later on.

At the ice sheet margin, be it land-based or at an ice shelf, it is appropriate to pose(2.10)In addition, we take the temperature at the inlet to the downstream channel to be that of the lake temperature, which we suppose known, thus(2.11)

As the lake drains in a flood, the margin points may move, and their determination then forms part of the problem to be solved. One condition that is associated with this is the lake refilling equation(2.12)where *V* is the lake volume. A second condition follows from a consideration of the ice dynamics over the lake. This is discussed further below in §2*c*. In this paper, we will make the simplifying assumption that the lake roof is flat, in which case and are geometrically determined functions of *V*.

### (b) A reduced hydraulic model

In order to provide a tractable model, we begin by non-dimensionalizing the Nye model. First, we write the hydraulic gradient in equation (2.4)_{2} in terms of the effective pressure. If the bed elevation is *b*, then the bed slope is , while the cryostatic pressure at the base is , where *s* is the ice surface elevation. In terms of these, we then have that(2.13)where the *basic hydraulic potential* (i.e. that constructed assuming basal water pressure equals ice overburden pressure) is given by(2.14)

To non-dimensionalize the model, we scale *Q* with , *S* with , *m* with , with , *t* with and *N* with , and we choose these scales to effect suitable balances of the dominant terms in the equations. The choice of is made to balance lake emptying rate with discharge, thus we choose(2.15)where is a measure of total lake volume (we make this definition more precise later). The other scales can then be determined in terms of , as follows:(2.16)in which we suppose is a representative ice thickness scale, and(2.17)Solving for , we find that the flux scale is given by(2.18)The dimensionless versions of the hydraulic equations are then(2.19)We non-dimensionalize surface and bed elevations *s* and *b* using the depth scale , with the result that the dimensionless hydraulic gradient is defined by(2.20)There are six dimensionless parameters in the model: and *γ*, and these are defined, after some algebra, by(2.21)

In order to simplify the model, we need to estimate the sizes of these parameters. To do this, we use the values , , , , , , , and (based on equation (2.3)) with *n*=3 and at 0 °C (Paterson 1994). As discussed after equation (2.5), we estimate , noting that lower values (by a factor of about 5) have been advocated by Clarke (2003). Suitable choices of length and depth scales for the Antarctic ice sheet along the flow line from Vostok to the coast are *l*=1500 km, , whence we find , and we suppose that , which is the inferred volume of Lake Vostok (Siegert 2005). For hypothetical floods of sub-Laurentide lakes, a choice of values is less certain. With Lake Ontario (volume 1600 km^{3}) in mind, we use Laurentide values of *l*=750 km and (cf. figure 1).

Using these values (and the Vostok length scales), we then find successively (using a Glen exponent *n*=3) that(2.22)and the dimensionless parameters are of typical sizes(2.23)The value of is based on the assumption that base flow rate of the outlet stream from Vostok is *Ml* m^{3} s^{−1}, and that this is equal to the product of the average basal melt rate and the catchment area for the Vostok stream.

We calculate the catchment areas below the ice sheet which will drain into Lake Vostok and into the Vostok stream based on a simple analysis of the Bedmap dataset (Lythe *et al*. 2001). We use ice thickness and bed elevation values to calculate hydraulic potential at the bed assuming that basal water pressure is ubiquitously at ice overburden pressure. We then derive flow paths across the hydraulic surface for the whole domain using ArcInfo GIS software. Based on this analysis, we calculate the ‘upstream’ area of the subglacial basin that would feed water in to Lake Vostok across these potentials, that is, we calculate the size of the area in which all the hydraulic flow paths lead into the lake, and we do the same for the hydraulic flow paths into the Vostok stream, which itself is a hydraulic grade line from Vostok to the coast. On this basis, we compute the catchment for Lake Vostok to be 28 000 km^{2} and the catchment to the Vostok stream to be 765 000 km^{2}.

To compute the basal melt rate , we use the formula(2.24)where *G* is geothermal heat flux, is basal shear stress, is basal sliding velocity, *k* is thermal conductivity, is the basal vertical temperature gradient, is ice density, and *L* is latent heat. We use values *G*=65 mW m^{−1}, , , , , , , based on a surface to base temperature difference of 50 K and a depth of 3000 m. The quantity *ϕ* is an enhancement factor which is one for pure (steady) conduction, but greater than one when advection is important. In accordance with thermal boundary layer theory, its size should be , where the reduced Péclet number *Pe* is given by , and using previous values together with a value for the thermal diffusivity *κ*=38 m^{2} yr^{−1}, we find . This is of course an upper bound. Substituting these values into equation (2.24), the constituent terms are a geothermal melt rate approximately 6.5 mm yr^{−1}, a frictional melt rate approximately 3 mm yr^{−1} and a thermal freezing rate in the range 3.5–14.5 mm yr^{−1}, although if the ice at the base becomes temperate, this freezing rate drops to zero. A plausible range for melt rate is thus 0–10 mm yr^{−1}, with the effect of advective cooling suggesting a low estimate within this range. If, for example, the basal melt rate is 2 mm yr^{−1}, then the normal Vostok stream discharge is or *Ml*=48 m^{3} s^{−1}.

#### (i) Simplified hydraulics

The parameters , , *δ* and are all relatively small, and we seek to ignore these where possible. An appropriate simplification was considered by Fowler (1999), who showed that the Nye model had solutions representing periodic jökulhlaups.

The neglect of is a regular perturbation (there is of course a rapid transient associated with the initial condition for *θ*, but this is not of concern), and we may safely put . We cannot, however, neglect the very small tributary flux term , because this term allows for the presence of a seal in between floods, as discussed later.

The temperature equation requires some discussion, partly because the effect of temperature has sometimes been ignored (Fowler 1999), although it can play an important part in quantitative solutions (Spring & Hutter 1981). Equation (2.19)_{5} really defines *m*; substituting this in equation (2.19)_{4}, neglecting the term in and using (2.19)_{3}, we have(2.25)We can see from this that *θ* relaxes to exponentially on a length scale of . Between floods, this distance is really very small as is the dimensionless temperature, and it is reasonable to suppose . During a flood, the equation above suggests that (where here and subsequently, *x* and (later) *X* subscripts denote partial derivatives) over a distance , and as a consequence the melting rate(2.26)will depend significantly on the water temperature, unless (discussed later) .

The small coefficient *δ* represents a singular perturbation and cannot be neglected in the vicinity of the lake. In fact, since the presence of a lake indicates a local hydraulic minimum, we must have near the lake, so that neglect of the term in *δ* in the momentum equation would imply *Q*<0, and the lake would never drain. Floods occur in this theory because of the existence of a seal, which allows backflow towards the lake. The position of the seal migrates slowly backwards towards the lake, and a flood is then initiated when it reaches the lake. Hydraulically, the lake overflows.

Away from the lake (where ), the limits , and allow an outer solution in which , and *N*, *Q* and *S* are related by(2.27)The value of is determined by matching this approximation to the flow near the lake, in which we cannot neglect the small term in *δ*. Specifically, we rescale the space coordinate as(2.28)so that(2.29)where(2.30)

The unknown marks the position of the seal. When it is positive, the lake is impounded, and gradually fills because of inflow at . As long as remains positive, channel outflow at the ice margin is small, but as soon as reaches zero, a flood is initiated, and the lake will drain. The ‘hydraulic’ length scale is about 135 km, comparable to dimensions of Lake Vostok, and we will generally suppose that also varies on this hydraulic length scale near the lake, but varies more slowly away from the lake. Thus, we suppose that as *X* becomes large away from the lake, , and is constant. In this case, given *S*, suitable conditions for *N* in equation (2.29) are(2.31)The prescription of the marginal effective pressure will be discussed later. The second condition as is a matching condition which ensures that the solution near the lake blends smoothly with the outer solution far from the lake. The extra condition on *N* thus determines in the outer solution as we suppose, and it is this which determines the seal position . From equation (2.19)_{4}, we have (neglecting terms of ) , where we define , and thus *S* is stepped forward by solving(2.32)where satisfies(2.33)if *Q*>0, and then(2.34)assuming the lake temperature to be at the melting point. We suppose that if *Q*<0.

The extra term, , arises because of the change of variable (2.28). It is negligible between floods. Numerically, it is irrelevant, because even if we allow the lake margin to move, we solve the hydraulic equations on a fixed grid, allowing the inlet portal to migrate along this grid.

### (c) Ice sheet flow and portal mechanics

In a previous study of Grímsvötn (Fowler 1999), the hydraulic model was self-contained. The physical basis of this model was the idea that the ice cover was broken at the caldera walls, thus allowing the lake to be open to the atmosphere, as is sometimes observed (Björnsson 1988). As a consequence, lake volume (and thus discharge) were directly related to the marginal effective pressure.

For sub-ice sheet lakes, we do not expect the water to be exposed in this way, and a proper consideration of the boundary condition for the lake discharge requires a further analysis of the over-lake ice dynamics. Figure 2 shows a cartoon view of the possible evolution of the ice surface during a flood. As the lake discharges, the ice is drawn down, and this requires an underpressure in the lake so that the ice can deform. If the aspect ratio of the over-lake ice shelf is large, then the appropriate approximate theory to describe this deformation is the viscous analogue of beam theory (Howell 1994). Viscous beam theory provides a relation between the surface velocity *w* and the effective pressure *N*. If, as shown, the ice and lake surfaces remain flat in the interior, then the surface deformation rate over the lake, where *z*=*h* denotes the lake roof, and(2.35)where is the lake roof area, a geometrically determined function of *h*. The relation between *N* and *w* provides a generalized version of the lake refilling equation which has been used in previous studies. A possible complication is that for the large floods with which we are concerned, the vertical collapse may be hundreds of metres in a matter of years. The stresses necessary to produce such rapid collapse may be much larger than the fracture strength of ice, and as a consequence, the ice may fracture near the lake margin. As collapse continues, a series of such ring fractures will form, allowing the central part of the ice to collapse passively, as indicated in figure 2. Such ring fractures are commonly seen on small ice cauldrons on Vatnajökull (e.g. Guðmundsson *et al*. 1997, 2004).

We will not provide the viscous beam theory here, although it has been done and will be reported elsewhere. A consequence of the analysis is that the open lake assumption still appears to be reasonable. In this case, we could suppose that the lake margin positions and remain fixed. Additionally, we suppose that as the lake empties, the marginal ice collapses or fractures in such a way that the collapsed ice blocks in the lake do not provide much resistance to the water draining towards the outlet portal (cf. figure 2). The interior lake water has the uniform dimensional pressure (as it did initially) assuming the lake roof is flat, and is in hydrostatic equilibrium with the water at the portal entrance. It follows from this that at , where the elevation of the lake roof is , the dimensional effective pressure is(2.36)and hence from equation (2.35),(2.37)This is the version of the lake refilling equation that is normally used. By our choice of volume flux scale, this can be written in the dimensionless form(2.38)where now represents dimensionless lake area (scaled with a typical value , for example that when full under basic hydraulic conditions), and we have chosen a lake depth scale to be(2.39) is, in fact, the order of magnitude of the drawdown of the lake during a flood. Thus, the ‘typical’ lake volume change during a flood is(2.40)A lake will empty during a flood if its volume is less than (approximately, since is simply a scale).

The parameter *ν* is then defined by(2.41)*ν* is the dimensionless inflow rate, which we take to be constant (and small). For , we have . This definition of allows us to talk sensibly about volume changes in an infinitely deep lake, in which floods will occur but the lake will not empty.

We call the above model a Grímsvötn type model. An alternative supposition, which we call a Vostok type model, is that as the marginal ice blocks collapse, there is severe impedance of the flow, and the channel must melt its way back with the lake margin. In this case, varies and we apply(2.42)As the lake discharges, we suppose that the central part of the over-lake ice remains flat, so that (and ) are determined kinematically in terms of lake volume. In order to solve equation (2.32) for the channel cross-section, we assume that *S* is continuously differentiable at .

#### (i) Ice sheet flow

Another consideration of importance is the recovery of the ice sheet between floods. The drawdown of the over-lake ice causes the ice to flow back in towards the lake, and in the (Vostok) case of a moving ice margin, it is this drawdown which causes the seal to reform, since it provides a basic hydraulic gradient back towards the lake. We have examined two ways of modelling the ice sheet flow. One is to take an isothermal shallow ice flow model (e.g. Hindmarsh 2001), which takes the form of a nonlinear diffusion equation. This model uses Glen's flow law together with the shallow ice approximation (Fowler & Larson 1978).

A simpler model for the ice flow assumes perfect plasticity (Nye 1951; Van der Veen 1999), when the ice is assumed to be shallow and perfectly plastic with a prescribed yield stress. In both models, the ice flow is suspended while a flood occurs, and the ice sheet is allowed to draw down kinematically as illustrated in figure 2. We have not found any serious differences between models in which either plastic or viscous ice flow is included.

#### (ii) Portal closure

A feature of the Vostok type portal model is that *N*=0 at the channel entrance. According to the Nye model (2.2), this allows indefinite growth of the channel cross-section there. Some further consideration of the reason for this is necessary.

Nye (1976) based his closure equation on Röthlisberger's (1972) earlier use of his borehole closure theory (Nye 1953). This assumes a circular void in an infinite expanse of ice. It is easy to modify this theory (Evatt in preparation) to allow for an externally imposed pressure at a finite distance, and even (with some licence) to allow an external surface at finite distance having no stress, but with a radial body force mimicking gravity. These two problems give similar results, and remove the unbounded growth that occurs in equation (2.2) when *N*=0. The modification which this leads to can be written in the form(2.43)where is the finite area of the enclosing circular ice. The dimensionless version (2.32) is modified to(2.44)where we have scaled with . The point now is that opening will indeed occur if , but when *S* reaches , it remains there.

## 3. Numerical results

It can be seen that there is a variety of ways in which we can choose to include ice dynamics and portal conditions in the model. We have implemented different versions of these, and conclude that the behaviour of the model is essentially similar no matter what choice we make. The exception to this statement is the effect of temperature, where the inclusion of the advective term, , in equation (2.32) can have a significant effect on flood discharge, if *γ* is not small. We do not report any results here in which the effects of thermal advection are included, ostensibly on the basis that the parameter *γ* in equation (2.33) *is* small.

### (a) Flood magnitude and frequency

Figure 3 shows the hydrograph of a typical flood from a sub-ice sheet lake resembling Lake Vostok. As we have commonly found to be the case for hypothetical subglacial lakes of such large magnitude, the floods are very large, and occur periodically with a very long period. As we discuss later, the period is essentially controlled by the rate at which the lake is refilled, and for the case of Vostok, the period may be of the order of 40 000 years, or even longer. The long period is due to the relatively small (28 000 km^{2}) catchment area of the lake. The peak discharge is about 10^{5} m^{3} s^{−1}, and the total discharge is some . In this simulation, the lake does not empty. With the current configuration of the ice sheet, the drainage path is expected to lie in the direction of Byrd Glacier in the Transantarctic Mountains. Figure 4 shows the present day hydraulic drainage paths, and figure 5 shows the ice surface and bed along the hydraulic flow line from Vostok.

The principal distinction between the simulated floods from a Laurentide-type lake (we use Lake Ontario as a representative) and a Vostok-type lake is in the assumed size of the catchment area. The basic hydraulic gradient is similar, but we suppose the catchment for Lake Ontario is much larger. Its current sub-aerial catchment is 64 000 km^{2}, and under an ice sheet this is likely to be larger, since the ice surface would have sloped towards Lake Ontario from some distance, and the basal topography is quite flat. We have used a value of 450 000 km^{2}. This value was computed by analysing the upstream contributing area of the hydropotential surface that should feed water to subglacial Lake Ontario. We take this as a typical value, while recognizing that for different lakes there will be a good deal of variation around this figure.

The result of the larger catchment is that the period is much shorter, some 800 years. The magnitude of the peak discharge is smaller, but the floods are of longer duration, some 10 years, as opposed to the year long floods from the Vostok type lake, and the total discharge is similar, about . Figure 6 shows a sequence of these floods, which appear as spikes because of the brevity of the floods relative to the period. Figure 7 shows the detail of the hydrograph of one of these floods.

## 4. Discussion

In our numerical studies of drainage from sub-ice sheet lakes, we have found that periodic flooding is the normal drainage mechanism, and this seems to be the usual case when the parameter *ω* is small. The magnitude and duration of the floods appear to be controlled by the channel hydraulics, but the peak discharge is essentially related to the lake volume; big lakes produce big floods. A question of interest is whether lakes will be emptied during a flood. This depends on the depth of the lake. In fact the drawdown during a flood is given by the quantity in equation (2.39), so that a first-order estimate for total flood discharge is(4.1)where is lake area and is effective pressure level in the channel. Incidentally, this suggests that the small floods apparently observed by Wingham *et al*. (2006) are associated with drainage through channels at low effective pressure.

The interval between floods is, fairly obviously, controlled by the refilling rate. This interval is simply the total flood discharge divided by the refilling rate, which itself is the product of the basal melt rate and the catchment area; thus a simple estimate for the period is just(4.2)where is lake catchment area and is basal melt rate. A small catchment, such as is thought to be appropriate for Vostok, leads to periods of tens of thousands of years. If we take , , then a basal melt rate in the range 1–4 mm yr^{−1} yields estimates for periodic discharge between 5000 and 21 000 years. However, we have also modelled flow paths during ice maximum conditions such as might result from lowered sea levels during glacial periods. The rearrangement of the ice surface under such conditions leads to water being forced in the direction of Law Dome, and a roughly fourfold extension of the Vostok basin to about 110 000 km^{2} (see figure 4). Under such conditions, refilling rates and flood frequencies would be likely to decrease accordingly. The larger catchments, which may also apply to putative Laurentide lakes, give shorter periods, but equally large total discharges.

An issue which we have not addressed, but which we must at least mention, concerns the rôle of melting and refreezing of the lake roof ice. In our discussion of the basal drainage system from the lake, we have assumed that inflow to and outflow from the lake are entirely due to supply and removal. In general, this is unlikely to be true, and a further source or sink is that due to melting or refreezing of the roof ice over the lake; indeed, this is known to be a dominant thermal process under ice shelves (Jenkins 1991; Holland & Jenkins 1999). A simple calculation of a Rayleigh number for a lake of depth 100 m subject to a geothermal heat flux of 50 mW m^{−2} gives a value of , so that thermal convection is inevitably vigorous, and melting and refreezing can be expected; indeed, this is known to occur under subglacial lakes Vostok and Concordia (Siegert 2005; Tikku *et al*. 2005). It is clear that a considered view of flood periodicity and magnitude requires an understanding of rates of subglacial melting and refreezing, and it is also clear that theoretical assessment of this process is of some difficulty.

## 5. Conclusions

Our central finding is that lakes below ice sheets will normally drain by flooding. The floods can be of large magnitude, may empty the lake, and occur at periods, which are of the order of the lake volume divided by the lake refilling rate. The refilling rate is (if we ignore melting and refreezing at the lake roof) simply the catchment area for the lake multiplied by the basal melt rate. This principally varies with the catchment area. For a central lake such as Vostok, the catchment area is relatively small, and the estimated period is of the order of 10^{4} years (it should be emphasized that this is a crude estimate).

As well as the present theoretical suggestion of massive floods from sub-Antarctic lakes, there is direct geomorphological evidence for their occurrence (Denton & Sugden 2005), there is direct observation of outlet floods occurring from Antarctic subglacially derived meltwater ponds (Goodwin 1988), and there is recent observation of rapid ice surface elevation changes which essentially imply internal subglacial flooding (Wingham *et al*. 2006).

Antarctica has mountainous basal terrain, and is littered with many small subglacial lakes (Siegert 2005). We conceive of the possibility that drainage from below the ice sheet occurs by many episodic drainage events, so that the basal drainage system resembles that of a Darcian flow in a porous medium, on a grand scale. As the lakes drain, the ice sheet surface will pump up and down, and on the time scale of the drainage events, the surface will display the appearance of a pan of bubbling milk.

If we shift our attention to the former Laurentide ice sheet, it is natural to enquire whether similar floods could have occurred there. Evidence for floods associated with the Laurentide ice sheet exists in such features as the Channelled Scablands (Bretz 1923), and in giant meltwater channels on Baffin Island (Kleman *et al*. 2001). However, the subglacial terrain of the former Laurentide ice sheet is relatively flat, and the presently existing large lakes (the Great Lakes, Great Bear Lake, Great Slave Lake) are the best candidates for large former subglacial lakes. It is likely that their floods would have been large scale, with outlet discharge direct to the North Atlantic or to the Gulf of Mexico. One reason for supposing that these lakes also existed subglacially is that their existence is closely related to the border between exposed crystalline basement rock and till-covered sedimentary rock (Clark & Walder 1994; figure 4). This dramatic association suggests that the transition between the two bed types was instrumental in excavating the lakes, and the obvious mechanism to do this would be associated with a transition from low basal water pressure upstream to high basal water pressure downstream. The consequence of this for basal sliding is to cause a change in basal stress at the bed, and consequent enhanced longitudinal stress. That such enhanced stresses should cause enhanced excavation seems plausible. In addition, the effect of such a basal water pressure jump would be to cause a damming effect because of the associated jump in the hydraulic head, which might be expected to lead to lake formation. Yet, very little is known about subglacial drainage under ice sheets, and there may be many surprises in store.

For the Laurentide ice sheet at the last glacial maximum (LGM), we have explored likely locations where water could pond to produce subglacial lakes. Using a digital elevation model of the current land and seafloor surface and elevation models of ice surface topography, we computed three-dimensional surfaces representing the hydraulic potential at the ice sheet bed, measured in pascals. We make the assumption that the ice sheet was wholly warm-based and that basal melting was uniformly spread. Meltwater should follow the maximum gradient of the hydraulic potential surface, and in locations of local minima in this surface, water should pond to accumulate as a subglacial lake. Meltwater routing algorithms within a GIS were used, and local minima in the hydraulic potential surface were ‘flooded’ until the lip was reached. The area of flooding at such locations becomes a prediction for subglacial lake position, and the depression of the local minima below the general surrounding surface we take as a guide to the likelihood of ponding in this location. Figure 8 displays predicted lakes for one experiment.

Given that lake prediction is sensitive to ice surface elevation and slope, and that we do not know this very well for a palaeo-ice sheet, we assessed the sensitivity of prediction by varying the ice surface. Three ice sheet geometries were used: the single-domed CLIMAP (1976) and Peltier's (2004) ICE-5G and our own multi-domed rendition approximating to Dyke & Prest (1987). These cover a range from single to multi-domed geometries and variation in thickness. Predictions were also computed for half-thickness versions of these ice sheets to assess sensitivity to lower-sloped and thinner ice sheets, a crude proxy for ice sheets during retreat. Three further experiments were run using different boundary conditions of effective pressure. The first assumed that effective pressure was uniform across the ice sheet bed with water pressure equalling ice overburden pressure. The second and third runs used varying effective pressures. For the area of the hard-bedded Canadian Shield, we presume subglacial drainage was predominantly by R-channels, and for the softer geologies beyond the Shield, drainage is presumed to be predominantly in canals within sediment (cf. Clark & Walder 1994). According to Walder & Fowler (1994), the effect of different drainage styles on effective pressure is to produce significantly higher effective pressures over the Shield. We conducted runs with a 1–20 and 1–50 bar difference between non-shield and shield areas, respectively. Lake predictions were made for each of the 12 different sets of boundary conditions as outlined above. Figure 8 is an example of one of these outputs. We note that we have not accounted for isostatic depression of the bed and the effect this would have on basal topography, and that predictions are based solely on hydraulic potential, so in no way does this allow for dynamic changes to drainage or the ice sheet surface, as drainage or indeed subglacial lakes evolve.

In general, the ICE-5G ice sheet produced the most abundant lakes, which is surprising given that it is the thickest of the ice sheets and hence has steeper surface slopes to drive basal water out. It seems to be its geometry that produces this effect, with the main ice dome in Keewatin, over Yellowknife. As expected, the half-thickness ice sheets produced more subglacial lakes due to the lower surface slopes, which implies lake occurrence should be more prolific during build-up and ice retreat, compared to the LGM. Introduction of different effective pressures for the Shield and non-Shield areas created a jump in hydraulic potential at the Shield boundary, which in places, where ice surface slopes were favourable, yielded many lake predictions. Such lakes are clearly evident in figure 8. It is apparent from the different runs that certain locations preferentially yielded subglacial lakes, and figure 9 is an attempt to represent this by presenting an amalgam of lake predictions from all 12 cases. In summary, lakes are predicted to occur in bed depressions that are deep enough in relation to ice surface slopes above them, and in locations whereby the hydraulic jump at the Shield boundary and its interaction with ice surface slope, traps water. Our representation of bed topography did not include the bathymetry of the major Canadian lakes systems (Great Bear, Slave, and the Great Lakes), and this is why they have not appeared as predicted subglacial lakes. However, as these are mostly deep lakes, it is sensible to presume that they also had a high potential for accommodating subglacial lakes. Lake Winnipeg might be an exception to this as it approximates a large puddle.

The suggestion that Laurentide lakes might have a short, millennial time scale for flooding has an intriguing possible consequence concerning Dansgaard–Oeschger events (Dansgaard *et al*. 1993; Bond *et al*. 1997, 1999), as they provide a perturbation in freshwater flux to the North Atlantic (Ganopolski & Rahmstorf 2001; Rahmstorf 2002), which is thought to be the cause of these events. In order for the theory to provide a suitable discharge and periodicity, we appear to need to have peak fluxes of the order of 0.1 Sv (10^{5} m^{3} s^{−1}), and that the flood periodicity be of the order of 1500 years. We see that the present flood theory can fit both these requirements (of order of magnitude).

Dansgaard–Oeschger events appear to have ceased at the end of the last ice age, although the 8200 year cooling event is thought to be associated with a massive discharge of the pro-glacial Lake Agassiz into the Hudson Strait (Clarke *et al*. 2004), but Bond *et al*. (1999) have suggested that pale imitations have continued into the Holocene. In particular, they suggested that the Little Ice Age between 1500 and 1800 might have been another event of the same type. If this was due to a perturbation in ocean circulation, then, following the logic of our thesis, we would suppose that such weaker climatic disturbances were driven either by floods from Antarctica, which may have rapid impact on Northern Hemisphere climate (Knorr & Lohmann 2003; Weaver *et al*. 2003; Richardson *et al*. 2005), or by floods from below Greenland. Lakes below the Greenland ice sheet have not been identified, but it is plausible that they exist, given inferred basal melt conditions (Fahnestock *et al*. 2001; Dahl-Jensen *et al*. 2003).

## Acknowledgments

We thank Duncan Wingham for several illuminating conversations. G.W.E. thanks W. H. R. Evatt for financial support. A.C.F. thanks the University of Limerick for their continued support of his research. C.D.C. thanks Monica Winsborrow for sharing her digital data of the Laurentide Ice Sheet.

## Footnotes

One contribution of 14 to a Discussion Meeting Issue ‘Evolution of the Antarctic Ice Sheet: new understanding and challenges’.

- © 2006 The Royal Society