Water Flow in Glaciers: Jökulhlaups, Tunnels and Veins

Abstract

The physics of water flow within and under glacier ice is examined with special reference to the periodic catastrophic outbursts of water (jökulhlaups) from the subglacial lake Grímsvötn, Vatnajökull, Iceland. The lake is sealed until it reaches a critical level which enables it to lift the glacier, helped by a hydrostatic cantilever effect. The differential equations for non-steady water flow in a subglacial tunnel are derived and applied to the 1972 Grímsvötn outburst. The discharge: time relation observed during the growth stage, and the abrupt ending of the flood, are both very well accounted for by a theory which is insensitive to the details of the subglacial tunnel system. The steady state, in which an intergranular vein or tunnel is simultaneously melted open by frictional heat and closed by plastic deformation, may be stable or unstable according to the conditions imposed at the ends. This explains why the flow of water in a vein does not normally increase unstably as in a jökulhlaup. An ice-dammed lake does not drain away through the vein system because the driving force on the vein-water is towards the lake rather than away from it.

Résumé

On examine la physique des écoulements d'eau à l'intérieur et sous les glaciers avec une attention spéciale pour les débâcles catastrophiques périodiques d'eau (jökulhlaups) dûes aux vidanges du lac sous-glaciaire de Grímsvötn, Vatnajökull en Islande. Le lac est obstrué jusqu'à ce qu'il atteigne un niveau critique qui le rende capable de soulever le glacier à l'aide d'un effet hydrostatique en encorbellement. Les équations différentielles pour l'écoulement de l'eau non-permanent par un chenal sous-glaciaire sont établies et appliquées à la vidange du Grimsvötn en 1972. Le débit: sa variation en fonction du temps pendant la période de crue, l'arrêt brusque du flot, sont deux phénomènes dont la théorie rend très bien compte, bien qu'elle soit insensible au détail du réseau hydrologique sous-glaciaire. L'état d'équilibre dans lequel une canalicule intergranulaire ou un tunnel sous-glaciaire est à la fois ouvert par la fusion née de la chaleur de frottement et fermé par la déformation plastique de la glace peut être stable ou instable selon les conditions imposées aux extrémités. Ceci explique pourquoi l'écoulement d'eau dans une canalicule ne croît normalement pas comme dans un jökulhlaup. Un lac fermé par de la glace ne se vide pas à travers le réseau de canalicules intergranulàircs parce que la force prédominante s'exerçant sur l'eau des canalicules la pousse vers le lac plutôt que vers l'extérieur.

Zusammenfassung

Die Physik der Wasserführung in und unter Gletschereis wird besonders im Hinblick auf die periodischen, katastrophalen Wasserausbrüche (Gletscherläufe) des subglazialen Grímsvötn-Sees am Vatnajökull in Island untersucht. Der See ist abgedämmt, bis er einen kritischen Stand erreicht, bei dem er den Gletscher mit Unterstützung eines hydrostatischen Ausleger-Effektes heben kann. Die Differentialgleichungen für eine nicht-stationäre Wasserströmung in einem subglazialen Tunnel werden hergeleitet und auf den Grimsvötn-Ausbruch von 1972 angewandt. Das Verhältnis von Ausfluss und Zeit, das während der Zunahmephase zu beobachten war, und das abrupte Ende der Flut lassen sich sehr gut mit einer Theorie erfassen, die unabhängig von der Detailstruktur des subglazialen Tunnelsystems ist. Der stationäre Zustand, bei dem eine intergranulare Ader oder Öffnung zugleich durch Reibungswärme aufgeschmolzen und durch plastische Verformung geschlossen wird, kann je nach den an den Enden herrschenden Bedingungen stabil oder instabil sein. Dies erklärt, warum die Wasserführung in einer Ader gewöhnlich nicht so instabil zunimmt wie bei einem Gletscherlauf. Ein eisgedämmter See fliesst durch das Adersystem deshalb nicht aus, weil die Druckkraft auf das Aderwasser eher gegen den See hin gerichtet ist als von ihm weg.

1. Introduction

The physics of the flow of water in glaciers contains a rich variety of interesting problems. The mechanics and the thermodynamics of water flow through the intercrystalline vein system (Nye and Frank, 1973), through a thin sheet at the bed of the glacier, and through a tunnel at the bottom of the glacier or elsewhere within the glacier, have much in common. At the same time these different flows can present apparent contradictions and delicate questions. For example, water flow through a subglacial tunnel can generate enough heat by friction to enlarge the tunnel in a catastrophic way; but water flow through the vein system, if it occurs, evidently does not normally do this. Why should this be so? During a catastrophic outburst from an ice-dammed lake does the water necessarily flow beneath the glacier within discrete channels or can it flood out into a basal sheet? When a catastrophic outburst is in full spate how can it end prematurely, before it has exhausted the supply of water? These are some of the questions we wish to discuss.

General mathematical treatments of the flow of water in glaciers have been given, in particular, byRöthlisberger (1972), Shreve (1972), and Wcertman (1972). This paper has a different emphasis because, although much of the treatment tries to be general, we shall continually refer to a special and unusual phenomenon: the periodic draining of Grimsvötn, the subglacial ice-dammed lake in the centre of the ice cap Vatnajökull in Iceland; see Björnsson ([1975]),Tómasson ([1975]) and the recent book by Thorarinsson (1974). The draining of the lake, which used to happen about every 10 years and now happens about every 5 years, is accompanied by an enormous outburst of water (a jökulhlaup) from beneath the ice edge to the sea, the outburst place being 50 km from the lake itself (Fig. 1). For example, the outburst of 1934 (Thorarinsson, 1953, 1957), which was a normal one for that period, carried at its climax 4 to 5 x 10⁴ m³ s^-1 of muddy water (about a quarter of the flow of the River Amazon, and more than that of the River Congo) bringing with it “icebergs as big as three-storeyed houses”. Almost the whole of the sandur, or outwash plain, some 1 000 km² in area, was flooded. But by the next morning the flow of the main outlet river, Skeiôarâ, had returned to normal. Any road along the south coast of Iceland had to cross the outwash plain and so it is not surprising that for many years none was built. Nevertheless, the construction of a road and bridges has been undertaken since the last outburst in 1972 (when the maximum How was 8 x 10³ m³ s^-I) in spite of the confident expectation of the next one in about 1977 or 1978. It is therefore practically important as well as scientifically interesting to understand the Grimsvötn phenomenon fully.

Fig. 1. Vatnajökull ice cap, showing Grímsvötn, the subglacial ice-dammed lake, and the presumed route (broken line) of the ice tunnel system, 50 km long, under the glacier Skeiðarárjökull. The road recently constructed across the outwash plain, Skeidaðarársandur, in the path of the outbursts is shown by the broken line between the glacier and the coast. Contour heights in metres. From Björmsson ([1975]).

