Fowler A. Mathematical Geoscience

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

490 8 Mantle Convection

8.4.1 Boundary Layer Theory

The Boussinesq equations describing thermal convection are written in the follow-

ing dimensionless form:

∇.u =0,

=−∇p +∇

u +RT k,

=∇

(8.132)

where u is velocity, p is pressure, T is temperature, and the Rayleigh and Prandtl

numbers are deﬁned by (8.12) and (8.13); k is the unit vector in the vertical direc-

tion.

By considering only two-dimensional motion in the (x, z) plane, we deﬁne the

stream function ψ by

u =−ψ

,v=ψ

; (8.133)

the vorticity is then (0,ω,0), where ω =−∇

ψ. Taking the curl of the momentum

equation, we derive the set

ω =−∇

ψ,

+ψ

−ψ

=∇

dω

=−RT

+∇

ω,

(8.134)

which are supplemented by the boundary conditions

ψ,ω =0onx =0,a, z=0, 1,

T =

on z =0,T=−

on z =1,

=0onx =0,a;

(8.135)

here a is the aspect ratio, and we have chosen free slip (no stress) conditions at the

cell boundaries.

Rescaling The idea is that when R  1, thermal boundary layers of thickness

δ  1 will form at the edges of the ﬂow, and both ψ and ω will be 1 in the ﬂow.

To scale the equations properly, we rescale the variables as

ψ,ω ∼

, (8.136)

8.4 High Rayleigh Number Convection 491

and deﬁne

δ =R

−1/3

. (8.137)

Rescaled, the equations are thus, in the steady state,

ω =−∇

ψ,

−ψ

=δ

∇

ω =

σδ

dω

(8.138)

In order that the inertia terms be unimportant, we require σδ

1, i.e. σ R

2/3

This assumption is easily vindicated in the Earth’s mantle, but is difﬁcult to achieve

in the laboratory. We also see that if internal heating is included, then we should add

atermδ

H to the right hand side of (8.138). For typical Earth parameters, this is of

order 10

−4

, and in practice negligible.

As in any singular perturbation procedure, we now examine the ﬂow region by

region, introducing special rescalings in regions where boundary conditions cannot

be satisﬁed.

Core Flow The temperature equation is linear in T , and implies T = T

(ψ) +

O(δ

). For a ﬂow with closed streamlines, the Prandtl–Batchelor theorem then im-

plies T

= constant (this follows from the exact integral



∂T

∂n

ds = 0, where the

integral is around a streamline, whence T



(ψ)



∂ψ

∂n

ds =0); it then follows that T

is constant to all (algebraic) orders of δ, and is in fact zero by the symmetry of the

ﬂow. Thus

T =0,

∇

ψ =0,

(8.139)

and clearly the core ﬂow cannot have ψ =ω =0 at the boundaries, for non-zero ψ.

In fact, ω jumps at the side-walls where the plume buoyancy generates a non-zero

vorticity. We examine the plumes next.

Plumes Near x =0, for example, we rescale the variables as

x ∼δ, ψ ∼δ, (8.140)

and denote rescaled variables by capital letters. At leading order, we then have

≈0, (8.141)

whence Ψ ∼v

(z)X, and to match to the core ﬂow, we deﬁne v

=ψ

x=0

as the

core velocity at x =0. Also

−Ψ

≈T

(8.142)

492 8 Mantle Convection

the latter of which integrates to give

ω =



TdX, ω



∞

TdX, (8.143)

where matching requires ω

to be the core vorticity at x =0. Integration of (8.142)

gives



∞

TdΨ=C, (8.144)

where C is constant, and it follows that the core ﬂow must satisfy the boundary

condition ωψ

=C on x =0. In summary, the effective boundary conditions for the

core ﬂow are

ψ =0onx =0,a, and z = 0, 1,

=0onz =0, 1,

=−C on x =0,ψ

=C on x =a,

(8.145)

and the solution can be found as ψ = C

1/2

ψ, where

ψ is determined numerically.

It thus remains to determine C. This requires consideration of the thermal boundary

layers.

Thermal Boundary Layers Near the base, for example, we rescale the variables

z ∼δ, ψ ∼δ, ω ∼δ, (8.146)

to ﬁnd the leading order rescaled equations as

≈0, (8.147)

whence Ψ ∼u

(x)Z, and −u

is the core value of the basal velocity. Then Ω

∼

determines Ω (with Ω =0onZ = 0, and Ω ∼ω

(x)Z as Z →∞, where ω

the core value of the basal vorticity), and T satisﬁes

