. 13
( 16)


the concept of any internal density boundaries apart from those associated with the
410 km and 660 km seismic discontinuities.
328 The In¬‚uence of Rheology: Asthenosphere to the Deep Mantle

13.4.2 Material properties
In the deep Earth the effect of pressure can be very signi¬cant and change
perceptions of material properties based on attainable laboratory conditions. Many
experiments and ab initio calculations refer to isothermal conditions rather than the
adiabatic state appropriate to seismic waves. The relation between the adiabatic
bulk modulus KS and the isothermal bulk modulus KT is (6.1.18)
KS = KT (1 + ±th γth T ), (13.4.6)
in terms of the thermal expansion coef¬cient ±th , Gr¨ neisen constant γth and
temperature T . The derivatives with respect to temperature are related by
‚KS ‚KT ‚±th γth
= (1 + ±th γth T ) + KT ±th γth + T . (13.4.7)
‚T ‚T ‚T
The thermal expansivity ±th decreases with increasing pressure but increases
with temperature. The Gr¨ neisen parameter γth will vary only slightly under
lower mantle conditions. The isothermal bulk modulus varies monotonically with
pressure and temperature with ‚KT /‚p > 0 and ‚KT /‚T < 0. At lower mantle
conditions (‚KT /‚T )V ∼ -0.002 GPa K’1 , ±th ∼ 2—10’5 K’1 , γth ∼ 1.5 with KT ∼
500 GPa, and so the second and third terms in (13.4.7) are signi¬cant. Indeed for the
derivatives at constant volume (‚K/‚T )V these extra terms can outweigh the ¬rst
at higher temperatures, and so reverse the sign of the adiabatic derivative compared
with the isothermal. However, the in¬‚uence of pressure on volume is such that this
will not occur for the derivatives at constant pressure (‚K/‚T )P, even though the
adiabatic temperature derivative will be somewhat smaller than the isothermal.

13.4.3 The lower boundary layer
Below a depth of about 2100 km the amplitude of seismic wavespeed heterogeneity
inferred from different styles of seismic tomography shows a steady increase
and near the core“mantle boundary the size of the variations is comparable to
those close to the surface (see, e.g., Figure 11.22). The increase in level of
heterogeneity is common to both the bulk modulus and shear modulus, even
though the major anomalies in the Paci¬c and beneath southern Africa near the
core“mantle boundary have opposite sign for the anomalies in bulk-sound and shear
wavespeeds. This anti-correlation of the different material properties has frequently
been used as an argument for the presence of chemical heterogeneity in the D layer
at the base of the mantle (e.g., Masters et al., 2000), since it cannot be achieved by
purely thermal means.
The discovery of the perovskite to post-perovskite phase transition (Murakami
et al., 2004), with an increase in density, has raised a range of new questions about
the nature of structures at the base of the mantle. The phase transition occurs for
pressures comparable to those at the core“mantle boundary, and it seems likely
that the presence of minor components brings the transition into the mantle regime.
13.4 The deeper mantle 329

The estimates of the Clapeyron slope for the pure Mg end member are of the order
of 8“10 MPa K’1 , more than twice the magnitude, and of opposite sign, to the
γ-spinel to perovskite transition near 660 km depth (γc ∼ ’3 K’1 ). As a result
the in¬‚uence of temperature on the topography of this exothermic transition will
be magni¬ed compared with the de¬‚ections of phase boundaries within subducted
slabs, discussed above in Section 13.3. A further complication comes from the
in¬‚uence of the iron (Fe) content of the perovskite and the admixture of oxide
components. The topographic undulations of the phase boundary are controlled by
the change in the phase transition pressure
‚p ‚p ‚p
δp = δT + δX = γc δT + δX, (13.4.8)
‚T ‚X ‚X
where, as in (13.4.5), X represents some additional chemical component. The
actual con¬guration will depend on the balance between thermal buoyancy and
compositional buoyancy. When the effect of including an additional component is
to increase the transition pressure (‚p/‚X > 0), then temperature and composition
work together and the phase transition moves to greater depth. However, if the
presence of a component is such that ‚p/‚X < 0 the depression of the phase
boundary will be reduced and could even lead to an elevated boundary in a region
of upwelling.
The presence of undulations in the phase boundary associated with mixed
materials leads to localized forces in the momentum equation (8.1.4), due to the
density differences arising from the phase change. An additional constraint has
also to be placed on the continuum equations to ensure conservation of the different
chemical species. The phase boundary can also be expected to involve some change
in viscosity that will in¬‚uence ¬‚ow patterns.
The complex mixture of components at the base of the mantle means that the
simple analysis for a single phase change has to be supplemented by a consideration
of mixture properties. Spera, Yuen & Giles (2006) have looked at the in¬‚uence of
Fe on the perovskite to post-perovskite transition and conclude that there is likely
to be a mixed phase region perched above a post-perovskite layer. For reasonable
molar proportions Fe/(Fe+Mg) ∼ 0.05“0.20, and temperatures appropriate to the
region around the core-mantle boundary (3800“4400 K), a perched layer 200“400
km thick would start about 200“400 km above the core-mantle boundary. The
mixed layer would be displaced upwards for regions with lower temperatures, such
as where former slab material reaches the D zone. The postulated mixed phase
zone occupies the region where seismic heterogeneity increase rapidly with depth
and the base might well be identi¬ed as a seismic discontinuity.
Mantle Convection

The Navier“Stokes equation introduced in Chapter 7 provides the basis for
understanding the behaviour of mantle convection. Because the mantle is highly
viscous, the dominant contributions to the momentum balance come from viscous
and buoyancy forces. This means that a scaling analysis from boundary layer theory
can be used to compute mantle convection velocities from estimates of mantle
viscosity and buoyancy forces. The peculiar geometry, or planform, of mantle
convection cells, with subduction zones separated by distances several times the
mantle depth, results from a strong increase in mantle viscosity with depth. The
internal mantle temperature distribution departs considerably from the adiabat,
because the mantle overturns slowly in a time comparable to the internal heat
production time scale and the rate of secular cooling. The topography associated
with cooling of the oceanic lithosphere and with hot spot swells gives a useful
reference to constrain the ratio of internal mantle heating relative to the heat ¬‚ux
coming from the core, although the topography associated with plumes must be
corrected for non-adiabatic effects. The Mesozoic and Cenozoic circulation of the
mantle can be inferred from plate motion histories, although a limitation is the lack
of information on initial conditions.

14.1 Convective forces
We start by providing a simple treatment of convection based on boundary layer
theory, and we will see that this provides velocities of the size associated with plate
motion (Figure 14.1).

14.1.1 Boundary layer theory
The Earth has much weaker heat diffusion than momentum transport, so that the
Prandtl number is large (∼1023 as discussed in Section 7.1.2). If we take the limit
of in¬nite Prandtl number the momentum equation represents a balance between
viscous forces and the forces associated with thermal and other buoyancies. This
insight puts us into a position to derive a a simple analytic model of convection

14.1 Convective forces 331


D d
T+ δT

Figure 14.1. Simple convection model with force balance associated with a downgoing
slab of lower temperature than its surroundings.

where the force balance is computed from estimates of thermal buoyancies and
the mantle viscosity. The approach was ¬rst used by Turcotte & Oxburgh (1967)
to show that it is possible to obtain a reasonable approximation to the velocity of
tectonic plates. Consider the weight of a subducting slab of length D:
FW = g Dd ρ±th ”T, (14.1.1)
where ±th is the thermal expansion coef¬cient (6.1.11), d is the thickness of the slab
and D the depth of penetration. ”T represents the average temperature difference
between the mantle and the slab, which is approximately ”T = ’T/2, where T is
the interior temperature. Hence
FW = ’g Dd ρ±th T/2. (14.1.2)
The thickness of the slab is
d = (kts)1/2, (14.1.3)
where k is the thermal conductivity and ts is the age of the slab material (see
Section 12.2). In terms of the distance from the ridge axis to subduction L and the
plate velocity v, the slab thickness
d = (kL/v)1/2. (14.1.4)
The viscous resistance to the penetration of the slab into the mantle is given by the
force per unit area,
σ = ·2v/D. (14.1.5)
We get the force per unit length by multiplying σ by the vertical length D:
FR = Dσ = 2·v. (14.1.6)
The force balance between the weight of the slab FW and the resistive force FR
yields the downgoing velocity as
v = ’gDdρ±th T/4·, (14.1.7)
332 Mantle Convection

