Fowler A. Mathematical Geoscience

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

750 11 Jökulhlaups

Fig. 11.2 Phase portrait for

N and Q,withν =0.1in

(11.34)

Fig. 11.3 Time series for Q

corresponding to Fig. 11.2

otherwise an unstable node): steady drainage from a lake is always unstable, in this

model.

Figure 11.2 shows a numerical solution of (11.34)inthe(N, Q) phase plane.

Clearly the spiral structure continues for (N, Q) away from the unstable ﬁxed

point. The time series corresponding to this diagram (Fig. 11.3) shows a sequence

of jökulhlaups of growing discharge, with long intervals (of O(1/ν)) between the

ﬂoods. In the following section, we comment on the growing amplitude of these

oscillations. If we focus on a single ﬂood, there is clearly no criterion to pick which

hydrograph will occur, and this is a drawback of the model. If we select initial values

and parameters to ﬁt the rising limb of the 1972 jökulhlaup, we ﬁnd (see Fig. 11.4)

that the peak discharge and decay curve do not ﬁt the data. The resolutions of these

conundrums are given in the following sections.

11.5 Periodic Oscillations

The main problem with the result of the previous section is that there is no limit

to the growth of the jökulhlaup amplitude with time. In particular, the approximate

model allows N to become negative between ﬂoods (although Q always remains

positive). This is unphysical, for if p>p

in the channel, then in reality leakage

would occur laterally along the glacier bed, and ﬂotation would occur.

11.5 Periodic Oscillations 751

Fig. 11.4 Hydrograph from

the 1972 jökulhlaup (crosses)

and a ﬁt obtained by solving

(11.34)

Physically, this implies that when the effective pressure N becomes zero, then the

closure equation (11.1) becomes inappropriate for determining S. We thus imagine

the following scenario. When N reaches zero, the channel can spontaneously open

(by lifting off the ice), and we might model this by imagining a very rapid ‘creep’

opening of S if p approaches p

. A suitable modiﬁcation of (11.1) could then be

∂S

∂t

−KSf(N), (11.39)

where f ≈ N

for N>0, but f →−∞very rapidly when N → 0. The non-

dimensional version of (11.33)

is then

Q =Q

5/4

−Qf (N), (11.40)

but now

Q  1ifN ≈ 0. In the limit, this simply squashes the trajectories in

Fig. 11.3 into N>0, and leads to a modiﬁed phase portrait in which, however,

limit cycle periodic behaviour is still not obtained.

There are a number of other possible weaknesses in this simple discussion of

the Nye model. The difﬁculty associated with allowing N = 0 also crops up at

the glacier snout, where the imposition of N =0 would cause unbounded channel

growth, according to (11.17)

. Other issues which merit discussion are our neglect

of thermal advection, and the assumption of a semi-circular channel.

11.5.1 Breaking the Seal

There is, however, a far more fundamental ﬂaw in the simple analysis given above,

and it is this: the existence of a subglacial lake in the ﬁrst place implies the existence

of a local minimum in the hydraulic potential, and thus implies that the basic hy-

draulic gradient Φ must be negative in the vicinity of the lake. Consulting (11.17),

we see that the neglect of ε, δ and Ω then implies that there is a water divide, or seal,

near the lake. At the lake, Φ<0, and water ﬂows back towards the lake, whereas

752 11 Jökulhlaups

further down the ice slope, Φ>0 and drainage is towards the ice margin. This situ-

ation cannot be maintained, because then N decreases at the lake and reaches zero,

at which point we suppose that a ﬂood is initiated.

But even this description is too simplistic, at least in the case of Grímsvötn. Nor-

mally the ﬂoods from Grímsvötn are initiated when the lake level is 60–80 metres

below ﬂotation, and N at the inlet is inferred to be in the region of 6–8 bars. The

initiation of channelised ﬂow must therefore normally begin by a mechanism other

than ﬂotation, when there is still a hydraulic barrier to the lake at the caldera rim.

