. 5
( 16)






2891 2740

6371 5153
5153 2891
2 4 6 8 10 12 14 16 18

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)

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)

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
ρ + ∇ · (ρ0v) = 0 (6.4.4)
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)

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)

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

The linearised equation of motion for the displacement u takes the form
‚2 ‚1
ρ ui = σ, (6.4.11)
‚xj ij

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)

and thus the equation of motion can be written as
‚2 ‚ ‚0
(σ1)ij ’ uk
ρ 2 ui = σ. (6.4.13)
‚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
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 .
Cijkl has the symmetries

Cijkl = Cjikl = Cijlk = Cklij. (6.5.3)

σ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)
+ 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

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

which is independent of the hydrostatic pressure.
In general cijkl will be anisotropic, but we can isolate the isotropic part by

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

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)

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
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.
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
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)
‚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)

+ ρgj = ρfj, (7.1.1)
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

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)
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)
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)
‚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
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)
(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
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
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)
We can now simplify the equation (7.2.3) to yield
‚ 1 ·
∇ 2v + 1 ∇ (∇ · v ) .
v + (v · ∇ )v = ’ ∇ p + (7.2.4)
‚t ρ ρUL
We introduce the Reynolds number
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)
‚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

L = 3—106 m, both the depth of the mantle and the size of an average tectonic
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
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:
+ ν∇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
ρ + ∇ · v = 0. (7.2.10)
ρ Dt
We can regard the medium as approximately incompressible if
|∇ · v| , (7.2.11)
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
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
The P´ clet number provides a measure of advective relative to diffusive processes
in the energy equation, and is de¬ned as
Pe = . (7.2.20)
With our set of mantle parameters, the estimate for the P´ clet number Pe is of order
4 . This value indicates that advective processes in the mantle dominate thermal
diffusion processes by four orders of magnitude, outside of thermal boundary

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
4 ‚p 3
Q= dx vz = ’ h, (7.3.7)
3· ‚z

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
= ’(ρs ’ ρl)g, (7.3.8)
so that the ¬‚ow rate would be 4 (ρs ’ ρl)gh3/·.
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
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],
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
ρguz = ’p + 2· at z = 0. (7.4.15)
Now, at z = 0,
‚2vx ‚2vx
= ’2·Bk3 sin(kx),
=· + (7.4.16)
‚x2 ‚z2
and so
p = 2·Bk2 cos(kx) on z = 0, (7.4.17)
= 0 on z = 0. (7.4.18)
The boundary condition at z = 0 is thus
ρguz = ’2·Bk2 cos(kx). (7.4.19)
The surface displacement is to be found from
= vz(x, 0) = Bk cos(kx), (7.4.20)
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
Nu = , (7.5.2)
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 .
ρ = ρ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)
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
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
‚2 1 ‚2 1
T + v1 = κH T + 2T , (7.5.11)
‚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
‚2 1 ‚2 1
’ P +· v+ v = 0, (7.5.12)
‚x2 x ‚z2 x
‚2 1 ‚2 1
‚1 1
’ P ’ ρ0g±th T + · v+ v = 0.
‚x2 z ‚z2 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)

We can simplify the subsequent analysis by taking free surface conditions at
z = 0, L so that the shear stress vanishes and
= 0 at z = 0, L. (7.5.14)
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 ,
‚t ‚x ‚z (7.5.16)
∇4ψ ’ ρ0g±th T = 0.
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,
T 1 = Ta sin(πz/L) sin(kx) eβt,
124 Fluid Flow






0 1 2 3 4 5 6 7

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
β + κH k + κH 2 Ta = ’kθψa,
· k + 2 ψa = ’kρ0g±th Ta.

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 + π2) ,
β= Ra 2 2 (7.5.19)
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)
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

Temperature (Kelvin)
800 800 1200
80 0 1200

1200 1600
120 12















80 00


1600 1600



200 12
0 0
1600 20 0 200 2000
1600 2400 00
2400 201600 2400
2000 00
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
Temperature (Kelvin)

1600 1600 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

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

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

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)

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)
|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
U2 U
©U, ν2 ©U, (7.6.8)
and thus
U ν
Ro = 1, Ek = 1. (7.6.9)

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
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


. 5
( 16)