but we must recall that the plate thickness d depends on the velocity v through
(14.1.4). The ¬nal expression for the velocity thus takes the form
gρ±th T k1/2 LD2 3
v= . (14.1.8)

In simple convection systems the aspect ratio is close to unity so it is reasonable
to take D ≈ L. We can obtain an estimate of the velocity v by using appropriate
parameters for the Earth™s mantle:
ρ = 4000 kg m’3 ,
D = 3000 km,
±th = 2—10’5 K’1 , T = 1400—¦ C,
k = 10’6 m2 s’1 , · = 1022 Pa s,
we ¬nd v = 2.8—10’9 m s’1 = 90 mm yr’1 . Even if we assume that slabs do not
penetrate further than half the mantle so that D ≈ L/2, we obtain v ≈ 30 mm yr’1 .
These two values bracket most of the inferred plate speeds.

14.1.2 Basic equations
We wish to focus attention on longer time scale features in ¬‚uid motion and
so need to suppress acoustic waves. We can do this by making the anelastic
approximation ‚ρ/‚t = 0, which yields the equation of mass conservation in the
∇ · (ρv) = 0. (14.1.9)
A further restriction to an incompressible medium leads to the vanishing of the
divergence of the velocity ¬eld (∇ · v = 0).
The equation of motion for ¬‚uid ¬‚ow in the absence of bulk viscosity takes the
form (7.1.4):
‚vj ‚vi
‚ ‚ ‚ ‚ ‚vk
’ 2 ·δij
ρ vj +ρ vk vj = ρgj ’ p+ · + .
‚t ‚xk ‚xj ‚xi ‚xi ‚xj ‚xk
If we allow a temperature dependent rheology the viscosity · will vary with
position through the variations in temperature and so we need to retain the spatial
derivatives of · in (14.1.10).
The ¬‚uid ¬‚ow equations are linked to the behaviour of temperature through the
transport of heat in convection. The temperature equation in the presence of ¬‚ow
is given by (7.1.15):
‚ ‚ ‚ ‚ ^^
ρ (CpT ) = ρvk (CpT ) + k T + 2·DijDij + h, (14.1.11)
‚t ‚xk ‚xk ‚xk
where k is the thermal conductivity, Cp is the speci¬c heat at constant pressure and
h is the heat generation per unit mass due to radioactivity; Dij is the deviatoric
14.1 Convective forces 333

strain-rate tensor. The dominant thermal behaviour is described by the thermal
diffusivity (7.1.17), κH = k/ρCp. For most of the mantle the conductive term will
dominate the viscous dissipation by a factor of 50 or more; however, signi¬cant
dissipation may arise in regions of strong shear such as the edges of thermal
boundary layers or rising plumes. In (14.1.11) we have retained the possibility
of spatial variations in speci¬c heat Cp and thermal conductivity k. We can expect
slow variations of these quantities with depth, but lateral variations induced by
temperature would normally be small. The additional effects come from the spatial
gradients of Cp and k and so will be concentrated in regions of sharp temperature

Reference state
To complete the picture we need to understand the relation of the density to
temperature and composition. We envisage a reference state in which the density
has an adiabatic radial pro¬le in hydrostatic equilibrium, so that

’∇p0(r) + ρ0(r)g(r) = 0. (14.1.12)

The variations in density associated with deviations from the equilibrium
temperature T0(r) pro¬les due to convective ¬‚ow can then be approximated by

ρ = ρ0(r) 1 ’ ±th [T ’ T0(r)] , (14.1.13)

where the thermal expansivity ±th (6.1.11) can be expressed as
1 ‚ρ
±th = . (14.1.14)
ρ ‚T P

14.1.3 Boundary conditions
The viscosity of the mantle is much larger than that of the core or the air and water
above. Thus the appropriate boundary conditions for the velocity are free-slip, i.e.,
zero normal velocity and zero shear stress as in (8.2.4), at the core“mantle boundary
(rc ), and the free surface boundary (re ) :

vr = 0, σrθ = σrφ = 0, r = rc , re .
at (14.1.15)

The corrections needed to allow for dynamic distortions of the boundaries are small
and generally ignored.
The temperature at the outer boundary will be dictated by the temperature at
the ocean bottom or the continental surface, whereas, at the core“mantle boundary
the temperature is practically isothermal and equilibrated to that of the core. If
core heat ¬‚ux has strong geographic variation, the isothermal surface and the
core“mantle boundary will not coincide.
334 Mantle Convection

14.1.4 Non-dimensional treatment
In our discussion of ¬‚uid ¬‚ows in Chapter 7, we have seen how considerable
insight into the physical character of the system can be obtained by working
with non-dimensional quantities, so that laboratory experiments and computer
simulations can be compared with conditions in the mantle through a set of
dimensionless numbers. The choice of non-dimensionalisation should be designed
to bring the main terms to comparable size.
A sensible choice for the length scale for mantle convection is the depth of the
mantle. LM = re ’ rc , i.e., 2890 km. There are two choices for the time scale
based on thermal diffusion or the advection time. We will here use the advection
time, Ta = LM/u, where u is the velocity of plates of the order 5 cm per year. The
corresponding ¬‚ow velocity UM = LM/Ta is around 10’6 m/s, and the effective
¬‚ow and diffusion time scales are then comparable. The set of main scaling factors
are thus:
L M = re ’ rc , Ta = , UM = . (14.1.16)
u Ta
Following scaling by the non-dimensional quantities we can write the equation
governing the behaviour of the convecting ¬‚uid in the form
‚ ‚
Re0 vj + vk vj =
‚t ‚xk
‚vj ‚vi
ρ0 ‚ ρ0 ‚ · 2· ‚vk
ρ0gj ’ p+ + ’ 3 · δij ‚x .
ρ ‚xj ρ ‚xi ·0 ‚xi ‚xj 0 k
where ρ0 and ·0 are the values for the reference distribution. The Reynolds number
Re0 is also constructed using these reference values
Re0 = = . (14.1.18)
·0 ν0
As we have seen in Section 7.2.1 the value of Re0 for the mantle of the Earth is very
small, around 10’20 and so the inertial forces on the left-hand side of (14.1.17) pale
into insigni¬cance relative to the viscous forces on the right-hand side. Viscosity
variations of even 1010 , equivalent to temperature changes of 1000 K, will not affect
this argument.
In consequence we have a Stokes ¬‚ow condition in the presence of gravitation
‚vj ‚vi
ρ0 ‚ ρ0 ‚ · 2· ‚vk
0 ≈ ρ0gj ’ p+ + ’ δij

ρ ‚xj ρ ‚xi ·0 ‚xi ‚xj ‚xk
in terms of the relative variations of density and viscosity.
The elliptic nature of the equation (14.1.19) has the consequence that the balance
of momentum is global and instantaneous. That is deformations in some place (i.e.,
14.1 Convective forces 335

plate motion) have the consequence of instantaneous stress variations throughout
the volume. This character of creeping ¬‚ow is important in a number of ways:
(i) for attempts to infer the ¬‚ow throughout the mantle from a knowledge of
plate motion (i.e. the insight of Hager & O™Connell, 1979),
(ii) for attempts to tackle the inverse problem of mantle convection where ¬‚ow
is inferred back in time (through the adjoint problem) knowing the history
of plate motion, and
(iii) for discussions of a decoupling layer in the asthenosphere; for reasonable
values of the low-viscosity channel, a factor of 100 or so relative to
the deeper mantle, the nature of the equation (14.1.19) prevents such
An alternative development for the non-dimensionalisation using the thermal
diffusion time, rather than the advective time, leads to a representation in terms of
the Rayleigh number Ra rather than the Reynolds number Re (Landau & Lifshitz,
1987, §57).

