Fowler A. Mathematical Geoscience

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

760 11 Jökulhlaups

and the lake effective pressure is deﬁned (as for a channel) as the cryostatic overbur-

den ice pressure minus the water pressure. The deﬁnition of δ

is identical to that

of δ in Chap. 10, but the symbol δ in this chapter has been reserved (see (11.18)).

The ﬁnal equation is the lake reﬁlling equation, which balances the net outﬂow

from the lake with its rate of loss of volume. In dimensional terms, this can be

written in the form

−Q

−

=−



wdS, (11.65)

where A

is the lake area. Its dimensionless form is given later.

Normally, we would rescale the equations into a form appropriate for an ice shelf,

but when rapid subsidence occurs, this is inappropriate. Rapid deﬂation of an ice

cauldron requires a positive effective pressure in the lake in order to suck the ice

down; for an ice shelf (see Sect. 10.2.6), we would suppose that a ﬂotation condition

applies.

The rescaling we choose balances different terms in the equations, and is analo-

gous to that used when scaling the equations of elasticity for a beam. We rescale the

variables as follows:

p,τ

∼Λ, τ ∼εΛ, τ

∼

w ∼W, u ∼

,x∼λ, t ∼

(11.66)

The length scale λ is prescribed (from the lake geometry), but Λ, μ and W must be

chosen appropriately. We expect Λ 1, and that W  1 is determined by the rate

of lake deﬂation. We take the deﬂation rate scale to be w

, so that (compare (10.33))

W =

[a]

, (11.67)

and from (11.65), we deﬁne w

(and thus W )by

. (11.68)

The parameter μ is chosen so that the new time scale is that of the channel ﬂow, t

in (11.6); this leads to the deﬁnition

μ =

. (11.69)

With these rescalings, the equations take the form

=0,

0 =−s

+ε