A related consideration is that mathematically, the neglect of δ in the momentum

equation (11.17)

is problematical, for the following reason. Neglecting ε,integra-

tion of the mass conservation equation implies

Q =Ω(x −x

∗

), (11.41)

where x

∗

is ﬁxed by the prescription that Φ(x

∗

) = 0 (assuming δ = 0). The mo-

mentum balance equation then gives S, so that the closure equation gives N, and

the variables are independent of t. In general, the lake reﬁlling condition is not sat-

isﬁed, and this suggests that the loss of δ is a singular perturbation, which requires

the consideration of a boundary layer in x.

We analyse such a boundary layer by writing x = δX, and allow for the possibil-

ity of reversed ﬂow at the seal, whose position, however, is no longer constrained to

be where Φ =0. (11.17) is written approximately as

∂S

∂t

|Q|

8/3

−SN

∂Q

∂X

=ω =δΩ, (11.42)

Φ +

∂N

∂X

Q|Q|

8/3

with the boundary conditions that

∂N

∂t

(0,t)=Q(0,t)−ν,

∂N

∂X

→0asX →∞.

(11.43)

The last of these is a matching condition, which enables N to be matched to a far

ﬁeld slowly varying proﬁle.

Between ﬂoods, or in conditions of normal drainage, Q  1, and (11.43) then

shows that N varies slowly, and (11.42)

implies that S relaxes to equilibrium (if

N =O(1))overatimescaleofO(1). Thus

S ≈



|Q|



3/11

Φ +

∂N

∂X

≈

8n/11

|Q|

2/11

sgnQ.

(11.44)

We can use this particular simpliﬁcation to understand the mechanism of ﬂood initi-

ation. Since we suppose Φ<0 at the lake, but Φ>0 far from the lake, it is natural

to choose the distinguished limit in which Φ varies over the same length scale X as

is appropriate in the boundary layer. Indeed, this is the case for Grímsvötn.

11.5 Periodic Oscillations 753

Fig. 11.5 N

versus X

∗

when Φ =1 −ae

−bX

,for

a =4, b = 2 (strong seal) and

a =4, b = 4 (weak seal). See

also Question 11.2

A Particular Example

In general, we need to solve the system (11.42), or (11.44), numerically. To gain

some analytic insight, it is instructive to simplify (11.44) by replacing the exponent

8n/11 in (11.44) by 1, and ignoring the denominator |Q|

2/11

. We thus have to solve

Φ +

∂N

∂X

=N sgnQ, (11.45)

where

Q =ω(X −X

∗

), (11.46)

and we require (assuming Φ(∞) =1) that

N →1asX →∞,

N =N

on X =0,

=−(ν +ωX

∗

(11.47)

The seal position X

∗

is unknown, and is a function of time. It is consistent with the

slow variation of N that X

∗

will vary slowly with t.

One quickly ﬁnds that it is necessary that X

∗

> 0 (there is a seal) in order that N

does not grow exponentially at +∞. The solution is then

N =N

−X

−



Φ(ξ)e

−(X−ξ)

dξ, 0 <X<X

∗

N =



∞

Φ(ξ)e

−(ξ−X)

dξ, X > X

∗

(11.48)

Continuity of N at X

∗

thus requires

∗

) =e

∗



∞

Φ(ξ)e

−|X

∗

−ξ|

dξ, (11.49)

and (11.47)

then provides an ordinary differential equation for the evolution of N

(or X

∗

): N

decreases with time.

It is clear from our supposition that Φ →1asX →∞that N

is positive when

∗

is large; our concern is then what happens as N

decreases in accord with

754 11 Jökulhlaups

(11.47). There are two possibilities, and these are illustrated in Fig. 11.5.Inthe

ﬁrst, N

(0)<0, and thus N

reaches zero when X

∗

is still positive, that is, while

the seal position is still in front of the lake margin. Thus ﬂotation occurs at the lake,

and a ﬂood ensues because of this. This appears to have been the situation in the

