Fowler A. Mathematical Geoscience

Подождите немного. Документ загружается.

Chapter 11

Jökulhlaups

In October 1996, a volcanic eruption underneath the ice cap Vatnajökull in Iceland

became an international news story. The eruption itself was spectacular enough,

sending clouds of ash high into the atmosphere, but its more impressive feature was

that the ﬁssure eruption took place under about 500 metres of ice. The erupted lava

melted the basal ice, causing ﬁrst a huge sag in the ice surface, and eventually its

complete collapse and the formation of an ice canyon, hundreds of metres deep,

within the ice.

It became evident that the subglacial meltwater was ﬂowing from the ﬁssure

towards the subglacial lake Grímsvötn, which lies in the caldera of a volcano under

the ice cap. As the weeks passed, the lake level rose (as could be inferred from the

uplift of the overlying ice), and eventually it reached ﬂotation level; that is to say,

the hydrostatic pressure of the lake water became sufﬁcient to lift the overlying ice,

and a jökulhlaup occurred.

Jökulhlaup is an Icelandic word, meaning literally a ‘glacier-burst’ (jökul-hlaup),

and it refers to catastrophic outburst ﬂoods from glaciers which occur in various

parts of the world, not only Iceland, but also Canada and the Himalayas, for ex-

ample. The one which occurred in 1996 was spectacular: icequakes at Grímsvötn

indicated a breach of the seal formed by the ice at the caldera rim (see Fig. 11.1),

and ten hours later a huge ﬂood emerged on the sandur plain in front of the glacier.

Its peak ﬂow was estimated as 45,000 m

−1

, and the ﬂood washed away a good

part of the road bridges across the sandur—which are built to withstand such ﬂoods.

What is of interest to us about this awesome ﬂood is that jökulhlaups occur reg-

ularly. Jökulhlaups from Grímsvötn over the last century have occurred at intervals

of 5–10 years, and the same cyclicity is evident in other such ﬂoods. Nor is this

evidential of cyclic volcanic activity. In fact, few of the Grímsvötn jökulhlaups are

associated directly with eruptions. Rather, the release of geothermal heat from the

caldera causes subglacial melting and a slow but regular rise in the lake level until a

ﬂood occurs. The water in a ﬂood ﬂows underneath the glacier, burning (somehow)

its way through the ice towards the glacier terminus. The fact that the 1996 ﬂood

took 10 hours from initial rupture of the seal at the caldera rim, until its emergence at

the glacier snout, indicates the travel time of the water over the intervening distance

of 50 km.

A. Fowler, Mathematical Geoscience, Interdisciplinary Applied Mathematics 36,

DOI 10.1007/978-0-85729-721-1_11, © Springer-Verlag London Limited 2011

741

742 11 Jökulhlaups

Fig. 11.1 The geometry of the lake and glacier, Grímsvötn and Skeiðarárjökull. (Figure courtesy

of Felix Ng)

What we want to do in this chapter is to build a model which will help us under-

stand how the ﬂood can make its way under the glacier, and perhaps what causes

the abatement of the ﬂood, usually after several weeks. We might also ask what

controls the recurrence time of 5–10 years, and we could consider other puzzles: in

a ﬂood, the lake is not drained completely. Nor (normally) is the ﬂood initiated (as

it was in 1996) when the lake is at ﬂotation level; normally it is about 60 meters

below this. We may not provide solutions to all these puzzles; but we will provide a

mathematical framework within which one can at least consider them.

11.1 The Nye Model

The basic fact that we need to know is that in temperate glaciers (those at the melt-

ing point, for example all Alpine glaciers), surface meltwater percolates to the bed,

and drains underneath the glacier in a subglacial hydraulic system which is usually

thought to consist of a network of channels. At the glacier snout, one or more emer-

gent streams often carry the discharge. The water of these pro-glacial streams is

usually milky in appearance, being highly charged with glacially derived sediments,

products of abrasion and comminution by the ice.

The classic Röthlisberger model of channel ﬂow (already introduced in Chap. 10)

conceptualises the subglacial streams as being cut into the ice, and existing at a

pressure p which is less than the overburden pressure p

of the ice. This is analogous

to the reduced pore pressure in soils, but it occurs for a very different reason. Firstly,

we expect the channel water pressure to be essentially that of the water at all parts of

11.1 The Nye Model 743

the bed—we visualise the boulder-strewn bed as being ‘leaky’, and allowing water

to migrate easily; if the water pressure were not locally uniform, rapid migration

would take place to ensure that it became so. But this pressure p must be less than

, otherwise ﬂotation would occur, and the glacier could advance catastrophically.

Now if p

>p, then the slow viscous creep of the ice will tend to cause channels

to contract, just as ﬂuid rushes into a cavity. Röthlisberger’s idea was that chan-