14.1.5 Computational convection
Having derived a mathematical representation of mantle ¬‚ow in the form of the
Navier“Stokes equations, we require numerical techniques for their solution. We
wish to focus on three-dimensional spherical convection and seek a discretisation
of the sphere. An obvious, but ultimately unsuitable, choice is a grid based directly
on longitude and latitude. Such grids are at a disadvantage, because there are many
more grid points near the poles than at the equator so that the polar regions are
oversampled, and because the convergent latitude and longitude lines result in two
singularities at the poles so that exceedingly small time steps are required to ensure
numerical stability. Collectively these dif¬culties are known as the pole problem.
We can avoid the pole problem if we discretise the sphere with a triangular mesh
derived from the icosahedron. This so called geodesic grid starts from projecting
the regular icosahedron, a platonic body with twenty equilateral triangles as its
faces, upon the sphere. The mesh is then constructed recursively by splitting nodal
distances in half and inserting new nodes at the midpoints. Each time the process
is repeated the lateral resolution in the mesh is doubled, and so a mesh of arbitrary
re¬nement on the sphere can be found. Recursive geodesic mesh re¬nements taken
from the three-dimensional spherical TERRA mantle dynamics code of Bunge &
Baumgardner (1995) are shown in Figure 14.2. The volume of the mantle can be
discretised by replicating the icosahedral mesh along the radial direction. The mesh
then subdivides a thick spherical shell into elements of almost uniform size.
Once the mantle volume is discretised, we need to ¬nd a local representation of
the variables and their derivatives. Many mantle convection studies adopt a ¬nite
element or ¬nite volume approach, because the local character of these techniques
is well suited for parallel computers and accommodates large local variations in
336 Mantle Convection

Figure 14.2. Successive re¬nements of icosahedral grid used to represent spherical be-

viscosity. In ¬nite elements one adopts a variational, or weak formulation of
the problem and represents the solution locally through low-order shape, or basis
functions. These functions are chosen in such a way that they vary between zero
and one within an elementary element of integration but vanish elsewhere, that is
they have a local support. A comprehensive survey of the ¬nite element method in
many circumstances in ¬‚uid mechanics is provided by Zienkiewicz et al. (2005).
We illustrate the basic ¬nite element procedure through the equation of motion,
and start by representing velocity and pressure locally through
v= Ψl vlj, j = 1, 2, 3 , (14.1.20)

where the Ψl are the ¬nite element basis functions. The weak form of the equation
of motion is derived by multiplying the momentum equation by appropriate trial
functions (φ) and then integrating the resulting system over the mantle domain ©,
‚vj ‚vi
‚ ‚ ‚vk
’ 2 ·δij
φ ρgj ’ p+ · + d© = 0, (14.1.21)
‚xj ‚xi ‚xi ‚xj ‚xk

so that we enforce the balance of momentum in an integral (weak) sense. A further
step is taken by using integration by parts on the viscous term in (14.1.21)
‚vj ‚vi
‚ ‚vk
’ 2 ·δij
φ · + d© =
‚xi ‚xi ‚xj ‚xk
‚vj ‚vi
‚φ ‚vk
’ 2 δij
’· + d©
V ‚xi ‚xi ‚xj ‚xk
‚vj ‚vi ‚vk
’ 2 δij
+ ·φ + dS, (14.1.22)
‚xi ‚xj ‚xk

which reduces the order of differentiation for this term. It is common to take trial
functions (φ) and ¬nite element basis functions (Ψ) as identical, in the Galerkin
approach to the ¬nite element technique. Upon inserting the local representation of
trial functions and velocity (14.1.20), the volume integral on the right-hand side of
14.1 Convective forces 337

(14.1.22) reads
‚Ψm ‚Ψl ‚Ψl ‚Ψl
vli ’ 2 δij
’ · vlj + v d© = Alimj vlj (14.1.23)
‚xk lk
‚xi ‚xi ‚xj

where A is the volume part of the viscous ¬nite element tensor operator for a
compressible ¬‚uid with variable viscosity
‚Ψl ‚Ψm 2 ‚Ψl ‚Ψm ‚Ψl ‚Ψm
Alimj = ’ · ’3 + δij d©,
‚xj ‚xi ‚xi ‚xj ‚xk ‚xk
for l, m = 1, . . . , n, n, i, j = 1, 2, 3. (14.1.24)
The expression (14.1.24) shows us that viscosity variations enter the local ¬nite
element operator through the nodal values of the viscosity ¬eld ·. In the special
case of an incompressible ¬‚uid with constant viscosity (14.1.24), reduces to the
scalar Laplacian operator
‚Ψl ‚Ψm
Slm = ’ · d©. (14.1.25)
‚xk ‚xk
We may make a number of choices for the ¬nite element basis functions (Ψ); these
include linear, quadratic or cubic functions. However, the integration by parts in
the weak form of the momentum balance has the effect of reducing the order of
differentiation for the viscous term from second to ¬rst order, so that it is suf¬cient
to consider piecewise linear basis functions in the elemental operator. The TERRA
code adopts piecewise linear basis functions for velocities and piecewise constant
basis functions for pressure. Its elementary element of integration d© on the sphere
is that of a spherical triangle.
We have seen that the momentum balance in the mantle is instantaneous and
global due to the Stokes ¬‚ow condition for ¬‚ow at extremely low Reynolds number.
We express this condition in the elliptic form of the momentum balance. The
discrete representation of the momentum equation is an algebraic sparse matrix
system, with entries clustered around the diagonal, where the sparse character
results from the local support of ¬nite elements. Numerical schemes exist to solve
such sparse matrix systems directly, for example through Gaussian elimination.
However, when ¬ne resolution is required the elliptic equation systems become
rather large and are best tackled through iterative schemes. The key to solve such
systems ef¬ciently lies in the multigrid procedure (see, e.g., Briggs, 1987). This
approach is optimal in the sense that the computational cost scales linearly in the
number of unknowns, that is it costs the same per grid point to solve small or large
matrices. This makes multigrid competitive with spectral methods based on the
Fast Fourier Transform (FFT), without incurring the penalties associated with such
global methods when handling large viscosity variations or in implementation on
parallel computers. Multigrid methods excel because they rely on a hierarchy of
nested computational grids, so that near- and far-¬eld components of the global
momentum balance are effectively solved at once. The nested structure of the
338 Mantle Convection

icosahedral grid is particularly well suited for multigrid. A number of mantle
dynamics models, such as the TERRA code, rely on multigrid methods for their
We also need to include the condition for mass conservation, which enters the
momentum balance as a constraint on the velocity ¬eld v. We take the algebraic
weak form of the momentum and the mass conservation equation in the anelastic
(14.1.9) limit:
GT ρv = 0, Av ’ Gp = f, (14.1.26)
where G is a gradient matrix and A is the full tensor operator built from the
operator elements A (14.1.22) that revert to the scalar Laplacian (14.1.25) for
incompressible ¬‚ow without viscosity variations. The coupled system (14.1.26)
is a saddle-point problem, and we seek to solve it for a velocity ¬eld v that satis¬es
the divergence-free constraint for incompressible, anelastic ¬‚ow. Pressure is the
Lagrange multiplier associated with this constraint. A standard approach to solve
(14.1.26) is through a conjugate gradient algorithm, often referred to as the Uzawa
method for incompressible ¬‚ow, where an inner loop solves the elliptic problem
by means of multigrid for a velocity ¬eld v, while an outer loop enforces the
divergence constraint on the velocity through a conjugate gradient search. This
method has been adopted in the TERRA code.
The parabolic nature of the energy equation (14.1.11) differs from the
elliptic form of the momentum balance, and ¬nite volume methods provide
a straightforward way for its solution. Remember that temperature variations
(‚T /‚t) arise from advection (’v · ∇T ) through a divergence term, and
from diffusion (∇2T ) through a scalar Laplacian (assuming isotropic thermal
conductivity) in addition to the heat source term. The ¬nite volume method, like
the ¬nite element method, is a technique for solving the differential equation as
a system of algebraic equations. The name ¬nite volume refers to small volumes
surrounding each node point on a mesh. Volume integrals that contain a divergence
term are converted to surface integrals, using the divergence theorem. The requisite
terms are then evaluated as ¬‚uxes at the surfaces of each ¬nite volume. The P´ clet
number of the mantle is large, so mantle heat transport is controlled primarily
by advection. Thus, outside of thermal boundary layers, the divergence term
dominates, and the ¬nite volume formulation is very effective. The ¬nite volume
method can also easily be adapted to unstructured meshes, such as the icosahedral
grid. Further, because the ¬‚ux entering a given volume is identical to that leaving
the adjacent volume, the ¬nite volume method is conservative. This approach has
been adopted in the TERRA code, with a simple ¬nite difference scheme used
to handle diffusion. The energy balance, and thus the temporal evolution of the
temperature ¬eld (‚T /‚t), is computed from an explicit second-order Runge“Kutta
marching method in time.
It is common to use parallel computers for mantle convection simulations.
The hardware in use includes traditional super computers, groups of individual
14.2 Convective planform 339

