660

760

410

660

4

2891 2740

2

6371 5153

5153 2891

0

2 4 6 8 10 12 14 16 18

K/p

Figure 6.3. The relation of the shear modulus G and bulk modulus K scaled by the pressure

p for the AK135 model. The depths corresponding to the different segments are indicated.

6.4 Incremental stress and strain

Small deviations from a pre-stressed or pre-strained state can be considered by a

perturbation treatment in either the Lagrangian or the Eulerian viewpoint.

For a state variable q, the perturbation can be written as

qE(x, t) = q0(x) + q1 (x, t), (6.4.1)

E

qL(ξ, t) = q0(ξ) + q1(ξ, t),

ξ ξ Lξ (6.4.2)

as a result of a small change in state from ξ, where q0 is the value in the initial

state. As in our treatment of linearised elasticity (Section 5.1), we can ignore

the distinction between x and ξ in the discussion of the perturbation q1. For an

incremental displacement u = x ’ ξ, the relation between the perturbations in the

material and spatial frames is

q1 = q1 + u · ∇q0, (6.4.3)

L E

where the advective term u · ∇q0 re¬‚ects displacements through the spatial

106 Continua under Pressure

gradients in the initial state. A consequence of (6.4.3) is that the linearised

conservation equation for mass transport can be written as

‚1

ρ + ∇ · (ρ0v) = 0 (6.4.4)

‚t

where v is the velocity ¬eld.

In such a ¬rst-order perturbation, the volume element after deformation

dV 1 = (1 + ∇ · u)dV 0, (6.4.5)

and similarly we can represent the changes in an areal element and its inclination,

through the unit normal n, as

dσ1 = (1 + ∇Σ · u)dσ0, n1 = n0 ’ ∇Σu · n0. (6.4.6)

6.4.1 Perturbations in stress

The Cauchy stress is de¬ned in the spatial frame, so that we can express the effect

of a small perturbation as

σ(x, t) = σ0(x, t) + σ1(x, t), (6.4.7)

and if we refer this stress perturbation back to the material frame

σ1 = σ1 + u · ∇σ0. (6.4.8)

L

For the Piola“Kirchhoff stress tensors in the material frame, introduced in Section

3.5, the perturbations may be obtained by linearising the deformation gradient

relations and so we ¬nd

σPK1 = σ1 + σ0(∇ · u) ’ (∇u)T · σ0, (6.4.9)

L

σSK1 = σ1 + σ0(∇ · u) ’ (∇u)T · σ0 ’ σ0 · ∇u. (6.4.10)

L

The linearised equation of motion for the displacement u takes the form

‚2 ‚1

ρ ui = σ, (6.4.11)

‚xj ij

‚t2

since the initial stress ¬eld σ0 is supposed to be in equilibrium. In terms of the

Lagrangian stress perturbation

σ1 = σ1 ’ u · ∇σ0, (6.4.12)

L

and thus the equation of motion can be written as

‚2 ‚ ‚0

(σ1)ij ’ uk

ρ 2 ui = σ. (6.4.13)

L

‚xk ij

‚t ‚xj

6.5 Elasticity under pressure 107

6.4.2 Perturbations in boundary conditions

The condition that there should be no slip between two solids requires that the

displacement u in any deformation must be continuous across the boundary

[u]+ = 0. (6.4.14)

’

In the case of a ¬‚uid“solid boundary, tangential slip along the boundary is allowed,

and the boundary condition [n0 · u]+ = 0.

’

For stress, the condition of continuity of traction [n·σ]+ = 0, can be transformed

’

0 · σPK]+ = 0. When we take account of the distortion