nels could be maintained open against this contraction by having the channel walls

melt—the source of the heat necessary to cause such melting being the frictional

heat released by the ﬂow of the water itself through the channel. This ingenious the-

ory is now generally accepted as being the mechanism whereby subglacial drainage

occurs, at least where there is evidence of distinct outlet streams, as in many valley

glaciers.

This theory is described below. Nye used it to allow consideration of time-

varying drainage, and thus developed a theoretical model for jökulhlaups, which,

in spite of its simplicity, is very successful. In his model, drainage is supposed to

occur through a single conduit, cut upwards from the bed into the ice, which has a

semi-circular cross section of area S. The two processes which affect S are the melt

rate m (mass per unit length downstream per unit time) and the viscous closure due

to the creep of the surrounding ice. The equation for S is thus

∂S

∂t

−KS(p

−p)

, (11.1)

where ρ

is the ice density, p is the channel pressure, and the second term represents

the creep closure due to a nonlinear ﬂow law ˙ε = Aτ

, where ˙ε is strain rate and

τ is stress. The constant K is proportional to A, but includes also some numerical

factors which arise through the exact solution of the ice creep problem. We see

that the equation for S introduces further variables m and N = p

− p, which is

(similarly to soils) called the effective pressure.

Now suppose that x represents a downstream spatial coordinate along the axis of

the channel (Nye used s, but we will reserve this for the ice surface, as in Chap. 10).

We assume the channel direction is slowly varying, which avoids the necessity of

including a complicated curvilinear coordinate system. We denote the volume ﬂux

of water through the channel by Q. Conservation of water mass in the channel then

requires

∂S

∂t

∂Q

∂x

+M, (11.2)

where m/ρ

represents the volumetric source due to side-wall derived melt, and M

is an additional source due to tributary ﬂow, surface meltwater supply, etc. We can

consider M to be prescribed.

We also require an equation describing momentum balance for the channel water.

We assume that the turbulent friction at the bed is given by a Manning correlation,

and we ignore inertial effects, which is apparently equivalent to assuming that the

Froude number is very small: at least during ﬂoods, it seems unlikely that this will

744 11 Jökulhlaups

be accurate, but in fact one can show (see Question 11.5) that in the present context,

inertial terms are always small. The Manning law then gives

g sin α −

∂p

∂x

=fρ

Q|Q|

8/3

, (11.3)

where f is a friction factor, and α is the mean bedrock slope (of the channel). The

relation of f to the Manning roughness factor n



is described below (in (11.19)).

The melt rate m is determined through an energy balance. The heat generated

by friction is Q[ρ

g sin α −

∂p

∂s

], and this is used to control water temperature, and

melt the side walls. If the average water temperature is θ

, while that of the ice is

, then an appropriate energy equation is



∂θ

∂t

∂θ

∂x





g sin α −

∂p

∂x



−m



L +c

(θ

−θ

)



. (11.4)

The left hand side is the material rate of change of water temperature with time, the

ﬁrst term on the right is the frictional source, and the second is the supply due to

the enthalpy change on melting. We can consider θ

to be known (for example, it is

the pressure melting point), but θ

must be further prescribed by a local heat trans-

fer condition at the ice wall across the thermal boundary layer there. An empirical

correlation for ﬂow in a cylindrical tube (the Dittus–Boelter correlation) gives



|Q|

1/2



0.8

k(θ

−θ

) =m



L +c

(θ

−θ

)



, (11.5)

where a

(≈ 0.2) is a constant, η

is the viscosity of water, and k is its thermal

conductivity.

Equations (11.1)–(11.5) give ﬁve equations for the ﬁve unknowns S,Q,m,p

and θ

. Certain initial and boundary conditions will also be appropriate; we will

come to these in due course.

11.2 Non-dimensionalisation

Equations (11.1)–(11.5) are non-dimensionalised by writing

Q =Q

∗

,S=S

∗

−p =N

∗

,m=m

∗

x =lx

∗

,t=t

∗

,θ

=θ

+θ

∗

(11.6)

The scales are chosen to effect the following balances in the equations. All three

terms balance in (11.1), thus

=KS

. (11.7)

Essentially, these give t

and S

. Nothing is balanced in (11.2), but we balance the

terms in (11.3) as follows. We write

g sin α −

∂p

∂x

=Φ +

∂(p

−p)

∂x

, (11.8)

11.2 Non-dimensionalisation 745

where

Φ =ρ

g sin α −

∂p

∂x

(11.9)

is the hydraulic gradient (if basal water pressure equals overburden ice pressure). If

z = b is the altitude of the base, then −∂b/∂x =sin α, and if the ice surface is at

z =s, then p

=ρ

g(s −b). Hence

Φ =−

∂

∂x