machines connected into so-called clusters and new multi-core machines.
Explicit message-passing in the MPI implementation is the standard programming
approach to accommodate these diverse platforms. Finite element, volume, and
difference techniques are all well suited for parallel computers because their local
support offers inherent parallelism. Typically one subdivides the computational
domain into smaller regions through an approach known as domain decomposition.
We illustrate this for the icosahedral grid, and show a subdivision of the grid for
4 and 16 processors in Figure 14.3. Most of the work within each sub-domain
is performed independently of the others, and boundary data are exchanged as
required among them.

Figure 14.3. Domain decomposition for the icosahedral grid. The spherical shell is shown
with one portion moved into the foreground. (a) The grid is partitioned into four sub-
domains and a greytone scheme shows mapping to individual processors, with different
tones used for processors 1, 2, 3, 4. (b) A recursion strategy allows for ¬‚exibility in the
number of processors employed. In a second decomposition step, each sub-domain itself
is split into a further set of four, to yield a total of sixteen subsidiary domains. In this way
one can accommodate parallelism with many thousands of processors.

14.2 Convective planform
We describe the geometry of convection cells through their planform. The term
originates from ¬‚uid dynamic tank experiments, where ¬‚ow patterns are imaged
by projecting collimated light (i.e., with parallel rays) through a tank ¬lled with a
transparent ¬‚uid such as corn syrup. Because the refractive index of light depends
on the local ¬‚uid temperature, light rays are bent away from hot regions and
converge toward colder regions. In the projected view upwellings show up as dark
and downwellings as bright areas.
Much of the original studies of planform for mantle convection were done in the
laboratory, and showed that the horizontal spacing of upwellings and downwellings
is comparable to the depth of the convecting ¬‚uid. This agrees well with the
theoretical predictions from linear stability derived in Chapter 7; see Whitehead
(1988) for a comprehensive review.
340 Mantle Convection

(a) (b) (c)

(d) (e) (f)

(g) (h) (i)

Figure 14.4. (a) A snapshot in time of the temperature ¬eld for a simulation of an incom-
pressible (Boussinesq) and only internally heated mantle ¬‚ow (zero Reynolds number).
The Rayleigh number RaH based on internal heating is 4—107 . The upper 200 km (roughly
the depth of the upper thermal boundary layer) is removed, to permit a view on the convec-
tion planform below the boundary layer. The isoviscous reference model shows point-like
downwellings from the upper boundary layer and relatively short wavelength convection
cells, quite unlike the Earth. (b,c,d) same as (a) but with the viscosity of the lower mantle
set thirty times larger than that for the upper mantle. (e) The superadiabatic temperatures
for the compressible, purely internally heated, isoviscous reference mantle convection cal-
culation with RaH = 1.1 —108 . Isolated downwelling plumes from the upper boundary layer
dominate the planform. (f) Same as (e) except for the addition of 38% bottom heating from
an isothermal core. Convection is in¬‚uenced by largely axisymmetric upwellings (plumes)
from the lower thermal boundary layer. (g) Same as (e) except for the addition of an endo-
thermic phase change at 670 km depth with γc = ’4 MPa K’1 . Downwellings pause in the
transition zone, before entering the lower mantle. (h) Strati¬ed viscosity case: same as (e)
except that the viscosity of the lower mantle has been increased by a factor of 30. The plan-
form is dominated by long, linear downwelling sheets. (i) combination of depth-dependent
viscosity and bottom heating.
The dominant features in (a),(e),(f),(g),(h),(i) are cold downwellings. Hot upwellings are
prominent in (b),(c). Case (g) has a very warm upper mantle.
14.2 Convective planform 341

However, in looking at the Earth it is clear that the characteristic horizontal scale
of the planform is much larger than the mantle depth, as evidenced by the spacing
between subduction zones. The difference is particularly striking for the Paci¬c
plate, where the plate width (12,000 km) exceeds the mantle depth (3000 km) by
a factor of four. In Section 13.1.5, we examined the linear stability of a ¬‚uid with
a low viscosity asthenosphere, and found that the dominant unstable convective
mode is shifted to longer wavelengths. This shift favours convection cells that are
much wider than they are deep.
Growing computer power means that convective planforms can now be studied
through numerical simulations. Such computer simulations are particularly useful
to understand convective systems where physical properties vary with depth, which
are dif¬cult to examine in the laboratory. The most pronounced variation in
mantle properties with depth is arguably the large viscosity jump between the
upper and lower mantle. But other parameters, such as thermal conductivity and
thermal expansivity are known to vary with depth. A wide range of geodynamic
studies on the in¬‚uence of these parameters has yielded the following conclusions:
an elongated mantle convection planform is associated with a reduction of the
Rayleigh number with depth. In other words, a low thermal expansivity or a
high thermal conductivity reduces the convective vigour of the deep mantle much
like an increase in lower mantle viscosity. We present a group of examples of
convective planforms for three-dimensional simulations of mantle convection in
Figure 14.4 to illustrate the differences between different modes of heating and
viscosity variations.
We take isoviscous, incompressible convection with pure internal heating as the
most simple representation of ¬‚ow in the mantle. This reference model, Figure
14.4(a), produces isolated, point like downwellings from the cold, upper thermal
boundary layer. The convective planform for this simple case is dominated by
short length scales, which agrees with laboratory results as we have noted above.
However, a dramatic change in the planform is produced by a single parameter
variation from the reference model, with an increase of the lower mantle viscosity
relative to that for the upper mantle. Sinking plumes in Figure 14.4(a) give way to
sheetlike downwellings much like the subduction dominated planform of the Earth,
Figure 14.4(b)“(d): the three panels show different aspects of the temperature
distribution for the same convection scenario.
Compressibility effects do not signi¬cantly alter the character of the ¬‚ow.
Compressible, isoviscous, purely internally heated convection, Figure 14.4(e),
produces point-like downwellings from the upper thermal boundary layer much like
in the incompressible reference model. The parameters that do affect the planform
include bottom heating from the core, and the buoyancy-associated mineral phase
transitions in the transition zone. Core heating adds a hot thermal boundary layer
at the bottom of the mantle. Because the surface area of the core“mantle boundary
is only a quarter of that at the outer surface of the Earth, the geometrical spacing of
342 Mantle Convection

the upwelling plumes increases the horizontal distance between upwellings and
downwellings, Figure 14.4(f). Otherwise the planform with additional bottom
heating remains similar to the case of convection with pure internal heating.
The introduction of phase transitions into the mantle produces minor planform
effects. In the model shown in Figure 14.4(g) we use a Clapeyron slope of γc = ’4
MPa K’1 for a phase transition at 670 km depth. This value is already probably too
large to be representative of the γ-spinel to perovskite phase change and results in
closely spaced, sheetlike downwellings that pause in the transition zone.
In Figures 14.4(h),(i) we show compressible convection with an increase in lower
mantle viscosity, with either pure internal heating or a mix of internal and bottom
heating. These cases are dominated by the increase in viscosity in the lower mantle.

14.3 Thermal structure and heat budget
The thermal structure of the mantle is characterised by steep thermal gradients
concentrated in narrow thermal boundary layers and a nearly adiabatic temperature
pro¬le in the intervening regions. The presence of internal heat sources combined
with the very long time it takes to cycle the volume of the mantle through
the oceanic plates (of the order of 1“2 billion years) results in departures from
adiabaticity, because the assumption of constant entropy employed to calculate
the adiabat no longer holds for internally heated convection. The amount of heat
entering the mantle from the core can be inferred from buoyancy studies over
mantle hot spots.