1996 Grímsvötn ﬂood. We term a seal arising from a hydraulic gradient such that

(0)<0astrong seal.

In contrast, a weak seal is one for which N

(0)>0. In this case, X

∗

reaches zero

while N

is still positive, i.e., the lake level is below ﬂotation. When the seal reaches

the lake margin, the lake starts to empty, and a ﬂood is initiated. This appears to be

the normal case in Grímsvötn, where the lake level is typically about 60 metres

below ﬂotation level at ﬂood onset, and the effective pressure is some 6 bars.

To summarise: we characterise a strong seal (in the context of (11.45)) as one

with a hydraulic gradient such that



∞

Φ(ξ)e

−ξ

dξ < 0, (11.50)

and a weak seal as one with



∞

Φ(ξ)e

−ξ

dξ > 0. (11.51)

For a strong seal, the seal is broken when the lake level rises to ﬂotation, but for

a weak seal, the drainage divide slowly migrates backwards as the lake reﬁlls, and

reaches the lake when it is still below ﬂotation; at this point the seal is broken and

the next ﬂood is initiated. The simple formulae (11.50) and (11.51) only apply to

the simpliﬁed model (11.45), but they do indicate the essential difference between

a strong seal and a weak seal, which is that the strong seal has a larger negative

hydraulic gradient near the lake. In the exercises, Question 11.2 calculates N

∗

)

for the representative choice

Φ =1 −ae

−bX

, (11.52)

from which we ﬁnd that a seal is weak if

1 <a<a

=b +1, (11.53)

and is strong if a>b+ 1. (This is the hydraulic gradient used in the illustrative

Fig. 11.5.) Question 11.6 extends the analysis to the better approximating equation

Φ +

∂N

∂X

sgnQ; (11.54)

for this equation (and a hydraulic potential gradient given by (11.52)) we can show

that a seal is weak if

a<a



(2/b),1



, (11.55)

where j



ν,k

is the kth zero of the Bessel function derivative J



(z). We can use the

asymptotic expansions of j



ν,1

for large and small ν,



ν,1

≈



ν +0.81ν

1/3

as ν →∞,

√

2ν as ν →0,

(11.56)

11.5 Periodic Oscillations 755

Fig. 11.6 N

(t) and the outlet discharge Q

out

(t) = Q

(Ω − ωX

∗

) with ω = 0.12 × 10

−3

ν = 0.89 ×10

−4

, a =2.8, b = 4.316 in (11.52). To convert the results to dimensional quantities,

we have used the scales Q

=0.18 ×10

−1

, baseﬂow ΩQ

=100 m

−1

, N

=32.37 bars,

= 0.0019 y, δ = 0.216, l = 50 km. These are our estimated parameter values, except that the

hydraulic gradient is less negative near the lake than it appears to be in reality. The (dimensionless,

with δl and t

) space and time steps used were 0.005 and 0.0005, respectively

to derive approximations for a

for large and small b:

≈



[1 +0.51b

2/3

]

as b →0,

b as b →∞.

(11.57)

Comparison of these approximations with the exact expression (11.55) show that

the small b approximation is very accurate for b<1, and quite accurate for b<4;

and in this range the small b approximation itself is close to 1.2 +b; thus a

and

are in fact quite close, and we might reasonably expect that this estimate for the

occurrence of a weak seal is quite accurate also for (11.44)

Actually, a better approximation allows for the smallness of Q by replacing

2/11

in (11.44)

by ω

2/11

, and then we ﬁnd that that the condition for a weak

seal is modiﬁed to the estimate

a<a



1/11



(2/ω

1/11

b),1



≈1.2 +0.44b (11.58)

for ω =0.12 ×10

−3

(see Question 11.6). This appears to be consistent with what is

observed numerically.

Figure 11.6 shows a numerical solution of (11.42) and (11.43), using the hy-

draulic gradient given by (11.52) with a = 2.8, b = 4.316.

With these values,

and the other parameters as estimated earlier, we ﬁnd ﬂoods with peak discharge

5032 m