There are many other examples of catastrophic outbursts from ice-dammed lakes and various different types are recognized; we may mention especially the accounts by Thorarinsson (1939) and Liestøl (1956) and the useful short summary in the book by Embleton and King (1968); for other references see the paper by Whalley (1971). We do not attempt a full treatment of the whole subject here, rather we try to analyse in some detail various aspects and puzzles of glacier hydraulics which are thrown into prominence by the j ökulhlaup phenomenon.

Section 2 outlines our interpretation of the Grimsvötn phenomenon and Section 3 deals with the physics of the lake seal in more detail. Section 4 sets up the general differential equations for non-steady flow in tunnels, and these are applied specifically to the 1972 Grimsvötn flood in Section 5 to explain the observed discharge hydrograph and the abrupt ending of the flood. Section 6 discusses the essential difference which can make a waterway behave stably, as in drainage through a vein, or unstably, as in a flood in a subglacial tunnel. Finally, we consider why Grimsvötn does not drain out through the vein system (Section 7) and under what conditions glacial lakes in general will drain in this way (Section 8).

2. An Outline Interpretation of the Grímsvötn Phenomenon

The essentials of the Grimsvötn phenomenon as interpreted by Helgi Björnsson and the author, building on the work of many previous writers, can be modelled as follows. In his separate paper, which is complementary to the present one, Björnsson ([1975]) has compiled the necessary detailed topographical information and has made numerical calculations based on the most accurate surveys.

A geothermal area at the base of the ice cap, shown schematically in Figure 2, causes melting and thereby creates a depression in the upper surface. Thus the run-off from the top surface generated by the ordinary processes of ablation no longer flows entirely towards the edges of the ice cap, as it would normally do, but part of it flows into the depression. Water accumulates in a lake (shown shaded), generated both by the surface run-off" and by geothermal melting. In Grimsvötn roughly 25% comes from run-off and 75% from geothermal melting. Grimsvötn is almost completely covered by a floating plug of ice some 220 m thick, which is essentially an ice shelf.

It is commonly supposed that the floods from ice-dammed lakes start by the lifting of the ice dam by flotation (Strom, 1938; Thorarinsson, 1939), and Thorarinsson (1953, 1957) suggested that this was the case with Grimsvötn—although, later on, incomplete topographic information led him to doubt whether flotation was possible in this instance (Thorarinsson, 1965). Let us examine the flotation hypothesis. Suppose the water level in the lake at a given time while it is filling is L (Fig. 2). We can construct the broken curves B, obtained by reflecting the ice-cap surface in L and magnifying its vertical scale by a factor ρ_i(ρ_w—ρ_i), where ρ_i and ρ_w are the densities of ice and water; the factor is about 10. To see the meaning of curves B, keep the top surface of the ice cap fixed and imagine the whole ice cap immersed in an ocean whose surface level is L. Curves B represent the bottom surface the ice mass would then need to have if it simply floated, like an iceberg, in isostatic balance in the water (the “iceberg” would need to have negative thickness outside points M and N).

Fig. 2. Schematic section through the ice cup. The lake (shaded) is seated along DE and FG.

We see that along DE and FG the ice is “grounded” if water existed outside D and G in hydrostatic communication with the lake its pressure would be more than enough to lift the ice. Thus the absence of any such hydrostatic connexion is vital to preserving the lake. In plan view the grounding takes place over an annulus, and it is this that seals the lake.

Figure 3 shows the critical part of the seal FG in magnified form and to scale, but with vertical exaggeration. The lake level now rises from L and so curve B also rises (II times faster).

Fig. 3. The critical part of the seat, to scale with vertical exaggeration 20, Full lines are the glacier bed and surface. Topography fromBjörnsson ([1975]) .

The critical condition is reached when B reaches B' SO that it just touches the bed curve at K; at this point the seal breaks. In three dimensions one has to visualize the inverted and magnified top surface moving upwards and eventually touching the bed surface; the point of contact of the two surfaces is the critical point κ, where the width of the sealing annulus is reduced to zero. (To visualize the topology consider the intersection of a horizontal plane with a curved surface near a saddle point, and then tilt both surfaces together.)