Λ[τ

−p

+τ

0 =−p

−τ

+ν

=ν



n−1



Aτ

n−1



n−1



Aτ

n−1

=ν

+τ

(11.70)

11.6 Cauldrons and Calderas 761

where

ν =

, (11.71)

and this will be small.

The boundary conditions become, on the surface z =s,

+(p −τ

=0,

+p +τ

=0, (11.72)

=μ



w −ν



and on the lake z =h,

+(p −τ

=−

+p +τ

=−

N, (11.73)

=μ



w −ν



and the lake hydrostatic condition is still (11.63). The lake reﬁlling equation (11.65)

becomes in dimensionless terms

−Q

−

=−



−

wdx, (11.74)

assuming a two-dimensional geometry; x

−

and x

are the upstream and downstream

positions of the lake margins.

Equations (11.70) describe the deformation of the viscous beam, with the bound-

ary conditions in (11.72) and (11.73) sufﬁcing also to determine the evolution of

the ice surface s and the lake roof h, providing the effective pressure N is known:

this is determined by the buoyancy condition (11.63). The equations and boundary

conditions include six parameters ν, ε, Λ, W , γ and μ; ε is the ice sheet aspect ra-

tio, deﬁned in (10.37), ν is the beam aspect ratio, deﬁned in (11.71), γ was deﬁned

in (11.62), μ in (11.69), and W in (11.67); this leaves the parameter Λ yet to be

chosen.

For Grímsvötn, the lake diameter is of order 5 km, while the Vatnajökull ice cap

is of depth ca. 500 m, and is some 100 km in extent. Thus ε ∼0.005 and λ ∼0.05;

hence ν ∼ 0.1 in this case, and generally we assume it to be small. Since ν is the

aspect ratio for the ice ‘beam’, the beam equation follows from the assumption that

it is small.

The parameter Λ is determined by requesting a balance in the constitutive law

for longitudinal stress, thus we choose

n−1

=1; (11.75)

A is the relevant value of the ﬂow rate coefﬁcient A, and is explicitly included

because the rate controlling value of A is the smallest, i.e., that at the surface, which

762 11 Jökulhlaups

in an ice sheet may be three orders of magnitude smaller than that at the base. We

also deﬁne a parameter

β =



n+1



1/n

. (11.76)

When written in terms of local dimensional quantities, this is

β =

, (11.77)

where the beam stress scale is





1/n

, (11.78)

with A

being the surface ﬂow law coefﬁcient (thus

A = A

/[A],cf.(10.33)) and

the surface deformation rate scale. For a subsidence rate appropriate to Fig. 11.9

of 12 metres an hour, w

≈10

−1

, and with A

=10

−2

bar

−3

−1

, d =500 m

and ν =0.25, this gives τ

≈11 bars and β ≈4. In other circumstances, the beam

stress is lower, and the value of β may be large.

From the equations and boundary conditions, we see that p +τ

∼ν

, and there-

fore we deﬁne

p +τ

=−ν

. (11.79)

Using (11.79), (11.75) and (11.76), the rescaled equations become

=0,

p +τ

=−ν

0 =−βs

+τ

+2τ

+ν

0 = σ

+τ

, (11.80)

=ν

Aτ

n−1

=Aτ

n−1

=ν

+τ

and the boundary conditions become, on z = s,

−



2τ

+ν



=0,

−σ

=0, (11.81)

≈μ



w −ν



and on the lake z =h,

−



2τ

+ν



=−γβNh

−σ

=−

γβN

=μ



w −ν

∗



∂

∂x

[s +δ

h −γN]=0,

(11.82)

11.6 Cauldrons and Calderas 763

taking W  1, but retaining the melt term m

∗

=m/W since it may be the increase

in basal melt rate which causes cauldron formation.

While all the parameters in the equations have been deﬁned, the choice of the

volume ﬂux scale Q

, just as in the Nye model, is still open. It is natural to balance

terms in the hydrostatic equation (11.82)

; bearing in mind (11.81)

, this suggests

μ ∼γ , and in fact our earlier choice in (11.27) corresponds to the choice

μ =

1 +δ

, (11.83)

and we will see that it remains appropriate to deﬁne μ in this way for closed subsur-

face lakes. The three other parameters involved in ﬁnding an approximate solution

to this set of equations are γ , β and ν. We have indicated that γ  1, ν  1 and

β  1. We therefore begin by seeking approximate solutions for ν  1, and then

subsequently considering possible choices for γ , β and μ. The limit ν →0isthe

classical approximation associated with beam theory.

With ν 1, we have w ≈w(x, t), τ ≈|τ

| and u

≈0, whence

≈A|τ

n−1

u ≈V −zw

(11.84)

where V(x,t)is to be determined. The coefﬁcient A has now been rescaled with

if this is not equal to one. From (11.84), we ﬁnd





1/n

−zw

(n−1)/n

. (11.85)

We now deﬁne the three quantities M, the bending moment, S, the shear force,

and T , the tension, as

M =−



2zτ

dz, S =



dz, T =



dz. (11.86)

We also deﬁne the secondary quantities U, the uplift force, and L, the lifting torque,

to be

U =



dz, L =



zσ

dz. (11.87)

Integrating the force balance equations and applying the boundary conditions on

z =s and z =h, we then ﬁnd that

∂T

∂x

+ν

∂U

∂x

=−γβNh

+β(s −h)s

∂M

∂x

−ν

∂L

∂x

+S =γβNhh

−



−h



, (11.88)

∂S

∂x

γβN

Subject to suitable boundary conditions at the margins x

, these equations will en-

able us to provide a closed solution which determines the surface depression rate

−w in terms of the underlying effective pressure N .

764 11 Jökulhlaups

Some care needs to be taken with this discussion. The rescaling of the stresses

in (11.66) is such that they may be large. Comparison with the original choice of

scales in (10.9) shows that the greatest stress is the longitudinal stress, which is of

order τ

as given by (11.78), and this can be some way larger than typical

stress magnitudes in glaciers. Since the yield strength of ice (before it fractures)

is generally thought to be of the order of bars, it is clear that the stresses in the

viscous ice beam may exceed this. In this case, we can expect the ice to fracture, so

that its effective rheology becomes, for example, plastic. Indeed, we see extensive

fracturing in Fig. 11.9.

Marginal Boundary Conditions

We forgo discussion of such plastic subsidence, and assume that the ice retains its

integrity. In order to motivate plausible boundary conditions, we consider the lake

margins x

to be ﬁxed, so that

w =0atx =x

. (11.89)

In addition, we compute typical viscosity scales η

for the ice sheet ice, and η

for

the overlake ice. These are taken to be

2[A][τ ]

n−1

,η

n−1

, (11.90)

where the ice sheet scales are those deﬁned in (10.33). We take [A]=0.2 bar

−3

−1

[τ ]=0.2 bar, which yields η

≈60 bar y, while for the overlake ice, we take A

−2

bar

−3

−1

, τ

= 11 bars, yielding η

= 0.4 bar y. On this basis, we see that

the overlake ice can be much less stiff than the ice sheet ice, and then it seems

appropriate to pose conditions of no horizontal velocity at the margins, thus

V =w

=0atx =x

. (11.91)

Because (11.88) are vertically averaged equations, we might expect not to be able

to satisfy both conditions in (11.91), but simply an integrated condition of no hori-

zontal mass ﬂux at the margins. Using (11.84)

, this condition can be written in the

form

V =

(s +h)w

at x =x

. (11.92)

In fact, we shall see that both conditions in (11.91) can be satisﬁed.

The boundary conditions (11.89) and (11.92) are applicable for soft lake ice

which may be appropriate for Vatnajökull, but in other circumstances, the lake ice

can be hard. This is probably the case in Antarctica, where the surface value of A

will be much lower, and the beam stress τ

is also lower. In this case it is appropri-

ate to apply conditions of zero longitudinal stress at the margins, i.e., τ

=0, and if

we suppose that this condition can be applied in a vertically integrated sense, then

M =T =0atx =x

. (11.93)

11.6 Cauldrons and Calderas 765

Summary

From the rescaled equations of motion (11.80), together with the upper and lower

boundary conditions (11.81) and (11.82), we have derived vertically integrated

equations (11.88), together with suitable lateral boundary conditions, either (11.89)

and (11.92) for soft lake ice (Vatnajökull), or (11.93) for hard lake ice (Antarctica).

We now need to use the beam approximation based on ν 1 to derive the viscous

beam equation.

The Case n =1

To illustrate the derivation of the beam equation, we take n =1. We then have, from

(11.85),

−zw

}, (11.94)

and thus

M =a

T =a

(11.95)

where

=−



−h



−h

(s −h), b

=−



−h



(11.96)

and we take A to be constant. Thus V and w satisfy the equations

∂

∂x

]=−

γβN

∂

∂x



γβNhh

−



−h





∂

∂x

]=−γβNh