−1

, and periods 7.4 years, very favourably comparable to observed peak

Reasonable estimates of Φ for Grímsvötn yield a =3.33, b =4.316, given our choice of Φ

and

. For these values, the seal is actually strong, but would be weak according to (11.58)ifa<3.1:

hence our modiﬁed choice of a so that the ﬂoods are weak. A likely modiﬁcation to the predicted

type of seal in (11.58) arises from the effects of water temperature near the lake; this is discussed

further in the notes.

756 11 Jökulhlaups

Fig. 11.7 The ﬂood

hydrograph, corresponding to

Fig. 11.6 (solid curve). Also

shown is the hydrograph of

the 1972 jökulhlaup (dashed),

with the discharge scaled by

0.7, and the time scaled by

2.19. The agreement is

evidently excellent. The

rescaling corresponds to

choosing

=0.26 ×10

−1

and

=0.00087 y = 0.32 d. See

also Question 11.7

discharges and periods (but it should be mentioned that these have varied over the

years). In addition, the ﬂood duration is comparable to that generally observed (it

is about twice as long), and the shape of the ﬂood hydrograph is now very similar

to that of the 1972 jökulhlaup (see Fig. 11.7). We can also see from Fig. 11.6 that

ﬂoods are initiated when the lake is below ﬂotation, with N

≈ 1.7 bars, about a

quarter the value observed, but this depends sensitively on the hydraulic parameters

a and b. The choice of a and b is made to illustrate the pre-ﬂotation ﬂood initia-

tion, but their values depend on the scales N

and Φ

. In view of the simpliﬁcations

made, which we discuss further below, this theory appears to give a good represen-

tation of the observations, and by tweaking the choice of parameters, it is of course

possible to improve the agreement.

The 1996 Eruption

One of the surprises of the 1996 eruption was that the seal did fail at ﬂotation. It had

been expected to fail at the usual 60–80 m below ﬂotation. We can guess that the

reason for this was that the eruption ﬁlled the lake in a relatively short time, weeks

rather than years, and therefore the equation for S in (11.42) could not be taken as

being in equilibrium. For example, if the lake reﬁlling were so rapid that ν  ω

(see Question 11.8), then N

decreases rapidly, and it is reasonable to expect N

to reach zero before the seal breaks, though this in fact requires a time-dependent

solution of a model such as (11.42) and (11.43). In normal circumstances, slow lake

reﬁlling causes a cyclic sequence of jökulhlaups when N

> 0, but if the reﬁlling

rate is dramatically increased for a short period, such as would occur following a

volcanic eruption, then there is a violent ﬂood which occurs when the lake reaches

ﬂotation.

Figure 11.8 shows a numerical experiment which illustrates this fact. The same

model with the same parameters as in Fig. 11.6 is solved, but the value of ν is

11.5 Periodic Oscillations 757

Fig. 11.8 Outlet discharge

for a post-eruption ﬂood;

parameters as for Fig. 11.6,

but ν =0.56 ×10

−2

for

1 <t<1.1 y, corresponding

to a reﬁlling rate of

1000 m

−1

for about a

month, which resembles the

conditions in 1996. The

magnitude of the ﬂood which

is obtained numerically

depends on the timing and

magnitude of the enhanced

meltwater pulse

changed to be 0.56×10

−2

between t =1 y and t =1.1 y. This corresponds to a ﬂux

of 10

−1

for about a month, which corresponds to the situation after the 1996

eruption, when 3 cubic kilometres of water entered the Grímsvötn caldera in the

month preceding the ﬂood. The peak discharge is now over 35,000 m

−1

, which

compares favourably with the inferred peak of 45,000 m

−1

(given our earlier

comment about Q

), and although the shape of the hydrograph is different (the 1996

ﬂood lasted only a day), it is very similar to hydrographs of earlier post-volcanic

ﬂoods, particularly those of 1922 and 1934, each of which lasted about ten days

(0.027 y), and which had similarly shaped hydrographs, with slow rise and fast