gs +(ρ

−ρ

)gb



, (11.10)

and we write

Φ =Φ

∗

, (11.11)

where

, (11.12)

and z = h

is a typical value of the Grímsvötn lake level (see Fig. 11.1). (At the

channel inlet, a value of water pressure p =p

corresponding to ﬂotation would give

g(h

−b) =ρ

g(s −b), so that the hydraulic head ρ

gs +(ρ

−ρ

)gb =ρ

there.) The balance of terms in (11.3) is now effected by having Φ ∼fρ

8/3

thus

fρ

8/3

. (11.13)

In (11.4), we balance the two terms on the right, thus

L, (11.14)

which ﬁxes m

, and we choose θ

by balancing the advective term with the source

term in (11.4):

. (11.15)

The ﬁve relationships in (11.7), (11.13), (11.14) and (11.15) determine scales for

, m

, t

, N

and θ

in terms of Q

, which for the moment we presume is pre-

scribed; l is given from the geometry. It is natural to prescribe Q

by balancing ﬂux

with source in (11.2), i.e. Q

=Ml; indeed, this is an appropriate choice for steady

subglacial drainage of Röthlisberger type, but as we shall see, it is inappropriate for

the violence of jökulhlaups, where the typical ﬂow rates are much larger.

Given Q

, we ﬁnd



fρ



3/8

=(Kt

)

−1/n

(11.16)

and the dimensionless equations (11.1)–(11.5) become (dropping the asterisks on

the variables)

746 11 Jökulhlaups

Table 11.1 Physical

parameter values

Symbol Value

4.2kJkg

−1

f 0.05 m

−2/3

g 9.8ms

−2

1.5km

k 0.56 W m

−1

K 0.5 ×10

−24

−n

−1

l 50 km

L 3.3 ×10

kJ kg

−1

n 3

1.8 ×10

−1

2 ×10

−3

Pa s

0.9 ×10

kg m

−3

kg m

−3

300 Pa m

−1

∂S

∂t

=m −SN

∂S

∂t

∂Q

∂x

=εrm +Ω,

Φ +δ

∂N

∂x

Q|Q|

8/3

, (11.17)

εS

∂θ

∂t

∂θ

∂x



Φ +δ

∂N

∂x



−m(1 +εrθ),



|Q|

1/2



0.8

=γm(1 +εrθ),

and the parameters ε, r, Ω, δ and γ are deﬁned, after some algebra, by

ε =

,δ=



1/4

11/8