14.3.1 Thermal boundary layers and the geotherm
The thermal diffusivity of the mantle is small, and so heat transport by conduction
is restricted to thermal boundary layers where temperatures increase rapidly with
depth. In the intervening well mixed regions the temperature distribution is nearly
adiabatic. The overall temperature pro¬le through the Earth with adiabatic and
boundary layer regions is called the geotherm; a comprehensive discussion has
been provided by Jeanloz & Morris (1986). Figure 14.5 shows the basic convective
geotherm for a whole-mantle convection model of Earth-like convective vigour.
Because of their large vertical temperature gradients, the thermal boundary
layers control, via thermal conduction, the in¬‚ux (through the bottom) and out¬‚ux
(through the top) of heat. Such boundary layers also dominate the dynamics of the
mantle through their large temperature and buoyancy contrasts, as we saw from the
boundary layer analysis in Section 14.1.1. Gravitationally unstable material from
the boundary layers must eventually either sink (if in the cold top boundary layer)
or rise (if in the hot bottom boundary layer), thus feeding upwelling or downwelling
convective currents.
In the mantle we identify the lithosphere and tectonic plates as the cold upper
thermal boundary layer. A hot thermal boundary layer also exists at the bottom of
the mantle, where it separates the mantle from the core.
14.3 Thermal structure and heat budget 343


Temperature [K]


e Geothe



Figure 14.5. Basic convective geotherm for a whole mantle convection model.

Adiabatic gradient
For a mantle that is well mixed in an isentropic state, the geotherm away from
thermal boundary layers is that of an adiabatic temperature gradient, which is given
‚T ‚T ‚S
= , (14.3.1)
‚p ‚S ‚p
S p T

from Maxwell™s relation
‚S ‚V
= , (14.3.2)
‚p ‚T
T p

and so the isentropic temperature gradient with respect to pressure is
‚T ‚T ‚V
= . (14.3.3)
‚p ‚S ‚T
S p p

We recall the de¬nition of the thermal expansion coef¬cient ±th (6.1.11), and the
speci¬c heat at constant pressure Cp (6.1.14):
1 ‚V ‚S
±th = , Cp = T . (14.3.4)
V ‚T ‚T
p p

The adiabatic derivative of temperature T with respect to pressure p from (14.3.3)
can then be recast as
‚T T ±th
= . (14.3.5)
‚p ρ Cp
344 Mantle Convection

For the Earth we have the hydrostatic equation (6.1.2)
= ’gρ, (14.3.6)
in terms of the gravitational acceleration g. For this hydrostatic case the adiabatic
temperature increase with depth is given by
‚T ‚T T±th dp
= = . (14.3.7)
‚r ‚p ρCp dr

Using typical values for thermal expansivity, density, and heat capacity yields an
adiabatic gradient of 0.5 K km’1 .

Heat sources
There is a persistent escape of heat from the Earth™s interior to the oceans and
atmosphere, and the total heat loss through the Earth™s surface is 44 TeraWatts
(TW). About 31 TW of this heat is lost through the oceanic plates; another 13 TW
comes from the continents. Roughly half the continental heat ¬‚ux (about 7 TW) has
its source in heat production from within the continental crust. The rest is supplied
via conduction from the mantle beneath. Thus the total mantle heat loss is about 37
Heat is produced in the mantle by the decay of radioactive isotopes. The most
important contributors are uranium (238 U, 235 U)), thorium (232 Th) and potassium
(40 K). Such isotopes are believed to be disseminated throughout the mantle, but are
expected to be concentrated in the thermal boundary layers at the top and bottom
of the mantle.
There is a tendency for U, Th, K to be segregated via melting processes and
to stay with one of the fractions. As a result the continental lithosphere that
represents a residue from multiple phases of melting has high concentrations of heat
producing isotopes. Crustal concentrations reach parts per million (ppm) for U and
Th. Potassium is much more abundant (of the order of per cent), but the radioactive
isotope (40 K) represents only 10’4 of total potassium and so the contribution to
total heat production is less than for uranium and thorium. Inside the mantle the
radioactive elements are thought to be about two orders of magnitude less abundant
than in the crust. Typically U, Th and K occur in similar proportions relative to one
other in the crust and mantle, even though their absolute concentrations can vary
greatly. The average abundance ratio Th/U is about 3.5“4, while the ratio K/U is
usually 1“2—104 .
A rough upper bound on the concentration of heat-producing elements in the
mantle can be made, if we assume the mantle is in a thermal steady state. In this
case the heat loss from the mantle must be balanced entirely by heat production
from radioactive decay. With a mantle mass of 4—1024 kg, we arrive at an upper
bound for heat production of 11—10’12 W kg’1 . Correcting for the heat produced
in the continental crust, the upper bound on mantle heat production is 9—10’12
14.3 Thermal structure and heat budget 345

W kg’1 . This upper limit is substantially larger than the heat production rate
inferred from geochemical studies of upper mantle rocks. For example, Jochum
et al. (1983) estimate that the heat production rate of upper mantle rocks is only
about 0.6“1.0—10’12 W kg’1 . If such heat productivity applied throughout the
mantle it would at most result in a surface heat ¬‚ux of about 4 TW, that is only one
tenth of the total surface heat loss.
We can account for some of the de¬cit by assuming that the mantle is not in a
thermal steady state, so that we assume that the mantle undergoes secular cooling,
so that the interior temperature pro¬le is progressively lowered by 50 to 100 K per
billion years. Mantle secular cooling of 70 K Gyr’1 would result in a surface heat
loss of 10 TW. Still it is not clear whether there is enough radioactivity inside the
mantle to account for the observed heat loss.

Jeanloz & Morris (1987) appear to have been the ¬rst to notice that secular cooling
and the presence of internal heat sources must act to move the mantle geotherm
away from adiabaticity. An easy way to see this is to focus on the effect of internal
heating, and to follow a mantle volume element on its way from the core“mantle
boundary to the surface. As the volume element rises there is cooling due to
adiabatic decompression, but at the same time the temperature increases in response
to the internal heat released from radioactive decay. The net result is to make the
radial temperature change smaller than the adiabatic gradient, so that the geotherm
is subadiabatic. A subadiabatic model geotherm derived from three-dimensional
spherical convection modelling is shown in Figure 14.6.
Analytic and computational studies suggest departures from the adiabat by as
much as 300“500 K in the mantle. The heat equation (14.1.11) allows us to
derive a simple scaling argument. Outside of thermal boundary layers, where
thermal gradients are necessarily large, we ignore the diffusive term k∇2T .
This simpli¬cation means nothing else than that heat is transported primarily by
advection, or alternatively that the P´ clet number in the mantle is large, as we have
seen before. We also ignore the effects of shear heating. With these assumptions
(14.1.11) reduces to static heating by internal heat generation:
DT h
= , (14.3.8)
Dt Cp
where DT /Dt is the total time derivative of temperature T moving with the
volume element. For internally heated convection all material must eventually cycle
through the upper thermal boundary layer (the lithosphere in mantle convection) to
lose its heat. Note that this behaviour differs from purely bottom-heated ¬‚uids,
where it is suf¬cient for material from the bottom boundary layer to cycle through
the upper thermal boundary layer. We can derive an estimate of the circulation
time scale (”t) from the geometry of plate tectonics. Taking the total length of the
oceanic spreading system as 65 000 km and the average plate velocity as 5 mm yr’1 ,
346 Mantle Convection

T [K]

0 b

(a) (b)

Figure 14.6. Convective geotherms with the reference adiabat removed, for (a) pure in-
ternally heated convection, and (b) convection with 50% core heating. The effect of core
heating is to make the geotherm closer to adiabatic.

some 3 km2 of ocean ¬‚oor is created each year. The same amount must enter
the mantle through subduction. With the assumption that slabs and the material
entrained with them are 200 km thick we infer that some 600 km3 material with
mass 2—1015 kg enters the mantle (with total mass 4—1024 kg) each year. At this
rate the entire mantle circulates through the upper thermal boundary layer on a
time scale ”tmantle of order 2 billion years (Gyr), which together with an assumed
internal heating rate h of 5“10—10’12 W kg’1 and a mantle heat capacity Cp of
1000 J kg’1 K’1 implies a non-adiabatic component of order 250“500 K in the
mantle geotherm.