into the material state as [n ’

of the boundary in the incremental deformation the boundary condition for the

linearised state, including allowance for displacement along the boundary, can be

written as

+

n0 · σPK1 ’ ∇Σ(un · σ0) = 0. (6.4.15)

’

The displacement term is needed in descriptions of seismic sources in terms of slip

on a fault, but otherwise is not signi¬cant. For a contact between a ¬‚uid and a solid,

the traction has to be normal to the distorted ¬‚uid“solid boundary

(n · σ) · (I ’ nn) = 0, i.e., njσij(δik ’ nink) = 0, (6.4.16)

and the equivalent linearised form is

+

n0 · σPK1 + n0∇Σ · (p0u) ’ p0(∇Σ · un0 · (I ’ nn) = 0. (6.4.17)

’

where p0 = ’n0 · σ0 · n0 = n0σ0 n0 is the initial pressure on the ¬‚uid side.

i ij j

6.5 Elasticity under pressure

As we have noted above we can equate the strain energy W for an elastic material

to the speci¬c Helmholtz free energy F/ρ because the thermal in¬‚uences of

deformation will be very small. In a pre-stressed state we can still use the

constitutive equation (4.3.11)

1 ‚W T

σ= F F, (6.5.1)

det F ‚E

in terms of the derivatives of the strain energy with respect to the Lagrangian

(Green) strain E. For ¬nite deformation we need to recognise the differences

in orientation between the Lagrangian and Eulerian triads re¬‚ected in the tensor

transformation. We recall that the quantity σSK = ‚W/‚E de¬ned in the

Lagrangian frame is the second Piola“Kirchhoff stress tensor (cf. Section 3.5).

Now make a quadratic approximation for W(E) about a pre-stressed reference

state

W(E) = W0 + σo Eij + 1 EijCijklEkl, (6.5.2)

ij 2

108 Continua under Pressure

The coef¬cient of the strain tensor itself is just the stress in the reference state σo .

ij

Cijkl has the symmetries

Cijkl = Cjikl = Cijlk = Cklij. (6.5.3)

Thus

1

σSK = σo + CijklEkl, FσSKFT .

σ= (6.5.4)

ij ij

det F

We now examine the behaviour for a small increment in strain superimposed on

an initial stress to de¬ne the elastic moduli for the pre-stressed state. We work with

the ¬rst-order incremental strain µ and rotation ω,

‚ui ‚uj ‚ui ‚uj

1 1

µij = + , ωij = ’ , (6.5.5)

2 2

‚xj ‚xi ‚xj ‚xi

and seek to ¬nd a representation for the incremental stress which is as close as

possible to that for linearised elasticity from a zero-stress state, but with explicit

representations of the in¬‚uence of the hydrostatic pressure and deviatoric stress. A

suitable choice is provided by setting

Cijkl = cijkl + 1 (σ0 δkl + σ0 δij ’ σ0 δjl ’ σ0 δil ’ σ0 δjk ’ σ0 δik), (6.5.6)

2 ij kl ik jk il jl

so that cijkl has the usual i ” j, k ” l, ij ” kl symmetries. The incremental

stress σSK1 is then related to the incremental strain µ by

σSK1 = cijklµkl + p0(2µij ’ µkkδij)

ij

+ 1 („0 µkk + „0 µklδij) ’ µik„0 ’ µkj„0 , (6.5.7)

2 ij kl kj ik

where we have separated the terms involving the initial pressure po and the

deviatoric stress „o. The corresponding incremental Lagrangian stress σ1 is given

L

by

σL1 = cijklµkl + 1 (’„0 µkk + „0 µklδij) + ωik„0 ’ ωkj„0 , (6.5.8)

ij ij kl kj ik

2

which is independent of the hydrostatic pressure.

In general cijkl will be anisotropic, but we can isolate the isotropic part by

writing

cijkl = (κ ’ 2 μ)δijδkl + μ(δikδjl + δilδjk) + γijkl, (6.5.9)

3

where the incompressibility K = κ and rigidity G = μ are de¬ned by

κ = 1 ciijj, 1 1

μ= 10 (cijij ’ 3 ciijj) so that γiijj = γijij = 0. (6.5.10)

9

Then cijkl ’ γijkl is the isotropic tensor which is the best approximation to cijkl in

a least-squares sense, and γijkl is the purely anisotropic component.

6.5 Elasticity under pressure 109

Seismic wavespeeds

As in (5.4.2) we take a plane wave representation

ui = Ui exp[iωpnkxk ’ iωt], (6.5.11)

and insert this into the linearised equation of motion (6.4.13). For high frequencies,

when we can neglect any gradients of the initial stress distribution, we ¬nd

p2 cijklnjnl ’ 1 („0 nlni ’ „0 njnk) ’ ρδik Uk = 0. (6.5.12)

2 kl ij

This is an eigenvalue equation, and as in (5.4.4) there are three roots for p2 for

each direction n, associated with polarisations speci¬ed by the eigenvectors Ur(n).

The equation (6.5.12) for the wave slownesses does not depend explicitly on the

initial hydrostatic pressure, but only on the deviatoric stresses „0. The plane wave

propagation is anisotropic because the slownesses pr and polarisations Ur depend

on the direction n. The anisotropy comes in part from the intrinsic anisotropy

γijkl (6.5.9), and the non-hydrostatic initial stress. In most situations in the Earth,

μ, ||„0||

||γ|| μ and then the independent wave-types will be a quasi-P wave

with wavespeed ± ≈ [(κ + 2 μ)/ρ]1/2, and two quasi-S waves with wavespeeds

3

close to β = [μ/ρ]1/2.

Observations of seismic wave propagation suggest that anisotropic effects are

concentrated in the lithosphere where non-hydrostatic stresses can be signi¬cant,

and in the boundary layers near the 660 km discontinuity and in the complex D

region at the base of the mantle.

7

Fluid Flow

The deformation of ¬‚uids plays an important role in geophysical phenomena, both

in the long term ¬‚ow of the Earth™s mantle and in the more rapid motions in the core.

We consider here the development of the equations describing the velocity ¬eld in a

¬‚uid and introduce the very important concept of non-dimensionalisation. In many

situations with relatively slow motion we can neglect the effects of compressibility

(at least to a ¬rst approximation). As a result we both simplify the equation system

and eliminate the rapidly travelling acoustic waves. We illustrate the development

of some simple ¬‚ows, including simple considerations of the onset of thermal

convection. Finally we consider the effects of rotation which are of sign¬cance

for the rapid ¬‚ows in the Earth™s core and, of course, for oceanic and atmospheric

phenomena.

We will make the assumption of a Stokes ¬‚uid, so that the stress σ is related to

the ¬‚uid velocity v by

‚vj ‚vi ‚vk

’ 2 ·δij

σij = ’pδij + · + , (7.1)

3

‚xi ‚xj ‚xk

where p is the hydrostatic pressure and · the viscosity. The use of (7.1) assumes

that there is no contribution from bulk viscosity.

7.1 The Navier“Stokes equation

The equation of motion for the ¬‚uid in the presence of body force g and

mass-acceleration f is (3.2.4)

‚σij

+ ρgj = ρfj, (7.1.1)

‚xi

and we can relate the acceleration f to the velocity ¬eld v using (2.4.7)

D ‚

f= v = v + (v · ∇)v, (7.1.2)

Dt ‚t

110

7.1 The Navier“Stokes equation 111

that can be written in component terms as,

D ‚ ‚

fj = vj = vj + vk vj. (7.1.3)

Dt ‚t ‚xk

Now, on combining the constitutive equation (7.1) with the equation of motion

(7.1.1) we require

‚vj ‚vi

D ‚ ‚ ‚vk

’ 2 ·δij

ρ vj = ρgj ’ p+ · + , (7.1.4)

3

Dt ‚xj ‚xi ‚xi ‚xj ‚xk

which is usually called the Navier“Stokes equation for the ¬‚uid.

If the viscosity · can be taken as constant over the ¬‚uid, then we can express the

Navier“Stokes equation as

‚2 ‚2

‚ ‚ ‚ 1

ρ vj + ρ vk vj = ρgj ’ p+· vj + 3 · vk, (7.1.5)

‚t ‚xk ‚xj ‚xi‚xi ‚xj‚xk

or in vector form as

‚

ρ v + ρ(v · ∇)v = ρg ’ ∇p + ·∇2v + 1 ·∇(∇ · v), (7.1.6)

3

‚t

which has explicit non-linearity through the convected derivative term.

If the viscosity varies with position (as may for example be associated with the

effect of temperature) we have to include an additional term,

‚vj ‚vi

‚ ‚vk

’ 2 δij

+ · + , (7.1.7)

3

‚xi ‚xi ‚xj ‚xk

on the right-hand side of (7.1.5). We will also need to consider the equations for the

distribution of temperature or any other factor which may in¬‚uence the distribution

of viscosity.

7.1.1 Heat ¬‚ow

We can link the equations of heat ¬‚ow and ¬‚uid ¬‚ow by returning to the equation

for the rate of change of the internal energy U, (4.2.11),

D ‚ ‚ ‚

ρ U=h+ k T + σij vi, (7.1.8)

Dt ‚xk ‚xk ‚xj

where h is the rate of internal heat production, T is the temperature and k is the

thermal conductivity. From the ¬rst law of thermodynamics

D D D

U= Q+ W, (7.1.9)

Dt Dt Dt

where Q represents the thermal contribution to the internal energy and W the

external work. We isolate the contribution from σij‚vi/‚xj in (7.1.8) due to the

pressure term

‚ pD

p vk = ’ ρ, (7.1.10)

‚xk ρ Dt

112 Fluid Flow

using the continuity equation (conservation of mass, equation 2.5.7). We will

assume that the remaining contribution to the internal energy can be represented

in terms of a thermal capacity C so that

D D pD

ρ U = ρ (CT ) ’ ρ (7.1.11)

Dt Dt ρ Dt

“ the particular form of the heat capacity C to be employed will depend on the

thermodynamic state of the material.

We recast the constitutive equation for the viscous ¬‚uid (7.1) as

^

σij = ’pδij + 2·Dij, (7.1.12)

in terms of the deviatoric strain rate

‚ ‚ ‚

^ 1

vi ’ 1

Dij = vj + 3 ‚x vkδij, (7.1.13)

2 ‚xi ‚xj k

which is constructed to have zero trace.

The stress-power contribution to the rate of change of internal energy is

‚ ‚ ^ ‚ vj,

σij vj = ’p vk + 2·Dij

‚xi ‚xk ‚xi

‚ ^^

= ’p vk + 2·DijDij (7.1.14)

‚xk

(since the deviatoric term vanishes when the indices i and j coincide). We have

already explicitly accounted for the pressure contribution in (7.1.11) so that we can

recast (7.1.8) into an equation for the temperature ¬eld in the presence of ¬‚uid ¬‚ow

‚ ‚ ‚ ‚ ^^

ρ (CT ) + vk (CT ) = h + k T + 2·DijDij, (7.1.15)

‚t ‚xk ‚xk ‚xk

which simpli¬es somewhat if C, k are constant.

7.1.2 The Prandtl number

A convenient measure of the ¬‚ow properties of the ¬‚uid is provided by the

kinematic viscosity

ν = ·/ρ, (7.1.16)

which can be regarded as a diffusion coef¬cient for linear momentum [with units

m2 s’1 ].

The equivalent measure for heat transport is provided by the thermal diffusivity

κH = k/ρC (7.1.17)

[with units m2 s’1 ]. The time required for thermal effects to propagate a distance l

is of the order of l2/κH .

7.2 Non-dimensional quantities 113

The relative signi¬cance of momentum and heat transport can be measured by

the ratio of ν to κH which is a dimensionless quantity Pr known as the Prandtl

number

Pr = ν/κH . (7.1.18)

A ¬‚uid with a high Prandtl number has much weaker heat diffusion than momentum

transport, and under this condition we will need to take account of the in¬‚uence of

spatial variations in temperature on the viscosity of the ¬‚uid since the time scales

of ¬‚ow are very long as, e.g., in convection in the mantle of the Earth.

For the Earth™s mantle:

ρ = 4000 kg m’3 , the average density of the upper mantle,

· = 1021 Pa s, the dynamic viscosity of the Earth™s mantle which can be

estimated from studies of post-glacial rebound,

κH = 10’6 m2 s’1 for the thermal diffusivity,

and so the Prandtl number Pr is approximately 2.5—1023 .

7.2 Non-dimensional quantities

A convenient representation of many aspects of ¬‚uid behaviour is provided by

working with dimensionless parameters. We have just introduced the Prandtl

number to describe the ratio of viscosity and thermal diffusivity, and will shortly

encounter the Reynolds number representing the relative signi¬cance of inertial and

viscous forces. In the context of convection, the Rayleigh number measures the

ratio of buoyancy forces to viscous forces and thermal diffusion. In the discussion

of convection it is also useful to look at the Nusselt number that indicates the

signi¬cance of convection in the transport of heat. In Chapter 14 we will also

encounter the P´ clet number providing a measure of the importance of advective

e

processes relative to diffusive processes. In the rotating ¬‚uids the Ekman and

Rossby numbers play an important role in representing the signi¬cance of interial

and viscous forces relative to the Coriolis force.

These non-dimensional parameters are extremely useful. They describe the

essential nature of the ¬‚ow, independently of the particular experimental or

natural setup. The Navier“Stokes equation can be put into a non-dimensional

form, with some of the terms being weighted by these dimensionless parameters.

By inspecting the non-dimensional form of the Navier“Stokes equation, we can

estimate which terms are of leading order and which terms play only a small,

perhaps irrelevant, role.

We de¬ne dimensionless variables, indicated by primes, through

x = Lx , y = Ly , z = Lz , v = Uv . (7.2.1)

Here L is a characteristic length scale, say the radius of the Earth or the depth of

a laboratory tank, and is assumed to be constant. U is a characteristic speed, also

114 Fluid Flow

assumed to be constant. This characteristic speed could be the mean velocity of

a tectonic plate or the average ¬‚ow velocity in a convective tank experiment. The

quantities L, U provide a non-dimensional representation of time t = t L/U. We

also need to transform the derivative operators, so that

‚ U‚ ‚ 1‚

= , = . (7.2.2)

‚t L ‚t ‚x L ‚x

In terms of the non-dimensional variables the Navier“Stokes equation (7.1.6) in

the absence of body forces (g = 0) becomes

U‚ 1

ρ Uv + ρ Uv · ∇ Uv =

L ‚t L

1 1 1

’ ∇ p + · 2 ∇ 2Uv + 1 · 2 ∇ (∇ · Uv ), (7.2.3)

3L

L L

We can now simplify the equation (7.2.3) to yield

‚ 1 ·

∇ 2v + 1 ∇ (∇ · v ) .

v + (v · ∇ )v = ’ ∇ p + (7.2.4)

3

‚t ρ ρUL

We introduce the Reynolds number

ρUL

Re = , (7.2.5)

·

and can then make a further simpli¬cation of the Navier-Stokes equation

‚ 1 1

∇ 2v + 1 ∇ (∇ · v ) .

v + (v · ∇ )v = ’ ∇ p + (7.2.6)

3

‚t ρ Re

We will normally assume that the various quantities are non-dimensionalised, and

henceforth drop the indicative .

The choices of the length scale L and reference speed U are guided by the nature

of the problem under consideration, but different authors may well make dissimilar

choices for the same problem. Thus, if we are interested in problems for the Earth™s

upper mantle, should we take the scale length to be just the upper part of the mantle

L ∼ 1—103 km or that for the whole mantle L ∼ 3—103 km ? We will change

the non-dimensional variables by a factor of three, without changing the physical

situation; there will also be consequent ¬‚ow-ons to other situations.

As a guide to situations involving ¬‚ow in Earth™s mantle we will make a

consistent set of choices throughout our development. The length scale is taken

from the depth extent of the mantle, and the characteristic ¬‚ow speed from plate

velocities:

L = 3—106 m, both the depth of the mantle and the size of an average tectonic

plate,

U = 6—10’10 m s’1 , based on the speed of an average tectonic plate of

20 mm yr’1 , since 1 year is approximately 3—107 s.

7.2 Non-dimensional quantities 115

7.2.1 The Reynolds number

The Reynolds number

ρUL UL

Re = = (7.2.7)

· ν

measures the ratio of inertial forces (which have the characteristic quantity of

ρL2U2) to viscous forces (which have the characteristic quantity of ·UL).

With the choices of physical parameters ρ, ·, scale length L and speed U

appropriate to the Earth™s mantle, we obtain an exceedingly small value for the

Reynolds number of Re =10’20 . From this we conclude that inertial forces in the

Earth™s mantle are small compared with viscous forces.

7.2.2 Stokes Flow

In the limit of very small Reynolds number, the Navier“Stokes equation is often

called the Stokes equation.

‚v ∇p

+ ν∇2v.

Re + (v · ∇)v =’ (7.2.8)

‚t ρ

As a result of the non-dimensionalisation the velocity terms can be expected to be

of order unity. The very small value of the Reynolds number preceding the inertial

terms allows us to discard these terms for ¬‚ow in the Earth™s mantle:

∇p

+ ν∇2v.

0≈’ (7.2.9)

ρ

This special class of ¬‚uid motion is termed Stokes ¬‚ow.

7.2.3 Compressibility

The continuity equation for ¬‚uid ¬‚ow (2.5.7) can be written in the form

1D

ρ + ∇ · v = 0. (7.2.10)

ρ Dt

We can regard the medium as approximately incompressible if

U

|∇ · v| , (7.2.11)

L

and thus

1D U

ρ . (7.2.12)

ρ Dt L

116 Fluid Flow

In a homogeneous ¬‚uid we can write the rate of change of the pressure in terms of

the density and the entropy per unit mass S as

D D ‚p D

p = φ2 ρ + S, (7.2.13)

Dt Dt ‚S Dt

ρ

where φ is the acoustic (bulk-sound) wavespeed, so that the condition (7.2.12) can

be written as

1D 1 ‚p D U

p’ S . (7.2.14)

ρφ2 Dt ρφ2 ‚S ρ Dt L

The requirement (7.2.13) will normally only be satis¬ed if each of the two terms

on the left-hand side is separately small compared with U/L. For the ¬rst term we

need

1D U

p , (7.2.15)

K Dt L

so that the changes in density due to pressure are negligible; we have here

introduced the bulk modulus K = ρφ2. For the second term on the left hand side

of (7.2.13), we can use thermodynamic identities, as in Section 6.1.2, to write

1 ‚p ±th T

= , (7.2.16)

ρφ2 ‚S Cp

ρ

where ±th is the coef¬cient of thermal expansion and Cp is the speci¬c heat at

constant pressure. We can then express the requirement for this contribution to be

small compared to U/L as

±th 1 U

V + ∇(k∇T ) , (7.2.17)

Cp ρ L

where V is the rate of dissipation of mechanical energy per unit mass of ¬‚uid due to

shear viscosity and k is the thermal conductivity of the material (see Section 4.2).

Variations in the density of the ¬‚uid due to internal dissipative heating or molecular

conduction of heat must be small:

±th ·U ±th θk

1, 1, (7.2.18)

CpρL CpρLU

where θ is a measure of the magnitude of the temperature differences in the ¬‚uid.

Both conditions (7.2.18) are readily satis¬ed in the mantle and so the primary

control is provided by (7.2.15).

A useful measure of the signi¬cance of compressibility is the Mach number M,

the ratio of the characteristic velocity to the bulk-sound speed

M = U/φ = U K/ρ. (7.2.19)

For the mantle, the bulk-sound speed for elastic waves φ is around 104 m s’1 and

so M is about 10’14 . For such small M the last term on the right-hand side of

7.3 Rectilinear shear ¬‚ow 117

(7.2.6) will be of no signi¬cance, and we can treat the ¬‚uid as if it is incompressible

(∇ · v ≈ 0). This has the great convenience that we can concentrate on the slow

time scales of the main mantle ¬‚ow and do not need to include the rapid ¬‚uctuations

associated with the passage of seismic waves.

7.2.4 The P´ clet number

e

The P´ clet number provides a measure of advective relative to diffusive processes

e

in the energy equation, and is de¬ned as

UL

Pe = . (7.2.20)

κH

With our set of mantle parameters, the estimate for the P´ clet number Pe is of order

e

4 . This value indicates that advective processes in the mantle dominate thermal

10

diffusion processes by four orders of magnitude, outside of thermal boundary

layers.

7.3 Rectilinear shear ¬‚ow

As we have seen for many ¬‚uids we can neglect the effects of compressibility (at

least as a ¬rst approximation) and then the continuity equation reduces to

∇ · v = 0. (7.3.1)

Consider then a simple situation with ¬‚ow in the z-direction induced by a

pressure gradient

vx = vy = 0, vz = vz(x, y, t), (7.3.2)

for which (v · ∇)v = 0 and the Navier“Stokes equations reduce to

‚2 ‚2 ‚2

‚ ‚ ‚

’ p = ’ p = 0, ’ p+· + vz = ρ 2 vz, (7.3.3)

‚x2 ‚z2

‚x ‚y ‚z ‚t

so that there are no pressure variations with x or y and the pressure p(z) is a

function of z alone. In the case of steady ¬‚ow ‚2vz/‚t2 ≡ 0 and so we have

‚2vz ‚2vz 1 ‚2vz

1‚ ‚vz ‚p

· + =· r +2 2 = = const, (7.3.4)

‚x2 ‚z2 r ‚r ‚r r ‚θ ‚z

where we have expressed the equation in both cartesian and cylindrical polar

coordinates. The nature of the solution of (7.6.4) will depend on the ¬‚ow ¬eld

vz(x, y) in the x’y plane.

We will adopt no-slip conditions and assume that a viscous ¬‚uid in contact with

a solid boundary must have the same velocity as the boundary. For a free surface

the traction must vanish and so the associated shear-stress components will be zero.

118 Fluid Flow

One of the simplest con¬gurations we can envisage is a two-dimensional

situation with solid boundaries at x = ±h and no dependence on y. In that case

‚2vz ‚p

· 2= , vz(h) = vz(’h) = 0, (7.3.5)

‚x ‚z

with solution

1 ‚p 2

(x ’ h2).

vz = (7.3.6)

2· ‚z

which is known as Couette ¬‚ow, with a velocity pro¬le in the form of a parabola

symmetric about the centre line (y = 0).

This simple model gives a reasonable representation of ¬‚ow through a volcanic

conduit in a ¬ssure eruption, e.g., Hekla in Iceland. The ¬‚ow rate per unit length of

¬ssure

h

4 ‚p 3

Q= dx vz = ’ h, (7.3.7)

3· ‚z

’h

on the assumption that the pressure gradient is constant.

Where would the pressure gradient come from in such a case? The upward ¬‚ow

of the magma is driven by the natural buoyancy of the lighter magma relative to the

denser surrounding rock. At a depth h the lithostatic pressure in the rock is ρsgh

for rock density ρs and the corresponding hydrostatic pressure for a magma column

is ρlgh for magma density ρl. If the walls of the conduit are free to deform, the

pressure gradient available to drive the magma is

‚p

= ’(ρs ’ ρl)g, (7.3.8)

‚z

so that the ¬‚ow rate would be 4 (ρs ’ ρl)gh3/·.

3

The solution here ignores the out¬‚ow condition at the top of the pipe and assumes

that we do not have to worry about magma solidi¬cation in the dyke. A more

thorough treatment requires a time-dependent solution with allowance for thermal

effects and a better treatment of surface conditions.

7.4 Plane two-dimensional ¬‚ow

In a situation with a very slow viscous ¬‚ow the acceleration terms ρDv/Dt can be

ignored, if we also have incompressibility ∇ · v = 0 as in (7.6.1) we are in the

Stokes ¬‚ow regime (7.2.9) and the Navier-Stokes equations reduce to

’∇P + ·∇2v = 0, (7.4.1)

where P is the pressure produced by ¬‚uid ¬‚ow

P = p ’ ρgz, (7.4.2)

7.4 Plane two-dimensional ¬‚ow 119

allowing for a gravitational body force. Since the divergence of the velocity

vanishes, the pressure P must be a harmonic function, i.e.,

∇2P = 0. (7.4.3)

For a plane two-dimensional ¬‚ow

vx = vx(x, z), vy = 0, vz = vz(x, z), (7.4.4)

we can satisfy the continuity equation by taking

‚ψ ‚ψ

vx = ’ , vz = , (7.4.5)

‚z ‚x

in terms of a stream function ψ(x, z). The curves ψ(x, z) = const represent

streamlines since v · ∇ψ = 0 and ∇ψ will be normal to the surfaces ψ = const.

The x and z components of (7.4.1) take the form

‚P ‚ ‚P ‚

= ’· (∇2ψ), = · (∇2ψ), (7.4.6)

‚x ‚z ‚z ‚z

and we can eliminate P to give

∇2(∇2ψ) = ∇4ψ = 0, (7.4.7)

which is known as the biharmonic equation.

Viscous ¬‚ow model of glacial rebound

We can illustrate the application of the stream function approach for a simple model

of the recovery of continental regions (e.g. Scandinavia) from the loading imposed

by the ice sheets of the last Ice Age, using the response of a viscous ¬‚uid to a

broadly distributed surface load.

We consider an initial vertical displacement

uz = w0 cos(kz), k = 2π/», (7.4.8)

where the maximum displacement is small compared with the wavelength ». Once

we have determined the response to (7.4.8), more general deformation shapes can

be generated by Fourier synthesis.

We note that the vertical velocity is to be found from ‚xψ and so look for a

stream function in the form

ψ = sin(kx) Z(z), (7.4.9)

and then

d4 2

2d

Z + k4Z = 0,

Z ’ 2k (7.4.10)

4 2

dz dz

with solution

Z = Be’kz + kDze’kz, (7.4.11)

120 Fluid Flow

where the growing exponentials have been excluded to give ¬nite values for Z as

z ’ ∞. Thus

ψ = sin(kx) e’kz(B + kDz), (7.4.12)

and the velocity components are

vx = sin(kx) e’kz[kB + (kz ’ 1)kD],

(7.4.13)

vz = k cos(kx) e’kz[B + kDz].

We will force the horizontal component of velocity to be zero at the surface, which

since uz » can be taken as z = 0, and so

D = B. (7.4.14)

The ¬nal constant B can be found by equating the hydrostatic pressure associated

with the topography to the normal stress at the upper boundary

‚vz

ρguz = ’p + 2· at z = 0. (7.4.15)

‚z

Now, at z = 0,

‚2vx ‚2vx

‚p

= ’2·Bk3 sin(kx),

=· + (7.4.16)

‚x2 ‚z2

‚x

and so

p = 2·Bk2 cos(kx) on z = 0, (7.4.17)

and

‚vz

= 0 on z = 0. (7.4.18)

‚z

The boundary condition at z = 0 is thus

ρguz = ’2·Bk2 cos(kx). (7.4.19)

The surface displacement is to be found from

‚uz

= vz(x, 0) = Bk cos(kx), (7.4.20)

‚t

and so the displacement at z = 0 is given by

‚uz ρguz guz

=’ =’ , (7.4.21)

‚t 2·k 2νk

which can be integrated to give

uz = w0 cos(kx) exp(’t/„c), (7.4.22)

where the decay time „c is given by

„c = 2νk/g = 4πν/g», (7.4.23)

7.5 Thermal convection 121

in terms of the kinematic viscosity ν and the wavelength of the initial disturbance

». The decay time is therefore (for long wavelength disturbances) inversely

proportional to the wavelength.

7.5 Thermal convection

When a ¬‚uid is heated there is normally a decrease of density due to thermal

expansion. A ¬‚uid with heating from below or internally with a cool upper surface

will tend to have cool dense ¬‚uid near the upper boundary and lighter hot ¬‚uid at

depth. Such a con¬guration is gravitationally unstable and so the hot ¬‚uid tends to

rise and the cooler ¬‚uid tends to sink to give a thermal convection regime.

7.5.1 The Rayleigh and Nusselt numbers

Two important dimensionless parameters are explicitly related to convection. The

Rayleigh number provides a dimensionless buoyancy ratio through the ratio of

buoyancy forces relative to viscous forces and thermal diffusion:

g±th ”T L3 ρg±th ”T L3

Ra = = . (7.5.1)

νκH ·κH

For the whole mantle using the scaling parameters U, L above and the additional

physical parameters

±th = 2—10’5 K’1 for the thermal expansivity,

g = 10 m s’2 for the gravitational acceleration, and

”T ≈ 1000 K as a representative temperature difference,

the estimate for the thermal Rayleigh number Ra is around 106 . As we shall shortly

see, this value is more than 1000 times the critical value for the onset of convection

and so the Earth™s mantle must be in a state of convection.

The Nusselt number is a normalized heat ¬‚ux, characterising the effect of

convection.

q

Nu = , (7.5.2)

qc

where q is the observed heat ¬‚ux across a convecting ¬‚uid, and qc is the heat ¬‚ux

for the same ¬‚uid that would be conducted in the absence of convection. Estimates

for the Nusselt number for the Earth are in the range 10“20.

7.5.2 The Boussinesq approximation

Because the principal driving force for convection comes from thermally induced

density variations we must include these effects in the buoyancy term of the

momentum equation, i.e., the gravitational force terms, but the variations will be

small enough that they can be neglected elsewhere. This assumption is known as

122 Fluid Flow

the Boussinesq approximation (Valentin Joseph Boussinesq 1842“1929) and leads

to considerable simpli¬cation of the equations for convection.

We may also assume that density variations ”ρ are linearly related to temperature

variations ”T . The constant of proportionality between density and temperature

variations is the coef¬cient of thermal expansivity ±th .

Set

ρ = ρ0 + ρ1, ρ1 ρ0, (7.5.3)

where ρ0 is a reference density and ρ1 is the thermally induced density change

ρ1 = ’ρ0±th (T ’ T0), (7.5.4)

where ±th is the volume coef¬cient of thermal expansion, and T0 the temperature

corresponding to the reference density ρ0. We can then write the Navier-Stokes

equation in the form

’∇P + ρ1g ^z + ·∇2v = 0,

e (7.5.5)

where P accounts for the hydrostatic pressure in the reference state

P = p ’ ρ0gz. (7.5.6)

The appropriate thermal equation for weak ¬‚ows can be derived from (7.1.15) by

neglecting the work done by the ¬‚ow and treating k, ρ, and C as constants to yield

‚

T + (v · ∇)T = κH ∇2T + h, (7.5.7)

‚t

in terms of the thermal diffusivity κH ; h represents the rate of internal heating.

7.5.3 Onset of convection

We will consider a layer of ¬‚uid of thickness L without internal heat sources (i.e.

h = 0) with upper and lower boundaries maintained at ¬xed temperatures Tu, Tl

respectively.

In the absence of convection when the temperature gradient is small, the

temperature equation reduces to

‚2 T = Tu, z = 0,

T = 0 with (7.5.8)

‚z2 T = Tl, z = L.

The temperature pro¬le is then

Tc(z) = Tu + (Tl ’ Tu)z/L = Tu + θz. (7.5.9)

When the lower temperature is raised suf¬ciently, the system will begin to convect

and we can look at the onset of convection by a perturbation analysis about the

non-convecting state in which T = Tc(z) and v = 0.

7.5 Thermal convection 123

Thus, we set

T = Tc(z) + T 1(x, z), v = v1(x, z), (7.5.10)

where T 1, v1 are ¬rst-order perturbations. The temperature equation then takes the

form

‚2 1 ‚2 1

‚1

T + v1 = κH T + 2T , (7.5.11)

z

‚x2

‚t ‚z

since the nonlinear terms v1‚xT 1, v1‚T 1/‚z are second-order quantities and can

x z

be neglected. The two-dimensional ¬‚uid ¬‚ow equations under the assumption of

incompressibility are

‚1 ‚

vx + v1 = 0,

‚z z

‚x

‚2 1 ‚2 1

‚1

’ P +· v+ v = 0, (7.5.12)

‚x2 x ‚z2 x

‚x

‚2 1 ‚2 1

‚1 1

’ P ’ ρ0g±th T + · v+ v = 0.

‚x2 z ‚z2 z

‚z

The boundary conditions are that the surfaces z = 0 and z = L are isothermal and

that there is no vertical ¬‚ow

T 1 = v1 = 0 at z = 0, L. (7.5.13)

z

We can simplify the subsequent analysis by taking free surface conditions at

z = 0, L so that the shear stress vanishes and

‚v1

x

= 0 at z = 0, L. (7.5.14)

‚z

There is no contribution from ‚v1/‚x since v1 = 0 for all x on the boundaries.

z z

As above, we introduce a stream function ψ such that

‚ψ ‚ψ

v1 = ’ , v1 = , (7.5.15)

x z

‚z ‚x

and then we can write (7.5.11), (7.5.12) in the form of two coupled partial

differential equations

‚2 1 ‚2 1

‚1 ‚

T + θ ψ = κH T + 2T ,

‚x2

‚t ‚x ‚z (7.5.16)

‚1

∇4ψ ’ ρ0g±th T = 0.

‚x

The solution of these two linear equations can be found by separation of variables

and to satisfy the boundary conditions (7.5.13), (7.5.14) we can take

ψ = ψa sin(πz/L) sin(kx) eβt,

(7.5.17)

T 1 = Ta sin(πz/L) sin(kx) eβt,

124 Fluid Flow

4000

Racr

3000

Unstable

2000

1000

Stable

0 1 2 3 4 5 6 7

kL

Figure 7.1. Stability condition for onset of convection in the presence of bottom heating in

terms of the critical Rayleigh number.

where the horizontal wavenumber k, the growth factor β and the constants ψa, Ta

are related by the coupled equations

2

2π

2 2

β + κH k + κH 2 Ta = ’kθψa,

L

(7.5.18)

π2

2

· k + 2 ψa = ’kρ0g±th Ta.

L

For β positive, any perturbation away from the reference state will grow in time and

the heated layer will be convectively unstable. For β negative, any perturbation will

decay and so the layer will be stable against convection.

On eliminating ψa and Ta between the two equations (7.5.18), we can express

the growth rate β in the form

k2L2

κH

’ (k2L2 + π2) ,

β= Ra 2 2 (7.5.19)

2)2

L (k L + π

where the dimensionless Rayleigh number

ρ0g±th θL4 ρ0g±th (T ’ T0)L3

Ra = = . (7.5.20)

·κH ·κH

7.5 Thermal convection 125

The growth rate β will be positive, leading to convective instability, if

(k2L2 + π2)3

Ra > Racr = . (7.5.21)

k2L2

If Ra < Racr, β < 0 and the layer will be stable against convection. The

behaviour is determined by the dimensionless wavenumber kL relating the vertical

and horizontal scales of the perturbation. The critical Rayleigh number has an

√

absolute minimum value when kL = π/ 2 and

min(Racr) = 27π4/4 ≈ 657.51. (7.5.22)

If the Rayleigh number is below this value, convection cannot occur and the layer

will remain stable against bottom heating (Figure 7.1).

From the de¬nition of the Rayleigh number (7.5.20) we can see that we

can interpret the condition Ra > Racr in a variety of ways. In order to

allow convection to begin we may require a suf¬cient temperature gradient or a

temperature-dependent viscosity to drop below a certain value.

A similar calculation may be carried out for no-slip boundary conditions or for a

¬‚uid heated from within, but in these cases the critical Rayleigh number has to be

found numerically.

Most studies of time-dependent convection processes have been addressed by

numerical solution of the governing equations, but that does not detract from the

utility of the linearised stability analysis in assessing the behaviour of different

parts of the horizontal wavenumber spectrum of the ¬‚uid ¬‚ow.

7.5.4 Styles of convection

There are considerable differences between the bottom heated case we have just

considered and the situation with purely internal heating that are manifest through

the temperature and velocity patterns. For the same vigour of convection the

Rayleigh number for internal heating needs to be a factor of ten larger than for

bottom heating, which already indicates that this is a less effective process for the

redistribution of ¬‚uid motion.

A comparison of numerical simulations for the bottom heated and purely

internally heated cases is shown in Figure 7.2 with the Rayleigh number chosen to

give similar convective vigour (Ra =106 for bottom heating). The regions shown

in the ¬gure are extracted from larger simulations to avoid any in¬‚uence from the

sides of the numerical domain.

The ¬rst case (Figure 7.2a) is for bottom heating with both top and bottom

surfaces maintained at constant temperature with free-slip boundary conditions.

Concentrated thermal boundary layers form at both the upper and lower surfaces.

The ¬‚ux of heat through the lower boundary for the case of bottom heating leads

to concentrated upwellings of lighter hot material. The hot ¬‚uid cools at the upper

boundary and becomes more dense and so descends, again in a localised fashion

126 Fluid Flow

(a)

Temperature (Kelvin)

1600800

400

1200

800 800 1200

1600

80 0 1200

2800

2400

1200 1600

120 12

0

80

20

00

0

2800

00

20

00

00

00

20

16

20

1600

00

2400

400

800

2400

800

1600

1600

1600

2400

1200

1200

00

2800

16

2800

80 00

800

1600

2400

1600 1600

12

00

2800

0

200 12

0 0

1600 20 0 200 2000

0

28

1600 2400 00

2400 201600 2400

2000 00

2400

2800

60 280 300 320 340 360 380 400 420 440 460 480 500 520 540 560 580 600 620 640 660 680 700 720 740 7

Distance

(b)

Temperature (Kelvin)

1600 1600 1600

1600

0

160

16

1600

00

1600

1600

1600

1600

1600

16

1600

00

16

00

1600

60 280 300 320 340 360 380 400 420 440 460 480 500 520 540 560 580 600 620 640 66 0 680 700 720 740 7

Distance

Figure 7.2. Comparison of convection from numerical simulations displayed through the

temperature distribution between (a) bottom heating and (b) uniform internal heating, with

the Rayleigh numbers chosen to give similar vigour of convection and a quasi-steady state.

Note that the temperature contrasts are both more con¬ned and of greater amplitude in the

bottom heated case.

(Figure 7.2a). The process leads to a set of regular cells outlined by rising hot limbs

and descending cold limbs with very strong temperature contrasts to the interior.

The second case (Figure 7.2b) has the same geometrical con¬guration but with

uniform internal heating with a cold upper boundary, and no heat ¬‚ux through the

lower boundary. In this internally heated case the coooling effect of the upper

boundary creates a thicker cool thermal boundary layer from which concentrated

downwellings emerge as the dense material becomes unstable. The pattern of

convection (Figure 7.2b) shows a number of irregularly spaced zones of descending

cold material, but only very weak hot upwelling. Because of the relatively low

temperature contrasts, the cold material warms noticeably in descent and may not

even reach the lower boundary. Upwelling occurs as a pattern of diffuse ¬‚ow rising

between the cool downwellings; without the injection of bottom heat there is no

driver for the generation of streams of light hot ¬‚uid.

7.6 The effects of rotation

In a rotating reference frame with angular velocity ω there are accelerations

associated with the use of a non-inertial reference frame. These forces are

7.6 The effects of rotation 127

commonly known as the centrifugal and Coriolis forces. The acceleration

Dv Dv

= + © — (© — x) + 2© — v. (7.6.1)

Dt Dt

I R

Here the subscript I refers to the inertial reference frame, whereas subscript R refers

to the rotating frame of reference. The term on the left-hand side, (Dv/Dt)I, is the

actual acceleration experienced by a ¬‚uid particle. The ¬rst term on the right-hand

side, (Dv/Dt)R, is the acceleration relative to the rotating reference frame. The

second and third terms on the right-hand side are, respectively, the centrifugal and

Coriolis forces.

We can expand the acceleration in the rotating reference frame in the usual way:

Dv ‚v

= + (v · ∇)v. (7.6.2)

Dt ‚t

R

We can drop the subscript R, because all velocities refer to the rotating reference

frame for the rest of this chapter. The equation of motion for a ¬‚uid in a rotating

reference frame takes the form

‚v ∇p

’ © — (© — x) ’ 2© — v + ν∇2v,

+ (v · ∇)v = ’ (7.6.3)

‚t ρ

where ν = ·/ρ is the kinematic viscosity.

7.6.1 Rapid rotation

In many cases the centrifugal force is not important. This arises because we can

express it as the gradient of a scalar quantity:

© — (© — x) = ’∇( 1 ©2d2), (7.6.4)

2

where the quantity d is the distance from the axis of rotation. Because of this

relation, we can absorb the centrifugal forces into the pressure term.

We replace the pressure term by an augmented pressure p = (p ’ 1 ρ©2d2). 2

This transformation reduces the problem to one that is identical in physical context,

but the centrifugal force does not appear. Note that this procedure is analogous to

the procedure of subtracting out the hydrostatic pressure to remove the effect of

gravitational forces.

The centrifugal force is balanced by a radial pressure gradient which is present

whether or not there is any ¬‚ow relative to the rotating reference frame, and which

does not interact with any such ¬‚ow.

The limitation to this transformation is the same as for the gravitational case.

First, the pressure must not appear in the boundary conditions. Second, since we

have moved ρ inside the differential operator ∇, the density must be constant.

Centrifugal force variations associated with density variations will give rise to body

forces that can alter or even cause a ¬‚ow.

128 Fluid Flow

The modi¬ed equation of motion for a rapidly rotating ¬‚uid is then

‚v ∇p

’ 2© — v + ν∇2v.

+ (v · ∇)v = ’ (7.6.5)

‚t ρ

7.6.2 The Rossby and Ekman numbers

In many geophysical ¬‚ows the effects of the Coriolis force are important to the

degree that they dominate the ¬‚ow. Let us take this situation to the extreme and

assume that the Coriolis force is large compared with the effects of inertia and

viscous forces. For steady ¬‚ow this means

|v · ∇v| |© — v|, (7.6.6)

and

|v · ∇2v| |© — v|. (7.6.7)

These assumptions imply that the following relations must hold for a

representative velocity (U), the rotational rate (©) and the length (L) scale in the

¬‚uid:

U2 U

©U, ν2 ©U, (7.6.8)

L L

and thus

U ν

Ro = 1, Ek = 1. (7.6.9)

©L2

L©

The quantities Ro and Ek are known as the Rossby number and the Ekman number.

The Rossby number represents the ratio of the inertial to the Coriolis forces, and

the Ekman number the ratio of the viscous and Coriolis forces.

Note: Viscous effects may become important, even when these simple scaling

considerations indicate otherwise, for example in boundary layers.

7.6.3 Geostrophic ¬‚ow

When both the Rossby number and the Ekman number are small the equation of

motion becomes

∇p

2© — v = ’ . (7.6.10)

ρ

Here the pressure has been modi¬ed to include the centrifugal pressure.

Flows in which Coriolis and pressure forces are in balance are called geostrophic

¬‚ows. An important property of geostrophic ¬‚ow is immediately evident. The

Coriolis force is always perpendicular to the ¬‚ow direction, as is the pressure. This

means that pressure is constant along a streamline, which is in marked contrast to

¬‚ow in non-rotating systems where pressure varies along a streamline.

7.6 The effects of rotation 129

7.6.4 The Taylor“Proudman theorem

An interesting property of the geostrophic ¬‚ow equation becomes apparent when

we take the curl of the geostrophic balance. We obtain

∇ — (© — v) = 0. (7.6.11)

Using vector identities, we can write

© · ∇v ’ v · ∇© + v(∇ · ©) ’ ©(∇ · v) = 0. (7.6.12)

© is not a function of position, so the second and third terms vanish. The fourth

term is also zero by virtue of the continuity equation for an incompressible ¬‚uid

(∇ · v = 0). We thus have

© · ∇v = 0. (7.6.13)

If © points into the z-direction, we require

‚v ‚v

© = 0, and so = 0; (7.6.14)

‚z ‚z

there is no variation of the velocity ¬eld in the direction parallel to the rotation axis.

This result is known as the Taylor“Proudman theorem.

In component form the Taylor“Proudman theorem is expressed as

‚u ‚v ‚w

= = = 0. (7.6.15)

‚z ‚z ‚z

If the ¬‚uid is contained by boundaries perpendicular to the rotation axis, so that

w = 0 at the boundary, then we have

‚u ‚v

= = 0, w = 0, everywhere in the ¬‚uid. (7.6.16)

‚z ‚z

The ¬‚ow is entirely two-dimensional in planes perpendicular to the axis of rotation.

One striking consequence of the Taylor“Proudman theorem is the formation of

features known as Taylor columns which are expected, for example, in core ¬‚ow.

7.6.5 Ekman layers

We assume a geostrophic ¬‚ow is present in a ¬‚uid rotating around the z-axis

bounded by a wall in the xy-plane. The ¬‚ow is uniform and uni-directional with

component u0 in the x-direction:

1 ‚p0 1 ‚p0

2©u0 = ’ , 0=’ . (7.6.17)

ρ ‚y ρ ‚x

There is a boundary region between the geostrophic ¬‚ow in the interior of the ¬‚uid

and the ¬‚ow just at the wall with boundary conditions u = v = w = 0, where

viscous forces are important. This is because in the boundary region a stress exists

between the boundary and the ¬‚uid. In the absence of other forces the stress can

130 Fluid Flow

be balanced only by the Coriolis force. The Coriolis force acts at right angles to

the motion. Hence we assume that ¬‚ow in the boundary layer is parallel to the wall

with non-zero components u, v,

‚u ‚u ‚v ‚v

= = = = 0. (7.6.18)

‚x ‚y ‚x ‚y