−



−h



(11.97)

together with the approximate kinematic conditions (if m

∗

=0)

≈h

≈μw, (11.98)

whence the depth H =s −h is a function only of space: for simplicity we take it to

be a constant. The consequent hydrostatic condition is, using (11.82)

and (11.83),

=μN

. (11.99)

Integrating this, we have

s =s

+μ



N −N

(t)



, (11.100)

where s =s

(constant, since w =0 there) and N =N

at x =x

. Differentiating,

we ﬁnally have

−

+w. (11.101)

766 11 Jökulhlaups

11.6.2 The Beam Boundary Layer

To complete the prescription for N , we need to ﬁnd w in terms of N by solving

the beam equations (11.97). The product γβ = N

/τ

is the ratio of the channel

effective pressure scale to the beam stress, and is generally of O(1) or larger. It is

therefore reasonable to suppose that γβ ν

. In this case, the solution becomes of

singular perturbation type, and supports boundary layers near the lake margins.

Note ﬁrst that since s

∼μ ∼γ ,theterm

γβN

is the largest of those on the right

hand side of (11.97). A leading order approximation (the outer solution) to (11.97)

is then simply

N ≈0, (11.102)

whence also

≈0, (11.103)

and thus

≈T, (11.104)

where T(t)is the tension. From this we have

(s +h)w

. (11.105)

Since s

≈h

≈0, then a

is approximately constant, and so

V ≈

(s +h)w

≈

(11.106)

since w ≈˙s/μ is independent of x. It follows from (11.101) that in fact

w ≈−

, (11.107)

and this together with (11.74) implies that the boundary condition for the channel is

−Q

−

−x

−

, (11.108)

which is essentially the same as we assumed earlier (cf. (11.29)). This relation is

essentially the punch line of our analysis, at least for describing subglacial ﬂoods.

The subsidence of the surface requires further discussion.

The outer solution for w does not satisfy the boundary conditions at the margins,

so that a boundary layer is needed near both margins. Both are similar, and we treat

that at x

. The choice of boundary layer scale depends on the size of the parameters.

We deﬁne

σ =

1/2

(γβ)

1/4

,ω=ν



γβ, (11.109)

and we will assume, as is likely, that ω  1. The boundary layer variables X and v

are deﬁned by

x =x

−σX, V =

, (11.110)

11.6 Cauldrons and Calderas 767

and Eqs. (11.97), (11.98) and (11.99) become

∂

∂X

[−a

]=−N +ω

∂

∂X



Nhh

−