14.3.2 Plates
Tectonic plates are the cold upper thermal boundary layer of the mantle. These
plates are the surface expression of mantle convection and accomplish the primary
mantle heat loss, that is they cool the mantle. In Section 12.2.1 we have derived
the thermal structure for such plates from a half-space cooling model. The plates
differ from the rest of the mantle ¬‚uid, because they are cooler and thus stronger
than the mantle, although, as we have seen in Section 12.2.4, their strength is not
unlimited and with suf¬cient stress they yield in the form of brittle fracture. Thus
the de¬ning feature of plates is their combination of strength and mobility.

Heat transport
Plates carry the majority of the mantle heat loss. The heat transport can be
estimated from equation (12.2.5), which relates the heat ¬‚ux to the plate age
kTM 1 a
Q(t) = ’ √ √ = ’√ , (14.3.9)
πκ ts ts

where we have set a = kTM/ πκ. The total sea¬‚oor area A is 3.1—108 km2 and
new sea¬‚oor is formed at a rate of 3 km2 per year. At this rate the entire ocean
14.3 Thermal structure and heat budget 347

¬‚oor would be replaced in about 100 Myr and plates at the time of subduction (ts)
would be on average 100 million years old. Taking the heat ¬‚ux (qs) for 100 Myr
old sea¬‚oor as 50 mW m’2 allows us to evaluate the integral heat ¬‚ux over the
whole sea¬‚oor from
ts √
Q= dt Sq = 2Sa ts = 2Aq. (14.3.10)

We ¬nd that the average heat ¬‚ux through the oceans is about 100 mW m’2 and that
the total heat ¬‚ux is 31 TW in agreement with observations. Thus plates account
for 75% of the total heat loss of the Earth and nearly 90% of the heat lost by the

In¬‚uence on mantle ¬‚ow
Plates affect the mantle through their mechanical strength, since they remain nearly
rigid over geological time and thereby stabilise the upper thermal boundary layer.
The plate effect is particularly evident in the subduction process. Whereas a ¬‚uid
boundary layer becomes gravitationally unstable through local Rayleigh“Taylor
instabilities, the high viscosity of plates ensures that downwellings enter the mantle
only at a few major subduction zones. The stabilising effect of plates is taken to
the extreme in the case of Venus and Mars, where there is only one single plate
with little apparent movement and where the planet convects in a so-called rigid-
lid regime.
The Earth is in the mobile-lid regime, that is it has a strong yet mobile lithosphere
(plates). Thus the long-wavelength planform of the mantle re¬‚ects both the
in¬‚uence of depth-dependent viscosity and the great mechanical strength of plates.
The essence of this effect has been captured by saying that plates organise the ¬‚ow.

(a) (b) (c) (d)

Angular Order

Figure 14.7. Near surface temperature (100 km depth) and the spectrum of thermal hetero-
geneity as a function of depth and angular order for three-dimensional spherical computer
models: (a) isoviscous, (b) strati¬ed viscosity convection models without the stabilising
effects of plate motion and a stiff lithosphere; (c) isoviscous, (d) strati¬ed viscosity with a
relatively stiff lithosphere
348 Mantle Convection

In Figure 14.7 we show results from the models of Bunge & Richards (1996) to
illustrate the plate effect for global incompressible, internally heated mantle ¬‚ow.
The planforms and the thermal heterogeneity spectra are shown for isoviscous
convection and convection with a high viscosity lower mantle, with and without
the stabilising effect of the plates. Both isoviscous convection with a plate
effect and strati¬ed viscosity convection without the plate effect are unable by
themselves to reproduce the observed planform of the mantle. The dominance
of long-wavelength heterogeneity in mantle convection results from a combination
of the two features: a signi¬cant increase in mantle viscosity with depth, and the
existence of strong surface plates. The large-scale character of mantle structure thus
arises from independently established physical mechanisms in otherwise simple
mantle convection models.

14.3.3 Hot spots and plumes
Up to this point we have encountered the cold, upper thermal boundary layer of
mantle convection in the form of plates. However, there must also be a hot, lower
thermal boundary layer at the bottom of the mantle associated with the conductive
heat loss from the core. Thermal instabilities arising from this hot lower boundary
layer are known as plumes, and constitute the other important mode of mantle
Tuzo Wilson introduced the concept of plumes in order to explain the anomalous
volcanic activity of hot spot regions, which are areas of extensive intra plate
volcanism not associated with plate margin volcanism. Later Jason Morgan
envisioned vertical conduits of hot thermal plumes rising through the mantle and
accomplishing a major component of heat transport through the mantle. Although
between 20 and 100 different hot spots have been proposed at different times,
only about some 40 prominent hot spots are identi¬ed by most authors, including
Hawaii, Iceland, R´ union, Cape Verde, and the Azores. Large continental volcanic
centres such as Yellowstone or Afar are also attributed to hot spots.
Our theoretical understanding of plumes stems primarily from careful laboratory
experiments on the basic ¬‚uid dynamics; a recent review is provided by Campbell
(2007). The various experiments have identi¬ed two prominent features in the
morphology of plumes: a long thin conduit which connects the top of the plume
to its base, and a large mushroom like head that expands in size as the plume
rises upward through the mantle. The plume heads form because hot material
moves upward through the conduit faster than the plume itself can rise through
the surrounding mantle. The head grows as it entrains material from the mantle.
Plume heads in the Earth can reach substantial size, exceeding 500—500—500 km
to produce a total volume of melt in excess of 107 km3 .
When plume heads encounter the base of the lithosphere they spread sideways,
and undergo decompression melting to form enormous volumes of basaltic magma,
called ¬‚ood basalts. Radiometric dating puts the durationof eruptions from the
14.3 Thermal structure and heat budget 349

plume heads at less than a few hundred thousand years. At this rate an area the size
of California or Germany could be covered by 4 km of basalt in less than 1 million
years to form a continental ¬‚ood basalt (if it erupts through continental crust) or an
oceanic plateau (if the eruption takes place in the oceans).
The location of prominent ¬‚ood basalts is shown in Figure 14.8. Examples of
continental ¬‚ood basalts are the Deccan traps and the Rajmahal traps in India, as
well as the Siberian traps of Asia. Other examples include the Karoo basalts in
South Africa, and the Parana basalts in South America which were linked with the
Etendeka basalts in Africa, prior to the opening of the South Atlantic. Oceanic ¬‚ood
basalts include the Ontong Java plateau of the southwest Paci¬c and the Maniheken
plateau of the Indian ocean. A comprehensive review of ¬‚ood basalt eruptions and
their possible relation to mass extinctions is provided by Courtillot & Renne (2003).





KAROO Reunion

Marion Kerguelen


Figure 14.8. Map of hot spots with well de¬ned tracks and ¬‚ood basalts at their origins
(indicated by small capitals).

Hot spot swells
Most hot spots are associated with distinctive bathymetric swells, where the ocean
¬‚oor is anomalously shallow for its age. Figure 14.9 shows the swells associated
with Hawaii and Iceland, through perspective views of ocean-¬‚oor topography. An
example of topographic effects on the continents is the elevated region around
Afar, in eastern Africa. Hot spot swells must be maintained by buoyancy forces
associated with the plume. If a strong plume is located beneath a fast-moving plate,
such as the Paci¬c plate, large amounts of anomalous topography can be created
each year. In the case of Hawaii, an area 1000 km wide is elevated by 1 km at a
velocity of 10 cm per year (the speed with which the Paci¬c plate travels across the
hot spot). Rigorous studies of hot spot bathymetry have been presented by Davies
(1988) and Sleep (1990) to place constraints on the amount of core heat entering
the mantle.
350 Mantle Convection

Figure 14.9. Hot spot swells around Iceland and Hawaii modify the bathymetry of the