In fact the flood occurs with the lake level still some 20 m below the level that this model would predict. According toBjörnsson ([1975]) this discrepancy is unlikely to be due to errors in the survey, and in the next section we suggest an explanation based on the fact that the ice exerts slightly less pressure on the bed than would be given by the weight of a free column. (We may remark here that Glen's mechanism (1954) would require a lake level 20 m above rather than 20 m below the prediction of the simple floating model.)

As soon as water penetrates the seal one possibility would be simply a slow subglacial overflow from the lake. There would be no jökulhlaup; a steady outflow through the seal would just balance the inflow to the lake. But we shall see (Section 6) that such a condition, although a steady state, is inherently unstable. This is the essential instability that is responsible for the jökulhlaup.

Accordingly, a subglacial channel system is opened up unstably by the flowing water, which moves towards the edge of the ice cap, the preferred route (Fig. 1), 50 km long, running beneath one of the outlet glaciers, Skeiöararjökull. Does the water flow under the glacier in one or more confined channels or tunnels, or does it possibly flow in the form of a continuous sheet at the glacier bed? The water is seen to emerge from beneath the ice front of Skeiöararjökull in two major and several minor channels, and collapsed “cauldrons” visible in satellite photographs (Björnsson, [1975]) suggest that for 10-15 km below the seal the water is in a channel; but from an observational point of view the water could be equally well in a fairly wide sheet or in a tunnel for the greater part of its course. Of course, the fact that the glacier itself does not surge forward during the flood sets some limits to the type and extent of the water sheet that one could contemplate. Detailed calculation shows that the hypothesis of one or several tunnels gives excellent agreement with the observations, and therefore we favour this model.

In a tunnel the water pressure p is less than the pressure p _i that would be given simply by the weight of the overlying glacier and so it tends to close by plastic deformation of the ice, but at the same time the frictional heat from the turbulent flowing water tends to enlarge the tunnel by melting (Röthlisberger, 1972). The melting rate dominates and the flood grows catastrophically in the way envisaged by Liestol (1956).

Now we have to ask what stops it. It might be expected to continue at least until the level L of the lake fell to the rim J (Fig. 3). In fact the flood stops after the water level in the lake has fallen only 100 m and is still 230 m above J (Björnsson, [1975]). The seal reforms—in spite of frictional melting by the torrential flow of water—and the suggested explanation, which differs slightly from Björnsson's ([1975]), is this. As the lake level falls, so does the pressure in the tunnel. Although this causes faster plastic closing, it is overwhelmed at first by the melting rate. Later on, however, when the flood is very great, the lake level is falling very rapidly and so is the pressure in the tunnel. We show later that, if τ measures the time remaining before asymptotic catastrophe, the pressure difference tending to close the tunnel increases at least as fast as I + (τ_o/τ)³, τ_o being a constant, and the associated strain-rate at least as fast as {I +(τ₀/τ)³)³, which is an extremely strong time dependence, approaching T^-9 Thus, quite suddenly, the plastic closure rate catches up with the melting rate and the flood rapidly abates. The seal is reformed.

This completes our qualitative outline of how the jökulhlaup starts and how it ends. We now examine the various assertions in more detail.

3. Hydrostatic Theory. Filling of the Lake

The hydrostatics of the seal can be presented in several different ways. It is instructive to look at it now from a slightly different point of view.

Let z_b be the height of the glacier bottom (Fig. 4) and z_s be the height of the glacier surface, taking as datum the level of the eventual water exit at the snout of the glacier. The thickness of the glacier is given by h = z _s-z _b Let us consider the situation between jökulhlaups, while the lake is filling. At the glacier bed we may suppose there is a thin film of water at a pressure p, say. If s is a coordinate measured down-glacier parallel to the bed, the effective negative pressure gradient driving the water towards the snout is

(1)

say, where the potential φ is given by

(2)

If we take p to be the ice overburden pressure pigh,

(3)

Fig. 4. Schematic section to illustrate the hydrostatics of the seal.

Figure 4 shows schematically a graph of φ(s) obtained essentially by drawing the bottom profile z _b of the glacier and adding to it 0.91 of the glacier thickness (it is similar to Röthlis-berger's “water-equivalent line”). Point F on the bed is the uppermost point of the water film, where the glacier goes afloat. Here the pressure in the water sheet must equal the hydrostatic pressure from the lake. Thus the lake surface, at height z₀, passes through v, the point on the curve of φ vertically above F. The lake is now sealed (with φ > ρ_wgz ₀) between F and G. If we assumed that as the lake level rose the glacier ice did not move vertically at all, the curves to the right of F would all remain fixed, the intersection point v would move to the right and the seal at κ would be broken when v reached M, the maximum of the φ curve vertically above κ. This would correspond to the hydrostatic argument of the last section. But it needs correction because the ice between F and κ, if not melted sufficiently by the advancing lake water, will be subject to a buoyancy force which will bend it upwards and thus reduce the pressure at the bottom of the still-grounded ice. The reduction of pressure occurs because the submerged ice does not immediately attain isostatic balance and so acts on the grounded ice as a buoyant inverted cantilever. The result is that p in Equation (2) will be a little less than pigh, φ will be correspondingly less and M will be lower.

The magnitude of this effect may be crudely estimated by assuming the ice to be perfectly plastic with yield stress k. The maximum vertical shear traction the cantilever can exert on the grounded ice is k. Suppose the vertical shear stress thus set up drops to zero over a length L’ of the grounded ice. Then, considering unit width, the force exerted on the bed from the weight of the ice in the length L', which is pighL’, is diminished by kh, and thus the pressure on the bed will be pigh-k(h/L’). If the transition length L = ½h, which is quite reasonable, this provides a reduction in pressure of 2k, which with a yield stress of I bar is just equivalent to the required 20 m of water. (I bar is an appropriate value for k because the lake takes several years to fill and so the rate of deformation is correspondingly small.) So the effect of the buoyant cantilever in prising open the seal gives a good reason why the flood starts when the lake level is still 20 m below the level the simple lifting model would predict. Note that by the construction κ lies, in plan view, between the summit N of the ice surface and the summit J of the bed, as we saw in Section 2 (these points are of course not necessarily true summits but only maxima in the profile drawn). The condition at the critical point is evidently dφ/ds = o. Ignoring the correction and using Equation (3), with h = z_s—z_b , this gives

or, in terms of the slopes of the upper surface α and the lower surface β, roughly

α = —1/10β.

That is to say, above κ the backward slope of the upper surface is one-tenth of the forward slope of the bed, a relationship that is already evident by inspection of Figure 3.

Before the seal is broken, the thin water sheet at the bed of the glacier has a divide at κ. Up-stream from κ it flows uphill to the lake; down-stream from κ it flows towards the snout. We return to this point when discussing the role of the intergranular veins in Section 7.

4. Differential Equations for Time-Dependent Water Flow in Tunnels in Glaciers

In this section we derive the general differential equations which govern- the non-steady flow of water in a conduit within a glacier, the conduit being simultaneously melted by the heat of friction and contracted by plastic deformation in the way discussed by Röthlisberger (1972). The conduit, which is assumed to be completely filled with water, can be thought of as a normal internal glacier water-course or as a large tunnel associated with a jökulhlaup. The change of specific volume on melting, the work associated with the change of diameter by melting and plastic deformation, the dependence of melting point on pressure, and the matter of temperature gradients across and along the tunnel all make it necessary to be careful with the thermodynamics. Readers who prefer to skip the details may like to go straight to the final equations (16)-(20). Related treatments have been given byRöthlisberger (1972), who was primarily concerned with steady How, Shreve (1972),Mathews (1973), and Weertman (1972).

Consider an inclined water-filled conduit (Fig. 5) of cross-sectional area S, which is a function of position s and time t. ∂S/∂t is made up of two parts, one due to melting and one due to plastic, deformation which we denote (∂S/∂t) Thus

(4)

where m(s, t) is the mass melted per unit length per unit time, and v _i = ρ_i ^-1 is the specific volume of ice. If Q(s, t) denotes the volume of water flowing through a cross-section per unit time, the discharge, we may write the continuity equation for the water flow as

(5)

where v_w = ρ_w ^-1 is the specific volume of water. Eliminating ∂S/∂t between Equations (4) and (5) gives

(6)

where Δν = v_w — v_i is the change in specific volume on melting, For the rate of plastic deformation we shall write

(7)

Here p_i is the ice overburden pressure, p is the pressure in the water and K_o is a constant depending on the shape of the cross-section. If the cross-section were circular and we used the Glen flow law for ice, this equation would be exact and n would be about 3 (Nye, 1953). It is plausible to assume that even with a non-circular section the rate will be proportional to a high power of the pressure difference.

Fig. 5. The thermodynamic system.

To apply the first law of thermodynamics we consider the system consisting of a short length l of the conduit, the end boundaries to the system moving with the water; thus l is a function of s and t. As the conduit contracts by plastic deformation, work is done by the external pressure. In time δt some ice will be melted and it is therefore convenient to place the cylindrical (or, more strictly, conical) boundary of the system so that it just includes the ice which is about to be melted. Thus the cylindrical boundary is fixed relative to the solid, but it moves because of plastic deformation. The mass within the system is fixed. The rate of work done on the system at the cylindrical boundary is then

The rate of work by the water pressure on the up-stream boundary is pQ and the total rate of work by both up-stream and down-stream boundaries together is

The rate of work done by gravity is Qρ_wg_sl, where g_s is the j-component of gravity. Thus the total rate of work done on the system per unit length is

(8)

We wish to know the rate of increase in internal energy of the system. The flow being turbulent, nearly all the water will be well mixed and at a temperature θ_w (s, t), say, but at the ice wall there will be a thin boundary layer across which the temperature falls θ_i (s, t) the melting point of ice appropriate to the water pressure. During time δt the water within the system consists of two parts: the water already present at time t and the new water being formed by melting. The total change of internal energy δU is thus the sum of two terms:

(9)

Here σ is the specific heat capacity of water, d/dt is the derivative following the element, Δu is the difference of internal energy per unit mass between water and ice at the melting point, and the water is taken to be incompressible. The first term is the internal energy needed to raise the temperature of the already existing water; the second term is the internal energy needed to melt the new water and raise it to temperature θ_w (the boundary layer is assumed to remain the same, gaining water on one side and losing it on the other). Note that in the first term the volume Sl is left outside the differentiation because this part of the water remains fixed in volume. We are neglecting any contribution from kinetic energy.

If we neglect conduction of heat across the ends of our moving system, the heat input to the system is zero; so the first law of thermodynamics enables us to equate the rate of work given by Equation (8) with the rate of increase of internal energy given by Equation (9). In doing so we make use of the general thermodynamical relation

(10)

L being the latent heat per unit mass. Hence

and substitution for (∂S/t)_pi from Equation (6) reduces this to

(11)

Equation (11) has the physical interpretation that if the system is taken to be only the entering water volume Q without the addition made by melting, Qρ_wg_s is the work done by gravity, —Q ∂p/∂s is the work done by the water pressure at the ends, ρ_w ^σS (dθ _w/dt) is the heat output used to raise the temperature of the water and the right-hand side is the heat output used in melting. However, to derive Equation (11) on these lines raises questions about the work done against the lateral pressure and the shear forces—questions which are avoided by the slightly longer, but we hope clearer, treatment we have given.

There is a further relation between the pressure gradient, the flow Q and the cross-section S of the conduit. Following Röthlisberger (1972) and Mathews (1973) we use the empirical Gauckler-Manning formula (Williams, 1970) for the mean velocity ü = Q/S for turbulent flow, thus

where R is the hydraulic radius (cross-section divided by perimeter) and n’ is Manning's roughness coefficient. This may be rewritten in the form

(12)

with N = (S/R²)⅔ρ_wgn^’2 . Since S/R² depends only on shape and not on size, Equation (12) holds, with N constant, for any shape of cross-section, provided it remains geometrically similar (for example, a semi-circle) and of constant roughness as it changes in size. For the special case of a circular section N = (4 π)^⅔ ρ^wgn’².

Now consider the heat-transfer problem (Mathews, 1973). When water at temperature θ_W flows turbulently down a tube of circular section whose walls are at temperature θ_i the rate of heat flow to the walls is given by an empirical relation (McAdams, 1951, p. 168) between three dimensionless numbers, the Nusselt number (Nu), the Reynolds number (Re) and the Prandtl number (Pr):

(Nu) = 0.023(Re)^0.8(Pr)^0.4 (13)

Here

where h is a heat transfer coefficient, D is the diameter of the tube, κ is the thermal conductivity of the fluid, ü is the mean velocity of flow, and η is the viscosity of the fluid. Equation (13) fits observations for (Re) between 10⁴ and 1.2 x 10⁵ and (Pr) between 0.7 and 100, with L₀/D > 60, L_o being the length of the tube. The rate of heat flow to the walls per unit length of tube is, in terms of h,

and in our case this is also

The Prandtl number is a material constant which for water is 13.5. Hence we find

(14)

with

(15)

The equations for non-steady flow in a conduit may now be collected together: Geometry and flow of ice:

(16)

Constinuity:

(17)

Flow of water:

(18)

Energy:

(19)

Heat transfer:

(20)

In the five equations (16) to (20) the course of the conduit, N and the ice overburden pressure p_i (s) (or p_i (s,t)), are taken as given, and so there are five unknowns Q, S, p, m and, θ_w which are all functions of s. and t. (θ _i is a known function of p.) Note that only the heat-transfer equation depends upon the assumption of a circular cross-section.

5. Application of the Equations to the Gímsvötn Jökulhlaup

(a) The main part of the flood

We make the hypothesis that the water in the Grimsvötn jökulhlaup flows in a tunnel or tunnels rather than as a sheet, and we look for a self-consistent solution of the equations corresponding to this case. If a tunnel is at the bed of the glacier, and we are not compelled to assume that it is, a crucial question will be whether in the putative solution p _i everywhere. For if this condition were violated the solution would not be valid; the water would Hood outwards from the tunnel, lifting the ice and forming a sheet.

First we simplify the general differential equations. It is found, not surprisingly, that the temperature difference θ_w —θ_i is a few degrees at most and so we neglect mσ(θ_w— θ_I ) in comparison with mL in Equations (19) and (20).

The role of advection in the problem can be confusing; the following explanation is based on calculations with Equations (19) and (20). If the tunnel were simply a pipe with perfectly insulating walls, the rise in temperature of the water as it lost potential energy in falling a distance z₀ from the lake to the exit would be gz_o/σ = 30 deg. On the other hand, if we allow cooling by the ice walls, assume no generation of heat within the tunnel, and imagine the tunnel to be fed with hot water at the top end, the water temperature will fall exponentially. At the high Reynolds numbers during the flood the distance for the temperature excess to decrease by a factor e is rather large: 14 km at the start of the flood and 90 km at the end (the total tunnel length is 50 km). From this point of view the transverse heat transfer is quite inefficient. But, of course, there is generation of heal within the tunnel, by friction, and it would be possible, in principle, for the water to get either hotter or colder as it travelled, or to remain at the same temperature. To take this last case, if the tunnel were fed with water at 0.9°C at the start of the flood (when (Re) = 1.4 χ 10⁷) or at 54^°C at the peak (when (Re) = 1.5 X 10⁸), there would be no gradient dθ_w/dt; the friction heat would just supply the amount of heat that can be transferred to the walls and used for melting, with none left over to heat the water any further. Of course, the water temperature θ_W has to exceed θ_i for the necessary heat transfer to the walls to take place, but 0.°C to 5.4°C is enough for the purpose (from this point of view the heat transfer is quite efficient). The outflow temperature during the 1954 flood was measured byRist (1955) as 0.05°C, but there is no information about the lake temperature. Thus it is not unreasonable to try taking dθ_w/dt = 0 in Equation (19). We shall do so because it leads lo excellent agreement with the other observations.

dθ_w/dt = 0 means there is no advection term in Equation (19), but it would be wrong to interpret it to mean that frictional heat generated at a given place is deposited locally, for we have seen that the characteristic travel distance is 14 to 90 km. (Mathews’ (1973) theoretical treatment of the advection effect, discussed qualitatively byLiestøl (1956), is the same as ours, but with the difference that in his case the water enters at 0°C and heats up by friction towards an asymptotic value. He points out that this case presents a problem because near the tunnel entrance there is then no thermal gradient to transfer heat to melt the walls. Gilbert's fuller solution (1971) allows the water to enter above 0°C.)

The four equations (16) to (19) now reduce to:

Geometry and flow of ice:

(21)

Continuitv:

(22)

Flow of water:

(23)

Energy:

(24)

for the four unknowns Q, S, p and m. Equation (20) merely serves to calculate θ_w — θ_i , if required. We recall that the heat-transfer formula (13) has only been tested up to Reynolds numbers of about 10⁵, while the Reynolds number during the flood is in the range 10⁷ to 10⁸; moreover, the formula is only for a circular section. Thus it is fortunate that Equation (20) is essentially independent of the other equations; it means that our results do not depend on the detailed applicability of Equation (13) to our problem.

The quantity (p_wg_s-∂_p∂_s) that appears in both Equations (23) and (24) is —∂φ/∂s, namely the negative gradient of potential φ defined as

(25)

where zt is the height of the tunnel at given s and t. the tunnel is at the bed of the glacier this is the same definition as Equation (2).

To obtain the solution appropriate to the growth of the jökulhlaup we now make the crucial assumption that melting of the tunnel is so fast that the rate of plastic closing is negligible in comparison, the justification being that this assumption leads to excellent agreement with the observed Q: t relation. Thus we neglect the last term in Equation (21). It follows from Equations (21) and (22) that

(26)

This means that, if there were no change of volume on melting, Q would be uniform down the tunnel at any given time; the melt water would not add anything to the flow because its volume would be needed to fill up the extra space made in the tunnel. In fact we see that the contraction associated with melting leads to a slight decrease in Q with S (later results show that this is less than 1% down the length of the tunnel).

Eliminating m and Q between Equations (21), (23) and (24) we have

(27)

where

(28)

In the expression for KI both N which depends on tunnel roughness and shape, and the potential gradient ∂ φ/∂s are functions of s and t. But the space average value of — ∂φ/∂s being simply the difference of potential ∂φ between the two ends of the tunnel, is p_wgz_o/l_o , where z₀ is the height of the lake level above the tunnel outlet and l_o is the length of the tunnel. z_o changes by only 8% during the flood and so we shall treat K_I as a constant, constructed from average values of N and - ∂φ/∂s Then Equation (27) integrates to give

(29)

where t = o is chosen as the time when S becomes infinite. S here is to be regarded as the average cross-section down the length of the tunnel.

The corresponding differential equation for the time development of Q is obtained from Equation (27) by using the Q: S relation in Equation (23), as follows:

(30)

where

(31)

Again, we have treated N and ∂φ/∂s as constants and we are treating Q as uniform with s. Equation (30) integrates to give

(32)

The most accurate discharge hydrograph (Q : t relation) made for a jökulhlaup from Grimsvötn is the one by Sigurjón Rist for the 1972 jökulhlaup (Rist, [1974], p. 57) which is reproduced here as Figure 6. It shows the separate discharges from the three main outlet rivers. The three discharges are summed in Figure 7 to give observed values of In Q: t which are compared with the predictions of Equation (32) with the parameters chosen to give the best fit (K ₂ = 7.45X10^-7m^-3/4s-3/4 with the asymptote t = o at 11.00 h on 30 March). The flood falls abruptly 6.5 d before Q would be theoretically infinite.

Fig. 6. Hydrograph for the 1972 jökulhlaup. The discharges of the three main outlet rivers Skeiàarà, Gigja and Sula nre shown separately by the broken curves. Measurements by Sigurjon Rist (reproduced from fig. 3 of Rist ([1974]) by permission of the editors of Jökull).

Fig. 7. The sum of measured discharges in Figure 6 (open circles and broken curve) plotted logarithmically and compared with Equation (32) [foil curve). The close fit of the rising limb with theory is remarkable.

The excellence of the fit between observation and theory up to the flood peak is rather astonishing. Although Rist's estimated error in the observations of Q is ±20% they fit the theoretical curve to within about 4%. The plotted points were obtained (by Björnsson) from his summation of Rist's graphs and they show no observable systematic deviation from the theoretical curve only random deviations. I suggest that, while the method of observation (from a helicopter, sighting the water level in the rivers against fixed marks and measuring the geometries of the channels after the flood was over) may have given a 20% systematic error, the random errors are much smaller. In any case the agreement with the predicted analytical form gives confidence in the various simplifications that have been made in the theory. Since the index 5/4 in Equation (30) is close to 1 the Q: t relation is not far from exponential, but nevertheless very distinguishably different, as can be seen by the departure from a straight line in Figure 7. In fact a plot of In Q: In (—t) gives an exponent in Equation (32) of —4-oo±o.o6. The constant N calculated from Equation (31) is 709 m⁸/³ kg, which for a single tunnel of circular section corresponds to a Manning roughness coefficient n = 0.12 m^-1/3s. According toRöthlisberger (1972, p. 181) n’ could be expected to range from about 0.01 m^-1/3 s for a straight smooth pipe to 0.1 m^-1/3 s for a meandering boulder-strewn torrent at the bed. Our calculated n’ lies at the rough end of the range; this is to be expected if there are several tunnels operating at once (as is observed at the outlet), and if their sections are non-circular. Moreover we have ignored the effect of the transport of the large amount of sediment and boulders.

After the main peak on 23/24 March there is a subsidiary peak which is associated with water having a different deuterium content. This suggests that, after the main outburst from Grimsvötn, some water was released from a different storage place (Björnsson, personal communication). We note that, provided this water used the same tunnels as the main flood from Grimsvötn, it makes no difference to the potential gradient or to N, and so the rising part of the Q: t curve is quite unaffected by it. The discharge up to the flood peak is determined by the potential gradient and by the growing tunnels, not by intermediate storage reservoirs.

K_I can be estimated from Equation (28) as 4.91X10-^-2m^-2/3s-1, and Equation (29) then gives the plot of In S: tshown in Figure 8.

Fig. 8. Theoretical cross-section S verms time, corresponding to the theoretical Q : t curve of Figure 7. Although the main part of the theory does not assume a single circular cross-section, the scale on the right shows what its diameter D would be.

We have seen that Q during the growth stage is very nearly uniform with s, and if N were uniform we could deduce that S and dø/ds were also uniform. φ would be linear in s. But non-uniformity in N, which is affected both by roughness and by the shape of the cross-section, means that we cannot treat ø as exactly linear in s. (This is a pity, because if ø were strictly linear in s we could go on to deduce all the details p(s, t).) Figure 9 (to scale) shows a possible curve, labelled ø_I, for ø(s) near the start of the flood (but after the melting rate has dominated plastic contraction). The top end is fixed by the lake level L_I above the position K the seal. Curve B in Figure 9 represents pwgzD^(s) and may be thought of as the bed of the glacier with the height scale translated into a scale of water pressure. If we assume that the tunnel, at any rate near the top end, is at the bed of the glacier (z_t - Z_b), the water pressure p(s) is simply the difference between curves φ_I and B. We want to compare it with the ice overburden pressure pi(s), which we shall assume is pigh. Accordingly we plot pigh upwards from the bed curve B as a base, to give curve G’ (Röthlisberger's “water-equivalent line”), the actual glacier surface being G. (Before the seal was broken the potential φ followed curve G.) If G' is above φ_I, the ice pressure pi exceeds the water pressure p and the tunnel tends to close plastically (although of course it is enlarging very much faster by melting). But if G’ chanced to be below ø_I the water pressure would exceed the ice pressure, the ice would lift from the bed, the water would spread out into a sheet and the tunnel model would have broken down. At the other extreme, if the bed curve B happened to lie above ø_I at some place, the water pressure would be negative relative to atmospheric pressure. At a negative pressure of about I bar cavities of water vapour would form.

We have drawn ø_I as shown in Figure 9 hearing in mind that, with uniform N ø_I would be linear with s and that the end points are known. ø_I cannot He too far below G’; for example, if ø_I were drawn 20 bars below G’ in the lower part of the glacier, pi—p would be unacceptably high; the melting rate would not be enough to overcome plastic closing.

As the lake level falls curves G’ and B remain fixed, by definition, but ø(s) must fall. ø ₂ shows a possible curve for ø(s) at the end of the flood. The top end is fixed by the new lake level L₂ and we have simply pivoted the curve about the end point D. The construction in Figure 9 shows that the value obtained for pi—p is very insensitive to the position of the bed. Or, put in another way, the value obtained for pi—p depends very little on whether the tunnel is assumed to be at the bed of the glacier or at some higher position.

(b) How does the jökulklaup stop?

According to Björnsson's calculations the flood abates when the lake level L₂ is still some 230 m above the level of the saddle point J (Figs 3and 9). At first sight this is surprising. Once the flood has started why does the lake not drain oui as far as the height of its rim will allow it? Björnsson notes that the end of the flood corresponds roughly to the moment when the lower surface of the plug of ice, 220 m thick, that floats on the lake falls to the level of J. I believe that this is a coincidence and that the flood would stop punctually as observed even if there were no ice plug on the lake. The reason is as follows.

Fig. 9. Longitudinal section of Skei∂arárjökull during the jökulklaup, to scale with vertical exaggeration x 20. The vertical scale is height translated into water pressure, B, glacier bed; a, glacier surface; G,’ water equivalent line; ø₁, ø₂ potentials at start and end of flood; L₁,L₂ lake levels at start and end of flood. Topography from Björnsson ([19751])

Although in our approximate calculation we have neglected the plastic contraction of the tunnel in comparison with the enlargement by melting, we can now estimate pi—p, for this is the difference between the fixed G’ curve and the changing curve for ø in Figure 9. During the flood the lake level falls 100 m. If we accept the lines drawn for ø₁ and ø₂ in Figure 9 the resulting increase in (pi—p) near the top of the tunnel would be of the order of 20 bars. This would represent an increase by a factor of about five. We cannot be precise about the magnitude of this factor because there are several possible complications that affect the course of the ø curves; therefore we do not insist on the fine details of Figure 9; the essential point is that the drop in lake level causes a very significant increase in the pressure difference that tends to close the tunnel. The effect is especially marked near the top of the tunnel, and it is worth noting that the collapsed cauldrons, reported by Björnsson as occurring 10 15 km below the seal, are in just this region.

If we make the very conservative assumption that the pressure difference trebles, the plastic contraction, which is proportional to the cube of the stress, increases by a factor of 27 and is now more than enough to overcome the enlargement by melting. When the two effects just balance, the flood is at its maximum. Then the tunnel contracts plastically and the flood comes to an end.

It is crucially important to notice that the effect of plastic contraction does not simply increase gradually with time as the lake level falls; on the contrary it grows extremely suddenly at the peak of the flood. There are two reasons for this. First, the fall in lake level is far from steady, If we write τ = -t, so that τ is the time remaining before the asymptote where Q = ∞ would be reached, we have seen that Q Q∝ ^-4 This means that, even if the lake had vertical sides so that the rate of fall of ils level was proportional to Q, its level would be given by

Z_o + z∞ — E_oτ ^-3 ,

where E_o is a constant and z∞ is the starting value of z_o at T = ∞. Second, because pi—p varies linearly with the lake level (Figure 9), it must be of the form

(33)

where P∞ and τ_o are constants. (It is preferable to keep a constant term in this expression, rather than having pi—p extrapolate to zero as T ∞, because, of course, the relation does not hold during the earliest part of the flood when plastic contraction may still be comparable with the melting rate.) Therefore the plastic deformation rate is proportional to

(34)

This time dependence means that the plastic deformation rate remains fairly constant and then increases very fast around τ = T ₀, becoming effectively like τ ^-9'. The fact that the lake has sloping sides makes the increase even faster.

This very abrupt growth in the rate of plastic contraction has two consequences. On the one hand it makes it legitimate to neglect plastic contraction during the Hood and still get excellent agreement with the observed Q : t relation. On the other hand, it explains the sudden termination of the Q ∝ T ^-4 relation without an observable transition period.

Figure 10, which is partly calculated and partly schematic, shows strain-rate S ^-I [dS/dt] versus time. M is the opening rate due to melting, while P is the contraction raie due to plastic deformation. During the main part of the flood, plastic contraction being small, M may be calculated from Equation (29); this gives the surprisingly simple relation S ^-1[dS/dt] = 3T^-I independent of K ₁, shown by the full-drawn part of M. Curve p is calculated from Equation (34) by making it intersect the calculated M at the observed time of flood maximum, because here the two rates must very nearly balance. The further assumption that (pi—p) trebles during the flood fixes curve p completely. If we had assumed, instead, a fivefold or even greater increase in (pi-p) and had taken into account the slope of the lake sides, curve p would start at a lower level, would remain low for longer and would turn upwards even more suddenly to cut M at the same point.

Fig. 10. The difference between the rate of opening by melting, M, and rate of closing by plastic deformation, P, during the 1972 jökulhlaup. Full curves are calculated on the assumption that P ≪ M broken curves are diagrammatic only.

The remaining broken parts of curves M and p are diagrammatic only. We may imagine an ordinary sub-glacial drainage tunnel existing up to point w, with M and p in balance. At this time the lake rises to the critical level and M begins to exceed p. A plausible extrapolation of the calculated part of M places w on 1 March, which is the day the local people sensed a sulphurous smell at the outlet. On 13 March they were in no doubt that a jökulhlaup had started, M was now much larger than P, and it remained so until very shortly before the flood maximum at midnight on 23/24 March. The details of the curves of M and P near their intersection are not certain, because we have no solution when P is comparable with M, and we must remember the spatial variations, but we could argue as follows. At the intersection point where M = P, dS/dt = o and S is a maximum. dø/ds is decreasing and so Equation (23) suggests that Q also decreasing. Equation (24) then implies that m is decreasing. Thus the maximum of the curve M in Figure 10 occurs before the intersection point, as indicated by the broken line (the full part being calculated on the assumption that P ≪ M) Eventually both M and P fall back to their original values.

To summarize, our theoretical derivation of the Q: t relation in Figure 7 during the growth of the flood has depended essentially on the conservation of energy, and on the Gauckler-Manning formula in the form of Equation (12) with N constant. Thus the argument has not depended on assuming any particular shape of cross-section or on whether the water is in one or in several different tunnels, but it has assumed that the cross-section or cross-sections remain geometrically similar as they grow. It does not depend on there being other storage places besides Grimsvöln, provided they drain through the same tunnel system. The argument for the sharp cut-off of the flood by plastic contraction of the tunnel is equally independent of the cross-sectional shape, or of whether there is just one tunnel or several; for it simply uses the fact that the plastic contraction rate is proportional to a high power of (pi—p). This insensitivity to details makes the close agreement of theory with observation, which is startling at first sight, rather less surprising.

6. Intergranular Veins and Jökulhlaups

If the intergranular vein system in a glacier is a connected network as proposed by Nye and Frank (1973), water, being the denser medium, will flow downwards. Viscous resistance to the flow will generate heat which will melt and tend to enlarge the veins. The realization of this led Lliboutry (1971) to ask why, in that case, the whole glacier does not melt itself away in a rather short time. He suggested that the answer is that the vein system is blocked. Another answer could be that it is the limited supply of water to the vein system that prevents catastrophe. Vein flow, assuming that it does occur, and jökulhlaup flow clearly have something in common. But what is the essential difference?

As Röthlisberger has shown (1972), when water flows through a tunnel, the rate of opening by frictional melting can be balanced by the rate of closing by plastic deformation. It is also possible that this is essentially what happens with the vein system in a glacier (Shreve, 1972), so let us trace the consequences of such a supposition. When the “tunnel” is a vein, its cross-section has three corners of molecular sharpness; it might be held (Nye and Frank, 1973) that these would produce virtually infinite stress concern rat ions and thus make plastic contraction possible with merely an infinitesimal underpressure. This is incorrect. Dr R. G. Oakberg (to be published) has now calculated the stress distribution that would be produced at a vein if the ice were isotropic and obeyed the Glen flow law. The stress near the vein corners is biaxially hydrostatic with no infinity, as conjectured by Nye and Mac (1972, p. 98). This means that, for its response to externally applied stress, the triangular vein may be replaced by an equivalent tube of circular cross-section.

The differential equations for the time-dependent water (low in a vein may be written down in close analogy with those for a tunnel, Equations (16)-(20). We neglect the small pressure difference across the curved vein surface and also the effect of surface energy. We assume Poiseuille flow (with the same equivalent circular cross-section) and perfect heat transfer so that ø_w = ø_iUsing the pressure: melting-point relation dø_i/dp = -C, where G is a positive constant, means that the convective term in the energy Equation (19) is

Thus Equations (16)—(19) are replaced by:

Geometry and flow of ice:

(35)

Continuity:

(36)

Poiseuille flow of water:

(37)

Energy:

(38)

y = pw/σC = 0.313 is a dimensionless constant whose appearance in the energy equation represents an effect noted by Röthlisberger (1972) in connexion with flow in ice tunnels. In Equation (35) we have used the form from Nye (1953) appropriate to a circular cross-section, A* being a constant in the Glen flow law.

With this model we can first ask how much underpressure in the vein is needed to provide enough plastic contraction lo compensate for the frictional melting. Putting ∂/∂t = o in Equation (35), substituting from Equation (38) for m and from Equation (37) for Q gives

(39)

We choose a vertical vein (g_s = g) and estimate ∂ p/∂s = pig; then

(40)

For a vein of circular section 20 μιη in diameter this gives a pressure difference between ice and water (pi-p) of 0.6 bar (6 m head of water). Thus the underpressure in a vein relative to the surrounding ice is not large and for many purposes we can neglect it. But if we neglect it altogether there is nothing to prevent the vein enlarging by melting and, if had access to an unlimited supply of water, growing indefinitely large.

So plastic contraction of a vein to balance the frictional melting makes possible a steady state of the Röthlisberger type. But is the steady state stable against perturbations? The answer is yes and no. If the pressure difference between the ends of the vein is kept fixed, the steady state is unstable (as plausibly asserted byWhalley (1971, p. 172) for Röthlisberger channels). But if the vein is fed from a reservoir it may be stable, as we shall see in a moment.

First we show the instability for fixed pressure by considering the thought-experiment of Figure IIa. A reservoir A, supplied with water from an inlet pipe c, is connected by a short length of ice tunnel with circular cross-section (the model vein) to another reservoir B, the tunnel being subjected to a fixed ice overburden pressure pi The slope of the tunnel is arbitrary. The object is to set up a Röthlisberger-type steady state. We do this by providing both A and B with very large capacity overflows at fixed heights, and we always maintain sufficient flow through C to ensure that the water levels remain at three heights. Thus the pressures pA and pB at the two ends of the tunnel are fixed. There will now be only one cross-sectional area S of the tunnel that will give steady conditions (pressure differences and gradients being fixed, the balance Equation (35) decides m\S, the energy Equation (38) decides Q/S, and then the Poiseuille flow Equation (37) decides S). Let the cross-section be chosen accordingly, and call it S ₀. The state is now steady.

Fig. 11. Thought-experiments to test stability of a tunnel or vein which is simultaneously opening by melting and closing by plastic deformation. At the outflow end the pressure is fixed. Arrangement (a) with the inflow pressure fixed is unstable; arrangement (b) with the inflow rate fixed is stable, provided reservoir A has an area less than a critical value.

In the general non-steady state, when pressures and conduit gradients are held fixed, Equation (37) implies Q ∝ S^z and hence, by Equation (38), m oc S¹. Thus in Equation (35) the contribution to ∂S/∂t from melting is proportional to S² while the negative contribution to ∂S/∂t from plastic closing is proportional to S. The steady state occurs at the intersection point R (Fig. 12) where the contributions balance. If S > S ₀, the melting rate exceeds the plastic closing rate and so the 'cross-section of the tunnel increases unstably (less water overflowing from A). Similarly if S < SS_o the cross-section diminishes unstably to zero (more water overflowing from A)

Fig. 12. [dS/dt] versus S far the unstable arrangement of Figure 11a.

We can achieve stability by changing the system to Figure 11b. This arrangement is identical to that in Figure 11a except that there is no overflow pipe for reservoir A. We start with the tunnel having cross-section S₀ and with the levels exactly as before—but in order to maintain these levels we have to reduce the inflow from c by the amount that formerly overflowed from A when we had the steady stale. We keep this inflow fixed. Then we have set up exactly the same steady state in the tunnel as before.

When the steady state is perturbed the level of water in A can change, and so to examine stability we must consider the combined system of reservoir and tunnel together. It is shown in Appendix A that for small disturbances the stability depends on the surface area of the reservoir A: below a certain critical area the system is stable; above this area it is unstable. (H. Röthlisberger had already surmised this behaviour and his letter to me on the question led to the calculation given in Appendix A.) The system illustrated in Figure IIa, which we have already seen to be unstable, is equivalent to having an infinite reservoir. The stability condition does not depend on there being a fixed inflow into reservoir A. The inflow represents a forcing term which merely serves to determine what steady state one has.

The differential equation for small perturbations of cross-section of the tunnel is second-order and linear. Mathematically, instability arises because the damping coefficient can be negative. At the critical reservoir area the damping is zero, and so at and around this area the system oscillates, stably or unstably. For small reservoir areas the system is over-damped and there is no oscillation. Similarly for large areas a disturbance grows unstably without oscillation. The critical area of the reservoir is very much greater than the cross-section of the tunnel or vein in all practical cases.

The analysis described is for a length of tunnel that is short enough for us to be able to regard the pressure gradient, discharge, cross-section, and the other variables as changing only in time and not in space. A study of the behaviour of spatial variations would be interesting but will not be attempted in this paper;Nye and Frank (1973) suggest that kinematic waves will travel along a vein. However, the present analysis goes far enough to show that the stability of a tunnel or a vein in a Röthlisberger-type balance depends on the end conditions. The pressure being held fixed at the outflow end, the system is unstable when the inflow is held at fixed pressure, but stable when fed from a reservoir of small area. The latter combination is normally appropriate to a vein in a glacier. I suggest that the stability of this combined system of vein and small reservoir explains not only why veins do not normally grow catastro-phically in size, but also why they generate just the right underpressure to balance melting by Plastice contraction. This ceases to be remarkable as soon as one realizes that it is merely a consequence of the tendency to seek a stable steady state. On the other hand, feeding a tunnel from a reservoir at fixed pressure is an appropriate model for draining a large lake and the instability of this configuration is responsible for jökulhlaups. If it were not for this instability, the lake level in Grímsvötn would remain at the critical level for lifting, and the lake would simply overflow gently through a subglacial drainage system, its outflow being balanced by its inflow.

When a tunnel is growing unstably (Fig. 10) there are two stages: at first the melting rate is still comparable in magnitude with the plastic closure rate, although of course it is faster. In the second stage the melting rate is considerably faster. In Section 5a it was the second stage that we analysed in detail, because we treated the plastic closure rate as negligibly small. We have not given a detailed mathematical analysis of the first stage (which would include all the details of the variations with s), but we have identified in this section the essential feature that causes instability.

7 Veins and Ggrímsvötn

In view of the argument just given why does Grimsvötn not drain out through the inter-granular veins long before the critical level for floating is reached? After all, it is essentially a reservoir at fixed pressure.

We have already noted (Section 3) that before the seal at κ (Fig. 3) is broken a thin water sheet at the bed of the glacier up-stream from κ is driven uphill towards the lake, not away from it. It is simple to extend this argument to include the intergranular vein water as follows.

FollowingShreve (1972) define the potential ø = pwg (z – z_o) + p,

(41)

where p is the pressure in the veins and z₀ is the height of the top surface of the lake. The effective pressure gradient driving the water is then — grad ø This definition of ø is the same as Equation (2), except that it now applies to veins in the interior of the glacier rather than only to the bed, and the zero level for potential has been changed.

Neglecting the underpressure in the veins, take p as the overburden pressure of the ice at the point (x, z) in question, pig(z_s — z). Then

(42)

The equipotential ø = o is

This is the curve ANC in Figure 13a, obtained by reflecting the surface of the glacier in the lake level and exaggerating the scale by a factor of 10; it is the same as curve B in Figures 2 and 3 and represents the bottom surface the ice would need to have to float in isostatic equilibrium. AF is the hypothetical bottom surface of the ice in contact with the lake. The other equipotentials in Figure 13a are obtained by simply translating the curve ø = o up and down. If the veins made the ice isotropically permeable the lines of water flow would be orthogonal to the equipotentials, but even taking account of anisotropy the only result we need is that the vein water will flow from higher to lower ø When the glacier bed is as shown, before the seal is broken, clearly the vein water is forced towards the lake.

Fig. 13. (a) Equipotentials for vein-water within the ice near an ice-dammed lake such us Grimsvötn. (b) Schematic equipotentials for vein-water when an ice-dammed lake dos not reach bedrock.

The effect of the underpressure in the veins is to reduce p, and therefore ø throughout the ice, thus lowering the potential barrier; the reduction could correspond to a few metres of water head.

8. Veins and A Lake Perched on the Ice

The foregoing analysis is readily extended to the situation of a glacier lake that does not necessarily reach bedrock. Figure 13b shows an idealized section through such a lake which has a horizontal ice bottom at z_t. Because a point (x, z) may now have both ice and water above it, the vein pressure/; in Equation (41) is to be taken as the overburden pressure due to the overlying ice and any overlying water that may be present. The equipotential ø = o and those above it are the same as before. To draw the equipotentials below it we need only observe that ∂φ/∂z = (pw – pi)g in ice and ∂φ/∂z = oin water.

If water is to escape from the lake by draining through the veins, the glacier bed must be everywhere below the equipotential φ = o because otherwise with an impermeable bed there would be no path that did not cross the forbidden region ANC of higher potential, N is a critical point. In terms of ice freeboard H, if the bed is less than 10Hbelow the lake level the lake will not drain through the veins. But if the bed lies below this critical depth the lake may drain through the veins, making use of the instability appropriate to a large reservoir of fixed pressure.

It appears that this is why glacier lakes, particularly marginal lakes where the glacier depth is small, are able to exist at all in spite of the porosity provided by the vein system in the ice; and this is also why lakes in the middle of temperate glaciers, where the ice is deep, are seldom if ever seen. It is interesting to compare this explanation for marginal lakes with that given byShreve (1972, p. 212); as far as I can see the explanations do not conflict, although they differ in emphasis.

A lake sealed in this way can empty suddenly if forward motion of the glacier establishes a path, say a crack or crevasse, that short-circuits across the forbidden region ANC. But we must remember that the region is only forbidden because veins have time to adjust their internal pressure to be about equal to the overburden pressure—and so a pathway must be opened quickly if it is to be effective. This, rather than a change of water level, may sometimes explain the sudden emptying often observed in marginal glacier lakes.

Acknowledgements

Helgi Björnsson visited Bristol University Physics Department for two years with the. aid of a European Council Scholarship (administered by the British Council). This paper has been able to focus on the classic jökulhlaup from Grimsvotn largely because of the detailed information on the phenomenon that he has assembled (Björnsson, [1975])- I am most grateful to him for his congenial collaboration.

I should also like to mention the field meeting of the International Glaciological Society in Iceland in 1970, organized by Professor Sigurdur Thorarinsson, which acted as a springboard for the work. Professor C. F. Raymond and Professor H. Röthlisberger both commented very helpfully upon the manuscript.