2γ



−h





2[−a

]=−ω



2γ



−h





, (11.111)

=μw, s

=μN

Hard Lake Ice Boundary Conditions

For simplicity, suppose that ω  1, and consider ﬁrst the ice sheet type boundary

conditions for hard overlake ice given by (11.93); these require

=w =0atX =0, (11.112)

i.e.,

=w =0atX =0, (11.113)

and the matching conditions to the outer solution are

v →σ





,w→−

as X →∞. (11.114)

Neglecting terms of O(ω), we ﬁnd that

=−

(s +h)w

XXXX

=−N,

(11.115)

where the ﬂexural viscosity D

is given by

,H=s −h; (11.116)

the second equation in (11.115) is the viscous beam equation.

(11.100) is still the integral of the hydrostatic equation, and thus s

= μw =

μ(N

−

), whence

w =N

−

; (11.117)

differentiating (11.115)

, we derive the viscous beam equation for w in the form

XXXXt

=−w −

, (11.118)

with the boundary conditions for w in (11.113) and (11.114). Exactly the same

boundary layer description applies at the left margin x

−

768 11 Jökulhlaups

Soft Lake Ice Boundary Conditions

Suppose instead that we take the boundary conditions (11.91) which are appropri-

ate for soft overlake ice, such as might be the case in the Vatnajökull cauldrons.

Neglecting O(ω), a ﬁrst integral of (11.111)

−a

, (11.119)

where T

(t) is an integration function. Using the deﬁnitions of a

and b

and the

kinematic conditions, this can be integrated subject to v =0atX =0 to obtain

v =−

(s +h)w



dX +T

X. (11.120)

Applying the matching conditions at X →∞, we ﬁnd T

=0 and







∞

2μ

dX. (11.121)

An entirely similar analysis at the other margin x

−

yields the comparable condition



−



=−



∞

2μ

dX, (11.122)

whence it follows (since each boundary solution is the same) that

σ(x

−x

−



∞

dX, (11.123)

and this determines the tension T . The beam equation (11.118)forw is derived the

same way, the only difference being that the boundary conditions are

w =w

=0atX =0,w→−

as X →∞. (11.124)

A Uniform Approximation for n =1

A uniform approximation to the solution can evidently be made by retaining only

the term in N in (11.97), or equivalently by formally assuming that γβ ∼ ν

.We

make this formal assumption, and additionally allow n = 1, thus we revert to the

equations in the form (11.88).

For the case of hard overlake ice, we have T =0 everywhere, and thus (11.85)

and (11.86) imply by inspection that

(s +h)w

, (11.125)

if we take A to be constant. We then calculate

M =

(n−1)/n

, (11.126)

11.6 Cauldrons and Calderas 769

where

(2n+1)/n

(2n +1)A

1/n

, (11.127)

and the beam equation takes the form

∂

∂x



(n−1)/n



=−

γβN

, (11.128)

subject to

w =w

=0atx =x

, (11.129)

and N is given by (11.101).

The case of soft overlake ice is somewhat impassable: the beam equation for

n =1 becomes non-local. If μ is small, then so is T , and the derivation for hard ice

works in the same way, leading to the same approximate equation (11.128).

11.6.3 Similarity Solutions

We would like to solve the beam equation (11.128) together with the buoyancy

condition (11.101), in order to trace the subsidence of a cauldron. A simple way to

do this is motivated by Nye’s pioneering paper of 1976 on jökulhlaups, where he

showed that the rising limb of the 1972 Grímsvötn ﬂood hydrograph could be well

ﬁtted by a power law in the form

Q ∝

(

t −t)

,S∝

(

t −t)

, (11.130)

which arises from the ﬂood model (11.42) on neglecting the closure term and taking

Q ∼S

4/3

. Since the lake reﬁlling equation implies

≈Q, this approximation also

implies N

∝

(

t−t)

, and this observation suggests how a useful similarity solution

to the present beam boundary layer model can be obtained, if we suppose that the

rising stage of a ﬂood can be described by a water ﬂux varying as in (11.130).

We consider the nonlinear model (11.128) and (11.101) together with the bound-

ary conditions in (11.129). The boundary layer model near x =x

for this equation

is obtained by writing

x =x

−σ

X, σ



γβ



n/[2(n+1)]

. (11.131)

In addition we deﬁne

w =−

+Υ. (11.132)

The bending moment is thus

M =

2/n

|Υ

(n−1)/n

, (11.133)