Following the style of analysis of Davies and Sleep, we make the simplifying
assumption that the origin of the buoyancy ¬‚ux is purely thermal, and that plumes
come from the core“mantle boundary. We envision the plume as a vertical conduit
with radius (r) and vertical ¬‚ow velocity (v), so that the buoyancy ¬‚ow rate (b) is
b = g ”ρ πr2 v, (14.3.11)
where ”ρ is the density difference between the plume (ρp) and the surrounding
mantle (ρm). The buoyancy ¬‚ux must balance the weight (W) of the anomalous
b = W = g (ρm ’ ρw)ws es uplate . (14.3.12)
where ρm ’ ρw is the density difference between mantle and overlying sea water,
ws and es are the width and elevation of the swell, and uplate is the plate velocity.
For Hawaii this yields a ¬‚ux of 7—104 N s’1 .
For a purely thermal origin we equate the heat ¬‚ux with the plume ¬‚ux, since
both depend on the excess temperature ”T = Tp ’ TM of the plume relative to
surrounding mantle, the resulting density difference
”ρ = (ρp ’ ρm) = ρm±th ”T, (14.3.13)
The plume heat ¬‚ux is also controlled by the excess temperature ”T ,
Q = πr2 v ρm Cp ”T, (14.3.14)
and is related to the buoyancy ¬‚ux b by
Q= b. (14.3.15)
14.3 Thermal structure and heat budget 351

The heat ¬‚ux is proportional to the buoyancy ¬‚ux weighted by heat capacity, gravity
and the thermal expansivity; yet the expression (14.3.15) does not involve the
excess temperature ”T . Thus, the same buoyancy ¬‚ux b could be maintained by
a small material ¬‚ux with large excess temperature or a large material ¬‚ux with a
small excess temperature.
For Hawaii, using typical values for the physical parameters (Cp = 1000
J kg’1 K’1 , ±th = 3—10’5 K’1 ) we obtain a heat ¬‚ow Q = 2—1011 W. This ¬‚ow
represents about 0.5% of the global heat ¬‚ux.
Similar studies of plume ¬‚ux have been presented by Davies and by Sleep for all
major hot spots. Although there are about 40 hot spots, all of them are weaker than
Hawaii and many are very much weaker. The best estimate of the total plume heat
¬‚ux derived this way is about 2.3—1012 W (2.3 TW), which is about 6% of the total
global heat ¬‚ow. In addition to the heat carried by plume tails, we must account for
the heat transported by plume heads. From the frequency of ¬‚ood basalt eruptions
this has been estimated at about 50% of the heat carried by plume tails. Thus the
total heat ¬‚ux transported by plumes would be around 3.5 TW, less than 10% of
the total global heat ¬‚ow rate, but still an important component of the Earth™s heat

Hot Spot reference frame
A common feature of many hot spots is a well de¬ned track made up of chains
of volcanoes and underwater sea mounts aligned in the direction of motion of the
overriding plate. The best example is the Hawaiian“Emperor chain, which extends
nearly 4000 km from the big island of Hawaii to the Aleutians near Alaska and
shows a clear progression with the age of submerged volcanoes increasing with
distance from Hawaii. The rate of development of the hot spot track is about 90
mm yr’1 for the past 40 Ma, which indicates the average velocity of the Paci¬c
plate relative to Hawaii during this period.
When we correct for the motion of the overriding plates, hot spots appear to be
relatively stationary with respect to each other and the mantle. Of course, there
must be some relative motion, due to mantle convection. For example, hot spots
in the Paci¬c seem to move relative to those in the African hemisphere. Detailed
studies for Hawaii have detected relative motion with respect to the Indian and
Atlantic hot spots of the order of a few mm per year. The clearest evidence for hot
spot motion has emerged for Hawaii, which seems to have originated signi¬cantly
further north from its present location, and subsequently drifted south. Tarduno et
al. (2003) provide a summary of the relevant palaeomagnetic data.
Nevertheless, the overall motion among hot spots appears to be smaller by about
an order of magnitude than the average motion associated with plate tectonics,
making hot spot ¬xity a useful ¬rst order approximation for time periods of 50-100
million years. This observation has been used to de¬ne a global reference frame
tied to the deep mantle through hot spots, the hot spot reference frame.
352 Mantle Convection

Plume excess temperature and global mantle heat budget
A longstanding observation about plumes is that they have low excess temperature
relative to normal mantle (e.g., Schilling, 1991). Most volcanic rocks associated
with hot spots have basaltic composition with major element chemistry similar to
that of mid-ocean ridge basalt. The petrology of hot spot lavas can be explained
with an excess temperature of the order of 200“300 K. This low value for the excess
temperature would appear at odds with independent considerations that suggest
a much larger temperature change across the core“mantle boundary of the order
of 1000 K. If plumes originate from a thermal boundary layer at the core“mantle
boundary the difference in temperature contrasts is dif¬cult to understand.
One way to reduce the excess temperature is through a chemically dense layer
at the core“mantle boundary. Such a layer could buffer the temperature contrast
between the core and the mantle. Yet, two other considerations are probably more
important. Firstly, one expects a lowering of the excess temperature associated
with plumes because the mantle has internal heat sources. As we have seen before,

(a) Adiabat (b) Excess Temperature

2250 K
Temperature [K]

1750 K


1250 K

750 K 200
min mean max
0 20 40 60 80 100 120 140 800 1300 1800 2300 2800
Pressure [GPa] Footing Temperature [K]

Figure 14.10. (a) Curves of constant entropy adiabats for footing temperatures of 750,
1250, 1750, 2250 K (the adiabatic projection of the mantle temperature up to the Earth™s
surface). The slope of the adiabat increases with footing temperature, and generally de-
creases with depth. Temperature jumps corresponding to the discontinuities in entropy of
the phase transitions are less evident in the colder pro¬les. (b) Total adiabatic temperature
variation from surface to core“mantle boundary plotted against footing temperature. There
is an increase in adiabatic temperature variation with footing temperature. Also shown are
the footing temperatures corresponding to mean, minimum and maximum radial tempera-
tures in a spherical convection model. Hot thermal upwellings (plumes) undergo a stronger
adiabatic cooling relative to surrounding mantle, so that their excess temperature decreases
systematically in the mantle from the bottom to the top.
14.4 Circulation of the mantle 353

internal heating in combination with the slow overturn of the mantle results in the
mantle geotherm departing from the adiabat by as much as 250“500 K. Plumes,
however, rise through the mantle relatively quickly, on a time scale of the order
of 100 million years. The plume temperature distribution is therefore expected
to be nearly adiabatic. We can express this result in a different way: the normal
geotherm decreases less rapidly than the adiabatic gradient in moving from the
deeper to the shallower mantle, whilst the temperature inside the plumes follows
the adiabat closely. The net result is a systematic loss of excess temperature in the
plume as it moves away from the core“mantle boundary.
A second consideration is evident from equation (14.3.7), which states that the
adiabatic temperature change depends on temperature. For example, from Figure
14.10 we see that an isentrope tied to a footing temperature of 2000 K, at zero
pressure, undergoes a temperature increase with depth nearly twice that of an
adiabat footed at 1000 K. For plumes the adiabatic temperature drop therefore
is about 200 K larger than that of normal mantle, due entirely to their higher
Our observations bear on the heat loss of the core, which as we have seen can
be constrained from buoyancy studies over hot spots. If we correct for the steeper
adiabatic gradient in plumes (about 250 K) and the effects of subadiabaticity in
the mantle geotherm (about 250 K) the excess temperature in a plume and the
associated buoyancy ¬‚ux at the core“mantle boundary could be larger by a factor
of three than the valued inferred from surface observations. The total heat ¬‚ux
transported by plumes could thus be in the range of 10“12 TW, which represents
about 30% of the total mantle heat budget. Together with a secular cooling
contribution of 70 K Gyr’1 (10 TW) and a radioactive heating rate of 4“6—10’12
W kg’1 amounting to 16 TW, the mantle heat budget would be balanced.

14.4 Circulation of the mantle
It is common to use the term circulation to describe the motion of the mantle, in
analogy to the general circulation in the oceans and the atmosphere. The direct
evidence we possess for mantle circulation comes from plate tectonics, because
the viscosity of the mantle is suf¬ciently large to ensure coupling between plate
motion and mantle ¬‚ow. Geodynamic calculations suggest upper mantle viscosities
as low as 1017 Pa s are required before plate motion and deep mantle convection
are effectively decoupled. Thus the kinematics of plate motion models provides
¬rst-order constraints on the dynamics of the mantle. A close linkage between
plate motion and mantle circulation must also be inferred from the requirement to
conserve mass, because mass transfer associated with plate motions at the surface
must be balanced by internal mass displacements in the mantle. Plate motions
change over time, and these variations reveal important temporal variations in the
pattern of the internal circulation.
The most important data used to construct plate motion models come from the
354 Mantle Convection