−Ψ

≈T

. (8.148)

In Von Mises coordinates x,Ψ , the equation is

−T

≈

∂

∂Ψ



∂T

∂Ψ



, (8.149)

and putting ξ =



(x) dx (so ξ marches from right to left in the direction of

ﬂow), this is just the diffusion equation

ΨΨ

, (8.150)

8.4 High Rayleigh Number Convection 493

with

T =

on Ψ =0,T→0asΨ →∞. (8.151)

A quantity of interest is the Nusselt number, deﬁned as

Nu=−



∂T

∂z

(x, 0)dx, (8.152)

and from the above, this can be written as

Nu≈





∞

T |

z=0

dΨ



x=0

x=a

1/3

. (8.153)

Notice that the plume temperature equation can also be written as (8.150), where ξ

is extended as



(z) dz,etc.

Corner Flow The core ﬂow has a singularity in each corner, where (if r is dis-

tance from the corner), then ψ ∼ r

3/2

, ω ∼ r

−1/2

, and (for the corner at x = 0,

z =0, for example) x,z ∼r. There must be a region where this singularity is allevi-

ated by the incorporation of the buoyancy term. This requires ω/r

∼1/δr, whence

r ∼ δ

2/3

. Rescaling the variables as indicated (x,z ∼ δ

2/3

, ψ ∼δ, ω ∼δ

−1/3

) then

gives the temperature equation as

−Ψ

≈δ∇

T, (8.154)

which shows that (since the ψ scale, δ, is the same as that of the boundary layers

adjoining the corner) the boundary layer temperature ﬁeld is carried through the

corner region. The corner ﬂow has T ∼T(Ψ), so that

∇

Ψ +T



(Ψ )Ψ

=0, (8.155)

with appropriate matching conditions. Jimenez and Zuﬁria (1987) claim that the

equivalent problem for the case of no-slip boundary conditions has no solution, but

do not adduce details. Their inference is that the boundary layer approximation fails:

this seems a hazardous conclusion.

Solution Strategy The temperature equation (8.150) must now be solved in the

four regions corresponding to the boundary layer at z =0, plume at x =0, boundary

494 8 Mantle Convection

layer at z =1, plume at x =a, with T being continuous at each corner, and

T →0asΨ →∞,

T =

on Ψ =0 [z =0, base],

∂T

∂Ψ

=0onΨ =0 [x =0, left],

T =−

on Ψ =0 [z =1, top],

∂T

∂Ψ

=0onΨ =0 [x =a,right];

(8.156)

in addition, T is periodic in ξ . Beginning from x =a, z =0, denote the values of ξ

at the corners as ξ

(x = 0,z=0), ξ

(x = 0,z=1), ξ

(x = a,z =1).Fromthe

deﬁnition of ξ ,wehaveξ

1/2

, where

are independent of C. Putting

ξ =C

1/2

ξ, Ψ =C

1/4

Ψ, (8.157)

then the problem for T(

ξ,

Ψ)is independent of C. If we can solve this numerically,

then



TdΨ=C

1/4



Ψ , thus

Nu ≈C

1/4





∞



1/3

, (8.158)

and lastly, C is determined from

C =





∞



4/3

, (8.159)

where the integral is evaluated at

. Since also −C

3/4



∞

Ψ at ξ =0, (8.158)

can be written as

Nu ≈2CR

1/3

. (8.160)

The necessary numerical results to compute C are given by Roberts (1979) and

Jimenez and Zuﬁria (1987). The results are slightly different, with the latter pa-

per considering Roberts’s numerical results to be wrong. For a = O(1),wehave

2C ≈0.2. Figure 8.3 shows the typical isotherm proﬁle for high Rayleigh number

constant viscosity convection.

No-slip Boundary Conditions For no-slip boundary conditions, the necessary

preliminary rescaling is ψ ∼1/δ

, ω ∼ 1/δ

, where δ = Ra

−1/5

. Thus the Nusselt

number Nu ∼R

1/5

. There is no longer parity between the thermal boundary layers

and plumes, as the former are slowed down by the no-slip conditions. The rescaled

8.5 Variable Viscosity 495

Fig. 8.3 Temperature

isotherms of a calculation of

constant viscosity convection

at Rayleigh number

0.9 ×10

. The thermal

plumes and boundary layers

are clearly indicated. Figure

courtesy Mike Vynnycky

equations are

