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

u

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

u

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

X T T

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.

14

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

330

14.1 Convective forces 331

L

T

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

2

gρ±th T k1/2 LD2 3

v= . (14.1.8)

4·

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

form

∇ · (ρ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+ · + .

3

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

(14.1.10)

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

changes.

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:

LM LM

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

(14.1.17)

where ρ0 and ·0 are the values for the reference distribution. The Reynolds number

Re0 is also constructed using these reference values

ρ0UMLM UMLM

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

3·

ρ ‚xj ρ ‚xi ·0 ‚xi ‚xj ‚xk

0

(14.1.19)

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

decoupling.

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-

haviour

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

n

v= Ψl vlj, j = 1, 2, 3 , (14.1.20)

l=1

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)

3

‚xj ‚xi ‚xi ‚xj ‚xk

V

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

3

‚xi ‚xi ‚xj ‚xk

V

‚vj ‚vi

‚φ ‚vk

’ 2 δij

’· + d©

3

V ‚xi ‚xi ‚xj ‚xk

‚vj ‚vi ‚vk

’ 2 δij

+ ·φ + dS, (14.1.22)

3

‚xi ‚xj ‚xk

S

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

3

‚xi ‚xi ‚xj

V

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

V

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

V

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

performance.

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

e

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

3000

2500

Temperature [K]

2000

rm

e Geothe

onvectiv

1500

C

1000

500

0

Top CMB

Depth

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

by:

‚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

S

344 Mantle Convection

For the Earth we have the hydrostatic equation (6.1.2)

dp

= ’gρ, (14.3.6)

dr

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

dp

= = . (14.3.7)

‚r ‚p ρCp dr

dr

S S

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

TW.

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.

Subadiabaticity

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

e

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

800

400

T [K]

0 b

a

-400

(a) (b)

-800

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)

0

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

mantle.

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)

Depth

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

convection.

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

e

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

NATB

Iceland

SIBERIAN

Yellowstone

CRB

DECCAN

RAJMAHAL

Hawaii

CARIBBEAN ETHIOPIAN

Afar

ONTONG JAVA

Galapagos

PARANA

KAROO Reunion

Louisville

TASMANIAN

Tristan

Marion Kerguelen

Balleny

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

oceans.

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

topography

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

Cp

Q= b. (14.3.15)

g±th

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

budget.

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

3500

2250 K

1400

3000

TCMB - TTOP [K]

Temperature [K]

2500

1000

1750 K

2000

600

1250 K

1500

1000

750 K 200

min mean max

500

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

temperature.

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.

(2002).

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.

Mueller.]

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)

39

33

30

19 39 29

14

11

8 18

60 3 59

46

64 50 52

26 54

67

68 75

270° 90°

77

79

91

142 122 119 94 90

104 83 97

136

116

126 113

118

130 100

139 110

158

168 134 151

162

153

170 155

199

196

173

176

194