decay.

11.5.2 Wide Channels and the 1996 Eruption

In fact, the detailed behaviour of the 1996 ﬂood was quite different from the normal,

viscously controlled ﬂoods, since the water hydro-fractured its way along the base.

It is possible to build an understanding of the discrepancy between the slowly rising

ﬂood hydrograph of the usual weakly sealed ﬂoods and the rapidly rising strongly

sealed ﬂoods, or eruptive ﬂoods, by consideration of the way in which the ﬂood wa-

ter propagates under the glacier. When the seal is broken by overpressuring of the

channel water, then the channel widens essentially as a crack does, and the channel’s

lateral extension must be described by hydrofracture (see also Sect. 9.5). The conse-

quent wide channel propagates downglacier by a hydrofracturing overpressured tip,

with the tail of the channel being a normal, underpressured, viscously closing wide

conduit. Observations of the 1996 ﬂood are consistent with this, since collapse of

the ice post-ﬂood indicated a wide channel near the lake of several hundred metres,

758 11 Jökulhlaups

and the ﬂoodwater burst out from the surface of the glacier several kilometres up-

stream of the snout (presumably along a thrust fault), indicating an overpressuring

of several bars.

11.6 Cauldrons and Calderas

One of the plausible assumptions about our implementation of the Nye model, but

one that may not always be realistic, is that the lake in the Grímsvötn caldera is open

to the atmosphere. Sometimes this may be the case, and then our direct relation of

lake volume to effective pressure at the lake margin is reasonable. But it is also rea-

sonable to suppose that the overlying ice is unbroken, forming a miniature ice shelf

over the caldera. In this case, we cannot simply relate water pressure to lake volume

during a ﬂood. As the lake empties, the overlying ice must deform to accommodate

the volume loss, and this deformation must be related to the effective pressure in the

lake. In this section we show how to compute this relationship.

The bowing of the ice surface results in the formation of a cauldron. Caul-

drons are common where subglacial volcanoes occur, and an example of a caul-

dron formed in the 1996 Vatnajökull eruption is shown in Fig. 11.9. The eruption

causes massive subglacial melting, and the outﬂow of this water (in this case, into

the Grímsvötn caldera) causes the subsidence that is seen.

11.6.1 Viscous Beam Theory

We can analyse the ice cauldron deformation using a viscous analogue to the beam

theory of classical elasticity. We begin by recalling the equations of two-dimensional

ice motion, scaled as for ice sheets in (10.38):

=0,

0 =−s

+τ

+ε

[−p

+τ

0 =−p

+τ

−τ



+ε



=Aτ

n−1

=Aτ

n−1

=τ

+ε

(11.59)

The ice surface is at z = s, and the lake roof is at z =h. The boundary conditions

are, on z =s:

+ε

(p −τ

=0,

+p +τ

=0, (11.60)

w =s

+us

−a,

while at the lake z =h,

11.6 Cauldrons and Calderas 759

Fig. 11.9 An ice cauldron forming during the Gjálp eruption under Vatnajökull in 1996. The

cauldron is about two kilometres in diameter and 100 metres deep, and ring shear fractures can be

seen, indicating yield of the ice. The subsidence rate of the ice surface was initially about 12 metres

per hour. Photo courtesy of Magnús Tumi Gudmundsson, University of Iceland, obtained from

http://www.hi.is/~mmh/gos/photos.html. For further information see Gudmundsson et al. (2004)

+ε

(p −τ

=−γNh

+p +τ

=−

γN

, (11.61)

w =h

+uh

−m,

where

γ =

, (11.62)

d being the ice depth scale, as in Chap. 10; γ is of O(1) or smaller. The terms a and

m represent surface accumulation rate and basal melting rate, respectively.

We suppose in addition that pressure in the lake is hydrostatic. This condition

(more precisely that there is no hydraulic gradient in the lake) can be written in

terms of scaled variables as

∂

∂x

[s +δ

h −γN]=0, (11.63)

where

−ρ

, (11.64)