KL(fρ

3/8



1/n

γ =





0.8

1/2



fρ



3/20

, (11.18)

r =

,Ω=

To estimate the sizes of the scales and parameters, we use the values in Ta-

ble 11.1. The value of the closure rate coefﬁcient K is given by K =2A/n

, and we

use Paterson’s (1994) recommended value of A =6 ×10

−24

−n

−1

at 0°C, with

11.3 Boundary Conditions and Lake Reﬁlling 747

n =3. The deﬁnition of the friction factor f in Manning’s roughness law as written

here is

f =n

2





2/3

, (11.19)

where n



is the roughness coefﬁcient, and R

is the hydraulic radius (= S/l,

where l is the wetted perimeter). For a semi-circular channel (S/R

)

2/3

=(2(π +

/π)

2/3

≈ 6.6, so that if n



= 0.09 m

−1/3

s, then f ∼ 0.05 m

−2/3

(and

0.01 m

−2/3

for n



= 0.04 m

−1/3

s). The value of Φ

follows from (11.13), and

we have anticipated the value of Q

which is found below; we then ﬁnd succes-

sively that

∼10

∼163 kg m

−1

,θ

∼3.6K,

∼0.6 ×10

s (0.68 day), N

∼3.2 ×10

Pa (32 bars),

(11.20)

and the dimensionless parameters are of typical sizes

γ ∼2.5,ε∼0.05,r∼0.9,δ∼0.22,Ω∼0.6 ×10

−3

, (11.21)

where for Ω we assume a base ﬂow rate of Ml = 10

−1

, which represents a

typical value of the discharge between jökulhlaups. We see that all the parameters

are O(1), which indicates that the scales we have chosen are sensible. Next we

need to ﬁnd a reason why Q

should be chosen as large as 1.8 ×10

−1

11.3 Boundary Conditions and Lake Reﬁlling

The equations in (11.17) require initial conditions for S and θ , and two boundary

conditions for Q or N, and one for θ . The boundary condition for θ is taken to be

θ =θ

at x =0 (11.22)

(at least when Q>0 at the lake). At the outlet, it seems we should prescribe

N =0atx =1, (11.23)

i.e., the water pressure becomes atmospheric

(as also does the ice pressure). At the

inlet to the channel, conservation of mass requires that

−Q(0,t), (11.24)

where V is the lake volume and m

represents the geothermal melt rate. Suppose

the lake level is at z =h, and open to the atmosphere.

We assume V = V(h), and

There are problems with this, however: see further discussion in the notes.

At Grímsvötn, this is normally the case. Apparently the high geothermal heat levels melt ice at the

caldera walls, so that water is usually present at the walls all the way to the surface. See Björnsson

(1988, pp. 70–73).

748 11 Jökulhlaups

in fact V



(h) =A

is the lake surface area (which may depend on h). Now the water

pressure at the inlet is ρ

g(h − b), where b is the grounding line elevation of the

bed, and this is equal to p

−N . Therefore (if b and p

do not vary)

=−

(0,t), (11.25)

and thus the boundary condition at the lake inlet is

−

∂N

∂t

−Q at x =0, (11.26)

where x = 0 is taken to be the position of the lake margin. We now ﬁnally choose

the scale Q

so that

, (11.27)

where we will take A

as constant. Using the deﬁnitions of N

and t

in (11.16), we

obtain







fρ



3/8



n+1



4/(3n−1)

. (11.28)

It follows that (11.26) can be written in the dimensionless form

∂N

∂t

=Q −ν at x =0, (11.29)

where

ν =

. (11.30)

We take A

=30 km

=3 ×10

, and use other values given previously; hence

we compute Q

∼1.8 ×10

−1

, as we assumed previously.

The reﬁlling rate m

can be estimated from the total rate of discharge, which is

about 5 ×10

kg y

−1

,or17m

−1

. This gives a typical value of ν ≈0.89 ×10

−4

We can guess that ν  1 controls the (large) length of the time period between

jökulhlaups. Over a catchment area of perhaps 2 × 10

, this reﬁlling rate

corresponds to a melt rate of 0.8 × 10

−7

−1

. With ρ

= 10

kg m

−3

and

L = 3.3 × 10

Jkg

−1

, the geothermal heat ﬂux required to provide this is about

2.4 W m

−2

; this compares with a typical (non-volcanic) value of 0.05 W m

−2

!The

heat is not, of course, delivered as a conductive heat ﬂux, but via the upwelling of

superheated geothermal ﬂuid (essentially sulphurous groundwater).

11.4 Simpliﬁcation of the Model

The parameters ε, δ and Ω are all relatively small. If we neglect them, it follows

that (supposing Φ and Q>0)

Q ≈Q(0,t), i.e. Q ≈Q(t),

S ≈Φ

−3/8

3/4

(11.31)

11.4 Simpliﬁcation of the Model 749

and hence that

m ≈

3/20

1/2

, (11.32)

and thus

γQ

1/2

∂θ

∂x

≈γΦQ

1/2

−θΦ

3/20

. (11.33)

We see from this that θ approaches a limiting value such that m = ΦQ over a di-

mensionless distance of order γQ

1/2

. Now since Q was scaled with a value which

is somewhat larger than a typical peak discharge, it is certainly clear that between

ﬂoods (when we expect Q ∼Ω 1) this distance will be very short; it might even

be short during ﬂoods. To accommodate this suggestion, we simplify the model by

assuming that m =ΦQ holds at all times, even though this may be inaccurate for

a short time during maximum discharge.

It then follows that N ≈N(0,t)=N(t),

and N and Q satisfy

Q =

11/8

5/4

−

N =Q −ν,

(11.34)

where

N =dN/dt. The model thus reduces to the solution of two ordinary differ-

ential equations!

If we take Φ to be constant, Eqs. (11.34) are easily studied in the (N, Q) phase

plane.

There is a ﬁxed point at Q = Q

∗

=ν, N =N

∗

=Φ

11/8n

∗1/4n

, and if we

write

Q =Q

∗

+q, N =N

∗

+Π, (11.35)

then linearised equations for Π and q are



˙q



≈



∗n

−

∗

∗(n−1)





, (11.36)

and solutions are proportional to e

λt

, where

−

∗n

λ +

∗

∗(n−1)

=0, (11.37)

which implies that the equilibrium is an unstable spiral (if

ν>

11(n+1)

2(3n−1)

(48n)

3n−1

≈0.58 ×10

−3

for n = 3; (11.38)

There are other considerations which suggest that a useful value of γ , at least during a violent

ﬂood, may be smaller than the value in (11.21); in this case the assumption of equilibrium of

(11.33) may indeed be realistic. See also the discussion in the notes.

One might wonder how the assumption that Q and N in (11.34) are functions only of t can be

squared with a hydraulic potential gradient Φ(x) which depends on distance downstream. In fact,

because the lake reﬁlling condition is applied only at the channel inlet, the value of Φ in (11.34)is

actually that at the inlet, i.e., Φ(0). The case where this is negative is considered below.