magnetic isochron record of oceanic plates. Thus reliable reconstructions of plate
motion are restricted to the past 100“150 Myr, the age of the oldest ocean ¬‚oor.
Early in the chapter we saw that the primary mantle heat loss occurs through the
oceanic lithosphere. This means that a close linkage between plate motion and
mantle circulation must also be inferred from thermal considerations, and that the
large scale thermal structure of the mantle must be dominated by past subduction.
The combination of convection models with reconstructions of past plate motion
allows us to explore the evolution of the Mesozoic and Cenozoic mantle in
mantle general circulation models. Such models solve the mantle convection
equations subject to boundary conditions on surface velocities taken from the
history of plate motion. Mantle circulation models explain some deep mantle
heterogeneity structures imaged by seismic tomography, especially those related to
past subduction. However, their lack of initial condition information is a signi¬cant
limitation. The problem of unknown initial conditions can be addressed using an
inverse approach to convection, based on the exploitation of the adjoint equations
for mantle convection.

14.4.1 Present-day and past plate motion models
Models of present-day plate motion are given in terms of rotation poles and angular
velocities. We start from Euler™s displacement theorem which states that rigid body
motion at the surface of the Earth can be expressed as a rotation about an axis
passing through the Earth™s centre. This was evident to Alfred Wegener when he
described the drift of North America relative to Eurasia as the rotation around an
axis passing through Alaska. Euler™s theorem gives us a compact and quantitative
way to describe the Earth™s surface velocity ¬eld. In particular, plate velocities can
be described by an angular velocity vector and then standard vector algebra can be
applied to combine motions. Thus if both the angular velocity of plate A relative to
B, A ωB , and that of B relative to C, C ωB , are known, we obtain the motion of C
relative to A from vector addition:

C ωA = C ωB + B ωA . (14.4.1)
Angular velocities can be deduced from a variety of observations that include the
orientation of active transform faults between two plates. Since relative motion
between two plates sharing a mid-ocean ridge is parallel to transform faults, the
fault arcs must lie on small circles. The rotation pole then lies somewhere on the
great circle perpendicular to the transform faults. With information from two or
more transform faults the position of the rotation pole can be determined (see,
e.g., Morgan, 1968). Another way to examine the relative motion is to map
spreading rates along a mid-ocean ridge, because the magnetic anomaly pattern
varies as the sine of the angular distance from the rotation pole. A third method to
compute the relative motion between two plates is the use of fault plane solutions
(focal mechanisms) from earthquakes along plate boundaries, although this method
14.4 Circulation of the mantle 355

provides just the location of the pole, and not the spreading rate. A global model
for current plate motions based on these types of data has been constructed by De
Mets et al. (1990, 1994).
Today plate motions are measured with great accuracy using space geodetic
methods (cf. Section 12.4.2). Very long baseline interferometry (VLBI) exploits
very distant quasars as source signals and terrestrial radio telescopes as receivers.
The difference in distance between two telescopes is then tracked over a period
of years. The use of Global Navigation Satellite Systems (GNSS) such as the
Global Positioning System (GPS) is another common method to measure plate
motions. A good review on the application of such space geodetic methods in
a geological context is provided by Dixon (1991). GNSS measurements allow
real-time monitoring of plate motions; a survey time of about 1“2 years is normally
suf¬cient to obtain a meaningful estimate of the rate and direction of motion at a
site. A plate motion model based only on GPS data was constructed by Sella et al.
Models of past plate motion are constructed primarily by matching magnetic
anomalies and fracture zones of the same age, corresponding to patterns of
palaeo-ridge and palaeo-transform segments at a given reconstruction time. The
plates then need to be restored to their palaeo-positions through a series of ¬nite
rotations, often referred to as stage poles. Bullard, Everett & Smith (1965) used
¬nite rotations based on similarities in the Atlantic coastlines to publish the ¬rst set
of computer-generated plate reconstructions describing the opening of the Atlantic.
Figure 14.11 presents a modern summary of reconstructions of plate motion over
the past 130 Ma, a period that has seen signi¬cant shifts in the relative positions of
the continents.
Palaeomagnetic plate reconstructions constrain the palaeo-latitude and angular
orientation of a plate, but leave its palaeo-longitude unknown. To minimise the
uncertainty in longitude one is required to specify a reference frame. One choice
is the selection of a slow-moving plate, such as Africa, as the common reference
against which the motion of all other plates is restored. Africa™s lack of elongated
hot spot tracks is commonly taken as evidence that the plate has moved by less
than 15 degrees over the past 100 Myr. An alternative choice is the use of the no-
net-rotation reference frame, where one assumes that the mean lithosphere velocity
integrated over all plates must vanish. The inherent dif¬culty with this reference
frame rests in the fact that it is not obvious whether the mean velocity of plates
should vanish. Thus a third and physically more plausible reference frame is
often used by relating plate reconstructions directly to the deep mantle through
hot spots in the hot spot reference frame, as we have seen before. Steinberger &
O™Connell (1998) have tried to improve on the accuracy of the hot spot reference
frame by accounting explicitly for errors associated with relative hot spot motion
using global mantle ¬‚ow models. In general, plate motion models based on a hot
spot reference frame result in minimal longitudinal motion of Africa (compared
356 Mantle Convection

Figure 14.11. Reconstruction of plate motions for the past 130 Ma. [Courtesy of D.

with most other plates), and thus con¬rm the relative ¬xity of the African plate.
Finite rotation poles have been published for most plates describing their late
Cretaceous/Tertiary history of motion. A widely used plate motion model in the
hot spot reference frame for the Cenozoic is that of Gordon and Jurdy (1986).
14.4 Circulation of the mantle 357

100 90 80 70 60 50 40 30 20 10 0
Time [Ma]

Figure 14.12. Subduction history for the last 110 Ma, where lighter tones are used for
progressively younger subduction. Note the signi¬cant differences in current and past sub-
duction in the northwest and southwest Paci¬c. [Courtesy of B. Steinberger.]

A consequence of the changes in the motion of the plates is that the con¬guration
of subduction zones has changed, as can be seen from Figure 14.12 where the
position of subduction is coded by age, with lightest tones used for the most recent
period. For the last 100 Ma, the perimeter of the Paci¬c has been the location of
most of the subduction zones. The current sites of subduction show up clearly in
seismic tomography (Section 11.4.4) for the upper mantle. Clear indications of the
former Farallon plate beneath the Americas and the subduction beneath South Asia
at the northern edge of the Tethys ocean are found in the upper part of the lower
mantle (see, e.g., Figure 11.23 at 1100 km).

14.4.2 Implications of plate motion models for mantle circulation
We can use plate motion models to make a number of direct predictions about
structure and dynamics of the mantle. For example, the temperature variations in
the mantle should be dominated by heterogeneities associated with past subduction
re¬‚ecting the fact that vast amounts of old oceanic lithosphere descended into the
mantle over the past 100 Myr (see Figure 14.12). In a classic paper, Chase &
Sprowl (1983) point to a rough correspondence between major geoid lows and
Cretaceous subduction zones, while Anderson (1982) noted the complementary
correspondence of the African geoid high with the ancient supercontinent Pangea.
Section 13.4.1 showed us that the geoid is dominated by these broad highs and
lows, and as a result the geoid spectrum peaks at spherical harmonic degrees l
of 2 or 3. The time lag between Mesozoic plate tectonics and modern geoid
358 Mantle Convection

180° 180°

(a) TPW every 5 Ma (b) TPW every 1 0 Ma
(1 0 Ma window) (2 0 Ma window)

19 39 29
8 18
60 3 59
64 50 52
26 54
68 75
270° 90°
142 122 119 94 90
104 83 97
126 113
130 100
139 110
168 134 151
170 155


. 13
( 16)