ω =−∇

ψ,

−ψ

=δ

∇

ω =

(8.161)

The core ﬂow is as before; the thermal boundary layers have ψ ∼ δ

, ω ∼ 1, z ∼

δ, so that vorticity balances buoyancy (an omission in Roberts’s 1979 paper, thus

precluding his similarity solution), and all three equations are necessary to solve

for T ; it is still the case that



Tdψis conserved at corners, but now in the plume

x ∼ δ

3/2

, ψ ∼ δ

3/2

, and T ∼ δ

1/2

. The initial plume proﬁle is effectively a delta

function, and the plume temperature is just the resultant similarity solution. The

remainder of the structure must be computed numerically, something which has not

been done.

8.5 Variable Viscosity

If we try and apply the above theory to the convection of the Earth’s mantle, we

would predict a surface velocity of order

u ≈C

1/2

2/3

, (8.162)

and with κ

= 10

−6

−1

, d = 3000 km, C

1/2

= 0.3, Ra = 7 × 10

, we ﬁnd

u ≈53 cm y

−1

. This is in remarkably good agreement with observed plate velocities

1–10 cm y

−1

. In fact, when the theory was ﬁrst proposed in 1967 by Turcotte and

496 8 Mantle Convection

Oxburgh, they used a depth d of 700 km, since it was then thought that only the up-

per mantle was of low enough viscosity to convect. Since u ∝d, the corresponding

estimate for u would be 12.5 cm y

−1

: essentially perfect! But even the whole-mantle

convection prediction is very good, given that effects of cell size, sphericity, and the

variability of parameters such as α with depth will all modify the result to some

extent.

However, as we described in the introduction to the chapter, there is a real prob-

lem with this theory: the viscosity of mantle rock is highly variable, so that where it

is cold, the rock is undeformable; and this is precisely in the thermal boundary layer

at the top surface. We now provide a high Rayleigh number boundary layer theory

for this situation. Essentially, for sufﬁciently high Rayleigh number, we will have

rapid boundary layer convection as before, except that this occurs below a stagnant

lid. Just as for no-slip boundary conditions, we will ﬁnd that the rigid lid causes an

1/5

behaviour in the Nusselt number.

8.5.1 Rheology of Polycrystalline Rocks

Many experiments on crystalline rocks lead to an expression for the viscosity of the

form

η =

2Aτ

n−1

exp



∗

+pV

∗



, (8.163)

where T is (absolute) temperature, p is pressure, τ is the second invariant of the

deviatoric stress tensor (2τ

= τ

), and the constants are a rate factor A,the

gas law constant R =8.3Jmol

−1

, the activation energy E

∗

, and the activation

volume V

∗

. Typical values of these constants are

A =10

MPa

−n

−1

,n=3.5,

∗

=533 kJ mol

−1

∗

=1.7 ×10

−5

mol

−1

(8.164)

In writing the equations in dimensionless form, we now have to choose repre-

sentative values of the absolute temperature, in order to have a meaningful viscos-

ity scale. Because the viscosity is so variable, it is not obvious how to do this. It

turns out that the right temperature to choose is the ‘rate-controlling’ value in the

asthenosphere, which is the region just below the (cold, rigid) lithosphere. It is rate-

controlling in the sense that the viscosity is minimal there, so that the velocity scale

is controlled by the asthenospheric viscosity.

However, we do not know the value of the asthenospheric temperature (although

we know what a reasonable value may be, i.e., 1500 K); and even if we did, we

do not know the viscosity as we do not know the appropriate scales for τ . So non-

dimensionalising the equations has to be done ‘blindly’, as it were, with the proper

choice of scales being determined after the fact.

8.5 Variable Viscosity 497

8.5.2 Governing Equations

We consider two-dimensional convection in a Cartesian box. A Boussinesq form of

the governing equations is

∂u

∂x

∂w

∂z

=0,

∂p

∂x

∂τ

∂x

∂τ

∂z

∂p

∂z

∂τ

∂x

−

∂τ

∂z

−ρg,

=2η

∂u

∂x

=η



∂u

∂z

∂w

∂x



=κ∇

(8.165)

Here, τ

(= τ

) and τ

(= τ

) are the longitudinal and shear components of the

deviatoric stress tensor. In addition, the viscosity is deﬁned by (8.163), the second

stress invariant is

=τ

+τ

, (8.166)

and we suppose the density is

ρ =ρ



1 −α(T −T

)



. (8.167)

We have ignored inertia terms, and also have put isothermal compressibility and

internal heating to zero.

8.5.3 Boundary Conditions

At the base, z =0, we prescribe

T =T

,w=0,τ

=0; (8.168)

at the top surface z =d,

T =T

,w=0,τ

=0; (8.169)

and at the sides x =0 and x =ad (say):

∂T

∂x

=0,u=0,τ

=0. (8.170)

In addition, the normal stress should be continuous. In practice this is used to

prescribe the uplift. For convection under a free surface at z =d, say, we prescribe

p +τ

, (8.171)

498 8 Mantle Convection

where p

is the surface loading (zero if atmospheric pressure; or hydrostatic pressure

if the mantle is sub-oceanic), and this extra condition will give  if p

is prescribed.

8.5.4 Boundary Layer Analysis

We begin by non-dimensionalising the terms as follows:

p −ρ

g(d −z),τ

,τ

,τ ∼

≡τ

,T∼T

η ∼η

,t∼

,(x,z)∼d, (u,w)∼

(8.172)

At this point we do not know either T

or η

: they must be determined later.

We introduce a stream function ψ satisfying

u =−

∂ψ

∂z

,w=

∂ψ

∂x

; (8.173)

then the resulting dimensionless equations are these:

∂p

∂x

∂τ

∂x

∂τ

∂z

∂p

∂z

∂τ

∂x

−

∂τ

∂z

−Ra(1 −T),

=−2η

∂

∂x∂z

=η



∂

∂x

−

∂

∂z



∂ψ

∂x

∂T

∂z

−

∂ψ

∂z

∂T

∂x

=∇

=τ

+τ

η =

n−1

exp



1 −T +μ{1 −z +Bp/Ra}

εT



(8.174)

where the parameters are given by

Ra =

αρ

,ε=

∗

μ =

gdV

∗

,B=αT

Λ =

2η

Aτ

n−1

exp



∗



(8.175)

8.5 Variable Viscosity 499

We expect that η

will be roughly the asthenospheric viscosity, and we proceed on

the basis that Ra  1, and also that ε  1, since if we take E

∗

= 533 kJ mol

−1

R = 8.3Jmol

−1

, T

=1500 K, then ε ≈0.023. The Boussinesq number B ≈ 0.05,

so we may neglect the term Bp/Ra in the viscosity. The other parameter μ takes

the approximate value 2.8 if V

∗

= 1.7 × 10

−5

mol

−1

, ρ = 3 × 10

kg m

−3

g =10 m s

−2

, d =3 ×10

m, and is clearly important; however, we will ﬁrst study

the simpler problem in which μ = 0, only adding some comments in the notes

(Sect. 8.8) about the possible structure if μ = O(1). Thus we take the viscosity

to be

η =

n−1

exp



1 −T

εT



. (8.176)

The structure we anticipate is this. There is a cold, rigid lid of thickness νs(x)

(say) adjoining the top surface, in which T<1 and η is exponentially large. Hence

ψ ≈ 0 there and ∇

T = 0. It is formally convenient to suppose ν  1 (with s =

O(1)), so that the temperature proﬁle is approximately linear with depth, and we

shall make this assumption.

Below the lid is a well-stirred, rapidly convecting region, in which T ≈ constant

(and thus T ≈1) and there are thermal boundary layers at the base and beneath the

lid, and plumes at the side. In these layers, T = 1 +O(ε) (otherwise η would be

exponentially small, and ψ exponentially large), and in particular this tells us (with

μ =0) that T

≈T

Suppose that the thermal boundary layer (sometimes called the delamination

layer) which joins the rigid lid to the rapidly convecting core is of dimensionless

thickness δ 1. We write

T =1 +εθ (8.177)

in this region; then continuity of heat ﬂux into the lid implies T

∼ 1/ν ∼ε/δ,so

we choose

ν =δ/ε. (8.178)

Next suppose that the plumes are of thickness δ

(this will also be the thickness

of the basal boundary layer). As for the isoviscous situation where the top condition

is no slip, we anticipate that δ

δ. Since the ﬂow below the lid has T =1 +O(ε)

everywhere, the same scales as for the isoviscous case should apply, and this implies

that in the core p,τ

,τ

,ψ ∼1/δ

, and therefore we choose

Λ =

2(n−1)

, (8.179)

and a balance of shear stress with buoyancy in the plume implies

=Raεθ

, (8.180)

if θ ∼θ

there (note also that ψ ∼1/δ

in the plume).