. 12
( 16)


to ¬nd models of seismic anisotropy that are compatible with both the sets of
observations suggest that the evanescent tail of the surface wave displacements
with depth leads to an underestimate of the in¬‚uence of anisotropy within the
asthenosphere. The presence of a zone of asthenospheric ¬‚ow may therefore be
a common feature between the continents.
The presence of such ¬‚ow de¬‚ected by the structure at the base of the lithosphere
may help to explain the complex pattern of hotspots tracks and volcanic edi¬ces
in southeastern Australia (W. Jason Morgan, personal communication). The
asthenosphere ¬‚ow is guided away from the deep rooted cratons in central Australia
towards the east by the set of lithospheric steps (Figure 13.2). As the material rises
there is consequent decompression and potential release of melt, which can reach
the surface in zones of pre-existing weakness.
The patterns of shear wave splitting in eastern North America with a distinct
change in the orientation of the fast direction of wave propagation in the New
England region over a relatively short distance require a change in the controlling
factor for the seismic anisotropy. Fouch et al. (2000) have suggested that this
arises from the de¬‚ection of asthenospheric ¬‚ow around the deep lithospheric keel
beneath the cratonic region that appears to have a re-entrant feature near 200 km
depth from surface wave tomography studies. Finite difference models of the ¬‚ow
pattern to be expected around a mantle keel indicate that a good ¬t to the observed
patterns can be achieved with a segment (˜divot™) missing from the base of the
keel. The mantle ¬‚ows around the deep lithosphere but tends to enter into the gap,
leading to rapid variations in the expected fast directions of shear waves over short
distance scales. Beneath the craton root the fast directions are parallel to the plate
motion, which is consistent with general patterns of observations.

13.1.5 The in¬‚uence of a low-viscosity zone
The presence of a zone of low viscosity modi¬es the way in which ¬‚uid ¬‚ow and
convective behaviour is organised. Busse et al. (2006) have provided an analytical
13.1 Lithosphere and asthenosphere 303

approach that allows the investigation of a region of lowered viscosity adjacent
to material with an extended zone of uniform viscosity. The con¬guration of the
system is sketched in Figure 13.6; it is convenient to take a symmetric scenario
with both upper and lower viscosity layers. This provides a simple approximation
to the presence of an asthenosphere and a lower mantle boundary layer.
We consider convection in a ¬‚uid layer of height h with constant kinematic
viscosity ν0 in most of the interior, but with a much lower value νb in thin layers
of thickness δh ( h) next to the top and bottom boundaries (Figure 13.6). For an
incompressible ¬‚uid (∇ · v = 0), the Navier-Stokes equation (7.1.5) is
‚vj ‚vi
‚ ‚ ‚ ‚
ρ + vk vj = ρgj ’ p+ · + , (13.1.1)
‚t ‚xk ‚xj ‚xi ‚xi ‚xj
with the linked thermal equation for material with constant thermal diffusivity κH
(7.5.7) ,
‚ ‚ ‚‚
+ vk T = κH T, (13.1.2)
‚t ‚xk ‚xk ‚xk
where we have neglected the work done by the weak ¬‚ow and assumed no internal
heat sources.
In the Boussinesq approximation, as in Section 7.5.2, we include thermally
induced density variations in the buoyancy term (ρgj) but neglect such variations
elsewhere. We can absorb the hydrostatic pressure in the reference state as in (7.5.6)
by introducing
P = p ’ ρ0gz. (13.1.3)

νb δh




Figure 13.6. Con¬guration of a convection system with zones of lowered viscosity, ¬‚ow
directions are indicated with open arrows.
304 The In¬‚uence of Rheology: Asthenosphere to the Deep Mantle

The buoyancy contribution is then ’ρ0±th T ’ T0 g(^z). e
It is convenient to work with non-dimensional quantities and we use the height
of the layer h as the length scale, the thermal timescale h2/κH , and the difference
between the top and bottom temperatures (T2 ’ T1) as the temperature scale. In
terms of non-dimensional quantities the Navier-Stokes equation takes the form
‚vj ‚vi
1 ‚ ‚ ‚P ‚ν
+ vk vj = ’ + Ra ˜(^z)j + + ,(13.1.4)
Pr ‚t ‚xk ‚xj ‚xi ν0 ‚xi ‚xj
where we have introduced the Prandtl number Pr (7.1.18) and the Rayleigh number
Ra (7.5.1)
±th g(T2 ’ T1)h3
Pr = , Ra = . (13.1.5)
κH ν0κH
The equation for the non-dimensional temperature ˜ has the simple form
‚ ‚ ‚‚
+ vk ˜= ˜. (13.1.6)
‚t ‚xk ‚xk ‚xk
We examine two-dimensional solutions of (13.1.5), (13.1.6) in the circumstances
when the Prandtl number is very large. As in Sections 7.4, 7.5.2 we introduce a
stream function ψ such that
‚ψ ‚ψ
vx = ’ , vz = . (13.1.7)
‚z ‚x
We can eliminate the term from (13.1.5) by taking the curl of the equation and
then from the x-component we ¬nd the following equation system for the stream
functions ψ0 in the interior and ψb in the boundary layers
∇4ψ0 ’ Ra = 0 in ’ 1 + δ < z < 1 ’ δ, (13.1.8)
2 2
∇4ψb ’ βRa = 0 in ’ 1 ¤ z ¤ ’ 1 + δ, 1 ’ δ ¤ z ¤ 1 . (13.1.9)
2 2 2 2
Here β is the ratio of the interior viscosity to that in the boundary layers ν0/νb.
With the assumption of stress-free boundaries at z = ± 1 we require

ψb = ‚2 ψb = 0 at z = ± 1, (13.1.10)
zz 2

where we have used a contracted notation for the partial derivative. We restrict
attention to one convection cell so
ψb = ψ0 = 0 at x = ±l. (13.1.11)
We require continuity of the tangential stress at the edge of the zone of lowered
viscosity, so that at z = ’ 1 + δ, 1 ’ δ
2 2

[‚2 ’ ‚2 ]ψb = β ‚2 ’ ‚2 ψ0.
ψb = ψ0, ‚zψb = ‚zψ0, (13.1.12)
zz xx zz xx

In the limit of very large Rayleigh numbers we will get a periodic array of
13.1 Lithosphere and asthenosphere 305

alternate rising hot plumes and descending cold plumes with a very narrow width.
We restrict attention to the interval ’l ¤ x ¤ l with a hot plume at x = ’l and a
cold plume at x = l. Away from these regions the horizontal derivative ‚x˜ will be
very small and can be neglected in (13.1.8), (13.1.9). Without loss of generality we
can set the mean temperature to zero and so the temperature boundary condition is
˜ = “1 z = ± 1.
at (13.1.13)
2 2

Stream function
We assume symmetry about the mid-plane of the layer z = 0 and concentrate
attention on the upper zone. We need to include the condition of continuity of
normal stress at the viscosity interface z = 1 ’ δ. To a good approximation the
effect of the normal viscous stress in the interior layer acts on the less viscous layer
as a compensating pressure term at z = 1 ’ δ,

‚xpb(x) = β 3‚3 + ‚3 ψ0 (13.1.14)
xxz zzz

In the thin low-viscosity layer this pressure drives a shear ¬‚ow speci¬ed by

∇2‚zψb = ‚xpb(x) = β 3‚3 + ‚3 ψ0 = βΨb. (13.1.15)
xxz zzz
z= 1 ’δ

Since the thickness δ of the layer is small we can replace the Laplacian by ‚zz and
integrate with respect to z to obtain
ψb = 1 βΨbz3 + f(x)¯, with z = z ’ 1 ,
¯ ¯ (13.1.16)
6 2

where f(x) is a function to be determined; this form for ψb satis¬es the continuity
conditions (13.1.10). The continuity conditions for tangential stress, (13.1.12),
require that at z = 1 ’ δ the interior stream function ψ0 satis¬es

ψ0 = ’δ3 6 βΨb ’ δc(x),

‚zψ0 = δ2 2 βΨb + c(x),

‚2 ’ ‚2 ψ0 = δβ 3‚3 + ‚3 ψ0. (13.1.19)
zz xx xxz zzz

We can extract f(x) from (13.1.17), (13.1.18) as
f(x) = ’ 1 1
ψ0 + ‚zψ0 , with d = ’ δ. (13.1.20)
2 2
δ z=d

The solution of the biharmonic equation ∇4ψ0 = 0 that meets the horizontal
boundary conditions (13.1.11) is

2 2
ψ0 = 2 C(x ’ l ) + cos(amx){Am cosh(amz)+Bm sinh(amz)} ,(13.1.21)
am = (2m ’ 1) . (13.1.22)
306 The In¬‚uence of Rheology: Asthenosphere to the Deep Mantle

The solution for Am, Bm exploiting the continuity conditions (13.1.17) “ (13.1.19)
is straightforward but algebraically complex (for details see Busse et al., 2006). To
a good approximation the solution for ψb can be written as

z3 z
¯ ¯
ψb = 2 “m cos amx ’3 2 + Gm(δ) + ... (13.1.23)
δ3 δ

where “m = (’1)mC/(a3 l) and

δ3βa4 d ’ 3 cosh2(amd)
Gm(δ) = (13.1.24)
3βa3 [cosh(a d) sinh(a d) ’ a d] ’
δ 2 cosh (amd)
m m m

The expression (13.1.23) can be well approximated by
z3 z
¯ ¯
2 2
ψb = 4 C(x ’ l ) ’3 2 + G1(δ) , (13.1.25)
δ3 δ
since terms with m ≥ 2 make only a minor contribution.

Thermal boundary layer
The stream function (13.1.25) with a parabolic pro¬le describes the dominantly
horizontal ¬‚ow in the upper layer with low viscosity, and can now be used to
consider the thin thermal boundary layer close to z = 1 . The temperature equation
(13.1.6) can be written in terms of the stream function as

∇2˜ = ‚zψb‚x˜ ’ ‚xψb‚z˜. (13.1.26)

With the introduction of the coordinates
ζ = 1 ‚zψb z,
ξ= d˜ ‚zψb(˜ ),
x x ¯ (13.1.27)

the temperature equation can be transformed into a diffusion equation
‚˜ 1‚ ˜
= 4 ‚ζ2 . (13.1.28)
A solution which satis¬es the transformed boundary conditions

˜ = ’1 ζ = 0, ˜ = 0 for ζ ’ ∞,
at (13.1.29)

can be found from the results in Section 8.4.2 as

ζ/ ξ
1 2
dy e’y .
˜ = ’1 ’ √ (13.1.30)
2 π 0

The temperature distribution (13.1.30) in the boundary layer does not link directly
to the rising plume at x = ’l, but numerical solutions indicate that this does not
represent a problem for positions close to the upper boundary.
13.1 Lithosphere and asthenosphere 307

The local Nusselt number Nuloc can be found from the normal derivative of the
temperature at z = 1 since the heat transport in the absence of convection is unity
‚ C(2 + G(a1) l’x
Nuloc = ’ ˜(ξ(x)) = . (13.1.31)
‚z πδ (2l ’ x)1/2
The regular Nusselt number is the average of (13.1.31) over the interval ’l ¤ x ¤ l
C(2 + G(a1) 2
Nu = . (13.1.32)
2 πδ
The heat transport arriving at the upper boundary has to equal that carried by the
rising plume
’l+ ’l+
2lNu = dx ˜ [’‚zψ]x=’l = lC[1 + H(a1)] dx ˜, (13.1.33)
’l’ ’l’

where the angle brackets denote the average over the z-dependence. The function

4 δ32βa2(sinh2(a1d) ’ 2(6 sinh(a1d) cosh(a1d) ’ a1)
H(a1) = 2 ;(13.1.34)
32βa3(sinh(a d) cosh(a d) ’ a d) + 3 cosh2(a d)
πδ 1 1 1 1

terms with m ¤ 2 make a negligible contribution. The constant C can be evaluated
by balancing the tangential stress on the plane x = ’l with the buoyancy of the
rising plume,
[2 + G(a1)]l
C= 2 Ra dx ˜ = Ra . (13.1.35)
4π[1 + H(a1)]2δ

Hence the Nusselt number can be expressed as
[2 + G(a1)]2l2
Nu = 2 Ra . (13.1.36)
2π[1 + H(a1)]δ2
For a constant viscosity layer the maximum heat transport occurs for an aspect
ratio 2l of order unity, whereas with the presence of the low viscosity layer the
maximum of the Nusselt number occurs for somewhat higher aspect ratios. For a
viscosity ratio of 1000 the horizontal cell size doubles. The approximate treatment
gives the correct dependence of Nu on l but somewhat underestimates the actual
This simple model indicates that a region of reduced viscosity in the upper
mantle as indicated from studies of glacial rebound (Section 12.5) has important
consequences, particularly with respect to the aspect ratio of possible convection
cells. The aspect ratios needed for a convective origin for plate motion are
far more consistent with the presence of a zone of lowered viscosity than the
constant-viscosity case.
308 The In¬‚uence of Rheology: Asthenosphere to the Deep Mantle


Island 4
Overriding Plate Subducting




Figure 13.7. Cross-section through a subduction zone system indicating the nature of ¬‚ow,
the stress state of the slab through observed earthquake mechanisms and the typical loca-
tion of the volcanic arc associated with subduction. (1) volcanic arc, small earthquakes
occur in the over-riding plate in the forearc; (2) seismogenic zone associated with thrust
earthquakes, e.g., the great 2004 Sumatra“Andaman event (Mw 9.3); (3) normal faulting
through the ocean lithosphere as in the 1997 Sumbawa event (Mw 7.7); (4) fracturing of
plate as it approaches trench and starts to bend; (5) dehydration reactions in former oceanic
crust, intermediate earthquakes near the top of the slab with largely down-dip tension; (6)
volatiles released from the downgoing slab; (7) deep earthquakes in slab under compres-

13.2 Subduction zones and their surroundings

The mechanical strength of lithospheric material is of particular importance in
subduction zones. We have seen in the analysis of forces acting on the lithosphere
in Section 12.3.2 that slab-pull is a major component of the system. The mass
circulation achieved by the recirculation of cold former lithosphere into the deeper
mantle forms a very important part of the mantle convection process, which we
address in Chapter 14. Here we concentrate on the way in which subducted material
enters the mantle and the consequent ¬‚ow ¬elds in their vicinity.
Figure 13.7 displays the general features associated with subduction in an
oceanic environment, with a deep ocean trench where the subducted lithosphere
enters the mantle and mantle ¬‚ow above and below the plate. Earthquake activity
generally occurs within the subducted plate except for great thrust events such as
the 1960 Chile (Mw 9.6) and 2004 Sumatra“Indonesia (Mw 9.3) earthquakes which
propagate along the interface between the subducting and overriding plate. Most
subduction is accompanied by volcanism in an island arc lying back somewhat
13.2 Subduction zones and their surroundings 309

from the trench on the overriding plate, linked to volatiles released from the former
oceanic lithosphere.

13.2.1 Con¬guration of subduction zones
The dominant behaviour in subduction is of a slab of oceanic lithosphere bending
at the entry into the mantle, and then descending to depth with a dip typically
in the range 35“60—¦ (Figure 13.7). The particular circumstances depend on the
relative motion of the subducting and over-riding plates. Although the dominant
behaviour can be extracted from two-dimensional models in which subduction
occurs perpendicular to the surface expression of the plate boundary, many systems
have a signi¬cant oblique component. For example, at the northern end of the
Sumatra“Andaman arc, along which the great 2004 earthquake (Mw 9.3) occurred,
the plate motion is only at an angle of 10—¦ to the boundary. The earthquake initiated
thrusting that propagated northward along the arc for some 1300 km (see Figure
11.26), with variations in character associated with changes in the strike of the
plate boundary and consequent modi¬cations of the physical state of the upper
plate surface.
Considerable variability can occur in the geometry of subduction, even for a
single plate boundary. Figure 13.8 illustrates the morphology of the subduction
zone in the Izu“Bonin“Mariana arc to the south of Japan, where the old Paci¬c plate
is subducting beneath the younger Philippine plate in a purely oceanic environment.
The shape of the subducted slab is inferred from seismic tomography studies. A
very substantial change in the geometry of the deeper part of the slab is inferred to
occur between the Izu“Bonin and Mariana segments. The near-horizontal segment
in the mantle transitions zone for the Izu-Bonin arc may well be a consequence


lat depth

Figure 13.8. Interpreted morphology and geometry of the subducting Paci¬c slab beneath
the Izu“Bonin“Mariana arc from P wave tomography images (Miller et al., 2005). The hole
corresponds to a region where there is change in seismic properties. Earthquakes from the
NEIC catalogue for events from 1967“1995 are superimposed on the slab con¬guration.
310 The In¬‚uence of Rheology: Asthenosphere to the Deep Mantle

of rapid slab roll-back in the evolution of the plate boundary. The change to the
much steeper dips in the Mariana arc occur where the Marcus“Nesker ridge enters
the subduction system, with at present the Ogasawara plateau at the trench. There
is a zone of much reduced contrast in P wavespeed at depth, that gives rise to the
apparent ˜hole™ in Figure 13.8. This zone may indicate the development of a tear
in the slab linked to the oblique component of subduction since it is marked by
pronounced earthquake activity.

Age Boundary
West Sumatra Bali Isl.



-0.5% 0.5%

Figure 13.9. Variations in the relative amplitude of bulk-sound speed and shear wavespeed
anomalies along the subducted plate in the Java“Sumatra subduction zone (Gorbatov &
Kennett, 2003). The cells with high-velocity anomalies associated with the subducted plate
were isolated and averages of the difference between the wave-speed variations from the
AK135 reference model were made through the thickness of the plate. The bulk-sound
speed values are subtracted from the shear wave-speed values, and the results are then pro-
jected on the vertical section. Thus large bulk-sound velocity anomalies appear as negative
values (patterned grey tone), while high shear wave speed is presented as positive anomaly
(solid grey tone). The points where the age of the sea¬‚oor changes are indicated and show
a strong correlation with the changes in the balance between the two wave types. Solid
black dots are earthquake hypocentres.

Such segmentation of the dip of a plate boundary is not an isolated occurrence.
For example, the Paci¬c plate descending beneath southern Chile has a number of
segments that progressively lead to the shallowing of dip to the north.
The physical properties of the lithosphere entering the subduction zone for a
single plate can also vary as a consequence of the age of the plate and its prior
history. Figure 13.9 shows such variability in the subducted Indo-Australian
plate along the Indonesian arc. The image is extracted from seismic tomography
undertaken with simultaneous inversion of P and S wave arrival times (Gorbatov
& Kennett, 2003), and indicates that contrasts in physical properties can persist
to signi¬cant depth. The dominance of shear wavespeed anomalies occurs for
lithosphere older than 85 Ma as it enters the trench, such oceanic lithosphere may
well have reached thermal maturity.
13.2 Subduction zones and their surroundings 311

These examples illustrate that the subduction process carries its own complexity
linked to the history of plate motions, so that there can be variations in the physical
properties of the subducted slab with segmentation on scales of 500 km or so.
Nevertheless, we can capture the dominant behaviour associated with the descent
of cool, dense material into the mantle with simple two-dimensional conceptual

13.2.2 Flow around the slab
Many of the characteristics of the ¬‚ow around a subduction zone can be represented
via a simple mechanical model (Figure 13.10) in which the stationary overriding
plate and the subducting plate both have thickness a, with convergence at speed
Vs. The subducting plate descends with dip δs, and the relative motion of the two
plates drives ¬‚ow in the mantle wedge above the subducting plate. The motion of
the subducting plate also leads to an underside ¬‚ow.
The ¬‚ow of mantle material outside the cold slabs is represented through the
behaviour of an incompressible viscous ¬‚uid, in a steady state, so the velocity of
the ¬‚uid must satisfy
·∇2v + ρ∇• = ∇p, with ∇ · v = 0. (13.2.1)
Here we have neglected spatial variations in viscosity associated with variation
in temperature, • is the gravitational potential and p is the pressure. For this
incompressible viscous ¬‚uid the vorticity = ∇ — v satis¬es ∇2 = 0.
The convergence of many subduction zones includes an oblique component, but
we can capture the main behaviour with the two-dimensional approximation that

δ vs
Over-riding Plate

Mantle Wedge

Figure 13.10. Streamlines for viscous ¬‚ow around a subducting slab, driven by the motion
of the sinking slab and its attached plate, using a corner ¬‚ow representation neglecting any
thermal convection effects. The lithosphere is taken as 50 km thick and the dip angle δs is
45—¦ .
312 The In¬‚uence of Rheology: Asthenosphere to the Deep Mantle

convergence is normal to the overriding plate. As in Section 7.4 we introduce a
stream function ψ such that the velocity ¬eld v in the ¬‚uid is expressed as
∇4ψ = 0,
v = ∇ — ψ ^⊥
e with (13.2.2)
where the unit vector ^⊥ lies perpendicular to the plane of convergence. It is
convenient to adopt cylindrical polar coordinates in each of the wedges so that
1 ‚ψ ‚ψ
v = vr, vθ, 0 = ’ , ,0 . (13.2.3)
r ‚θ ‚r
We want to ¬nd velocity ¬elds that have ¬xed values on radial vectors from the
wedge corners and seek solutions in the form (Batchelor, 1967)
ψ(r, θ) = rf(θ), (13.2.4)
so that from (13.2.2) we require
‚4f ‚2f
+2 2 +f = 0, (13.2.5)
r3 ‚θ4 ‚θ
with general solutions of the form
f(θ) = A sin θ + B cos θ + Cθ sin θ + Dθ cos θ; (13.2.6)
the constants A, B, C, D have to be found from the speci¬c boundary conditions.
In the mantle wedge zone we require
v = 0 on θ = 0; vr = Vs, vθ = 0 on θ = δs, (13.2.7)
where we measure θ anticlockwise from the base of the overriding plate. The
conditions are satis¬ed when (McKenzie, 1969)
rVs[(δs ’ θ) sin δs sin θ ’ θδs sin(δs ’ θ)]
ψw = , (13.2.8)
δ2 ’ sin2 δs
the associated shearing stress
· 2Vs· [δs cos(δs ’ θ) ’ sin δs cos θ]
σw = +f = ; (13.2.9)

δ2 ’ sin2 δs
r r s
the singularity at r = 0 is associated with the arti¬cial conditions at the wedge
corner where a sharp change in direction is required.
The stresses above and below the subducting plate act to oppose the motion
of the lithosphere and thus work has to be done to keep the plate moving. The
mechanical energy is then converted into heat by viscous dissipation in the mantle.
As demonstrated by McKenzie (1969), the stress heating due to viscous dissipation
in the ¬‚ow Hw = σ2 /·, and so the contours of the stress heating function are

the same as those for the shearing stress. The high stresses near the wedge corner
would lead to a modi¬cation of the viscous constitutive relation, but the solutions
(13.2.8) and (13.2.9) provide a reasonable representation for r > 50 km.
13.2 Subduction zones and their surroundings 313

δ vs
Over-riding Plate
2.5 20
Mantle Wedge




Figure 13.11. Shear stresses in MPa associated with the ¬‚ow in Figure 13.10, for a conver-
gence rate of 100 mm/yr and mantle viscosity of 3—1020 Pa s. The forces act counter to the
descent of the slab material and are larger beneath the subducting plate than in the mantle

For the corner ¬‚ow beneath the subducting plate, the boundary conditions are

vr = ’Vs, vθ = 0 on θ = 0; vr = Vs, vθ = 0 on θ = δs, (13.2.10)

with θ measured clockwise from the horizontal base of the subducting plate. The
corresponding stream function is

rVs[(π ’ δs ’ θ) sin θ + θ sin(δs + θ)]
ψb = ’ , (13.2.11)
π ’ δs + sin δs

with associated shear stress

2Vs· [cos δs + cos(δs + θ)]
σb = ’ . (13.2.12)

r π ’ δs + sin δs

The streamlines of constant ψ and the contours of shear stress can be mapped out
by sweeping across the relevant span of angles. The resulting patterns are illustrated
in Figures 13.10, 13.11 for a dip angle of 45—¦ , a convergence velocity of 100 mm/yr
and a lithosphere thickness of 50 km. The mantle viscosity in Figure 13.11 is taken
as 3—1020 Pa s both in the mantle wedge and below the subducting plate.
The streamline pattern in Figure 13.10 shows that the sinking slab drags ¬‚uid
down as it descends, as might be expected. Yet, in the mantle wedge, the path of
much of the ¬‚uid heads upwards towards the corner before turning to descend with
the slab. The stress ¬eld is contoured in Figure 13.11 and we see that the stresses
beneath the subducting slab are much more concentrated than in the mantle wedge.
314 The In¬‚uence of Rheology: Asthenosphere to the Deep Mantle

13.2.3 Temperatures in and around the subducting slab
The descending subducted slab carries with it the temperature distribution imposed
before subduction. This idea forms the basis of the conveyor belt approach
employed by Frohlich (2006) to characterise the temperature distribution in terms
of the age of the oceanic lithosphere at the entry to the trench tage with an additional
time component in subduction tsub = h/Vz, where h is the current depth and Vz
is the component of the velocity of the subducting slab in the vertical direction
(Vz = Vs sin δs). Frohlich employs the thermal half-space model for the oceanic
lithosphere discussed in Section 12.2.1, and switches from an oceanic boundary
condition to a mantle boundary condition once the slab starts to descend. This
highly simpli¬ed approach, based on conduction processes alone, leads to an
estimate for the temperature at a distance z measured perpendicular to the top of
the slab of
z z
¯ ¯
T (¯) = TM + (TO ’ TM) erf √ ’ erf √
z , (13.2.13)
2 κtsub
2 κttot
where TM is the mantle temperature, TO is the oceanic temperature, and the total
ttot = tage + tsub = tage + h/Vz. (13.2.14)
The coolest temperature in the slab occurs at
h h ¦
zc = 2κtage 1+ ln 1 + ,
¯ (13.2.15)
¦ ¦ h
where the thermal parameter for the slab ¦ = tage Vz. The temperature contrast
between this coldest point (Tc) and the mantle temperature (TM) is controlled by
¦/h, and thus the residence times in contact with ocean and mantle since
= . (13.2.16)
h tsub
The temperature contrast is approximately linear in log(¦/h) (Frohlich, 2006)
Tc ’ TM ¦
≈ 0.18 + 0.31 log , (13.2.17)
TO ’ TM h
a result which can be applied even if the slab does not have a simple geometry.
Compared with numerical estimates of the temperature distribution (Kirby et al.,
1996) this simple model including just conduction effects tends to overestimate
temperature in the slab at shallower depths and to underestimate below 300 km.
A comparable approach can also be built from the plate thermal model (12.2.14)
as used by McKenzie (1969), but the series representation is much less amenable
to analytic simpli¬cation.
The previous results ignore the details of the behaviour of the ¬‚uid in the
mantle wedge with imposition of a constant temperature condition on the top of
the subducting slab. The corner ¬‚ow will in fact maintain an advective boundary
13.2 Subduction zones and their surroundings 315

layer on the top of the slab that will control the temperature at the interface
between the slab and the wedge (England & Wilkins, 2004); this temperature Ts
is depth-dependent.
Although the dominant direction of ¬‚ow in the mantle wedge above the slab is
down the subducting plate, there is a component v⊥ perpendicular to the slab that
can advect heat. In terms of the coordinate z perpendicular to the top of the slab,
the temperature equation is

d2T dT
κH 2 + v⊥(¯)
z = 0, (13.2.18)
z z
d¯ d¯
in z < 0 above the top of the slab. The thickness of the advective boundary layer
La is determined by the requirement that the time scale for the diffusion of heat
through this layer, L2 /(π2κH), is comparable to the timescale for heat to be carried
the same distance by the ¬‚ow, La/U, where U is a characteristic velocity. Thus we

L2 U π2κH
∼ 1, La ∼ ,
so (13.2.19)
2κ L
π Ha U
and so the P´ clet number introduced in Section 7.2.4 is given by
∼ π2
Pe = (13.2.20)
England & Wilkins (2004) examine the case of corner ¬‚ow in the mantle wedge
in some detail and derive approximate expressions for the thickness of the boundary
layer in terms of the dip angle δs and convergence speed Vs. For the corner ¬‚ow in
the mantle wedge

La ≈ , (13.2.21)
2 πVsξ
where r is the radial distance from the wedge corner and ξ is a slowly varying
function of δs with a value around 0.55. Similar power-law dependences could
be expected in the more complex real situation. For typical slab parameters, the
advective boundary layer is 10“20 km thick, with a temperature contrast of more
than 600—¦ C between the wedge and the slab, corresponding to a heat ¬‚ux of around
0.1“0.2 W m’2 .
The temperature in the advective layer is approximately
T (¯) = Ts + (Tw ’ Ts)erf
z , (13.2.22)
where Ts is the temperature at the top of the slab and Tw the maximum temperature
in the mantle wedge on a perpendicular pro¬le to the slab. At the interface between
the mantle wedge and the subducting plate, conservation of heat ¬‚ux requires the
316 The In¬‚uence of Rheology: Asthenosphere to the Deep Mantle

temperature gradients in the slab and in the advective layer to be equal (under the
assumption of common thermal conductivity), so that
√ √ ’1
πLaTM a πLa
1+ √
Ts ≈ Tw + ,
erf (13.2.23)
2a 2 κtsub
2 κr/Vs
where a is the plate thickness.
These simple models are in good general agreement with the results of full
numerical modelling (England & Wilkins, 2004). Dissipation can also make a
signi¬cant contribution to heating near the wedge corner. Once the slab has
penetrated 100“200 km into the mantle, the dominant in¬‚uence on temperature
structure comes from advection by ¬‚ow in the mantle wedge.
The corner ¬‚ow model describes the dominant mantle motions but does not
represent the undoubtedly important movement of volatiles from the slab to the
surface where they help to create the volcanic arc. This secondary process has
much smaller mass-transport. The temperature at the top of the subducting plate
does not vary strongly with subduction parameters and so the negative correlation
of the distance of the volcanic front from the trench with convergence rate and slab
dip suggests that the positions of the volcanoes in the island arc are controlled by
a strongly temperature-dependent process taking place in the wedge (England &
Wilkins, 2004).
The situation on the underside of the subduction zones may well be somewhat
more complex than in the simple corner ¬‚ow models. For a number of subduction
zones, seismic tomography results reveal a concentrated zone of lowered velocities
beneath the subducted slab at depths around 400 km that is not easily explained,
but might involve some counter ¬‚ow.
Although two-dimensional numerical simulations of the evolution of a
subduction zone can provide a good indication of behaviour, the details have a
strong dependence on the rheology assigned to the different materials. Frequently
a visco-plastic formulation is employed for the plate material, but this does not
provide an adequate description of fault failure and the initiation of subduction.
Most current subduction zones have some component of oblique convergence, and
in some cases, such as the Andaman Arc, the oblique component can be dominant.
This means that three-dimensional effects cannot be ignored in the understanding
of subduction systems.

13.2.4 Subduction and orogeny
Since the recognition of current-day subduction zones and their related features,
plate tectonic concepts have been extensively employed in the explanation of past
patterns of orogeny. The gravitational drag of the slab has a tendency to induce a
retreat of the hinge where the slab bends. Such slab retreat appears to have played
a signi¬cant role in some regions where there has been major realignment of plate
13.2 Subduction zones and their surroundings 317

35.1 ± 0.9 Ma

Pressure [GPa]


Depth [km]
3.4 cm/yr


32.9 ± 0.9 Ma
1.6 cm/yr
31.8 ± 0.5 Ma
0.5 cm/yr

29.9 ± 1.4 Ma
200 400 600 800 28 30 32 34 36
T [°C] Age [Ma]

Figure 13.12. Pressure“temperature and depth“age paths estimated for ultra-high-pressure
rocks in the Dora“Maira unit in the Western European Alps. The age estimates were
obtained using in situ U-Th-Pb dating on titanite samples. The pressure and temperatures
were estimated using a variety of different methods. [Courtesy of D. Rubatto & J.

boundaries, such as the Izu“Bonin region near Japan and the Tonga“Kermadec
arc. For the trench line now near Tonga, the line of subduction has pivoted about
the North Island of New Zealand over the last 40 Ma. The hinge retreat has led
to a rotation of the Fiji Island block by 90—¦ , determined from palaeomagnetic
measurements. The rapid slab retreat has enabled asthenospheric upwelling to
produce back-arc spreading centres in the Lau basin behind the Tonga volcanic arc.
Images from seismic tomography indicate that in regions where rapid slab roll-back
has occurred, ˜stagnant™ slabs lie horizontally on top of the 660 km discontinuity
with a connection to the present subduction. Such a situation can be seen for the
Izu“Bonin arc in Figure 13.8.
Within the generally convergent regimes associated with subduction and the
consumption of oceanic lithosphere by passage into the mantle, offers the
possibility of an extensional regime. Some time can be expected to elapse after
the initiation of subduction before roll-back becomes signi¬cant. The insertion of
extensional episodes in a convergent environment helps to reconcile some of the
puzzling geological structures in former orogens.
Within many mountain belts there are isolated occurrences of ultra-high-pressure
minerals whose characteristics require that material has been carried to great depth,
318 The In¬‚uence of Rheology: Asthenosphere to the Deep Mantle

presumably as a result of subduction, and then brought rather rapidly to the surface
so that the high-pressure phases survive. The example in Figure 13.12 (Rubatto
& Hermann, 2001) is based on analysis of nodules from the Dora“Maira unit
in the Western Alps. The presence of coesite, a high-pressure form of quartz,
indicates that material has reached depths of more than 100 km. The rates of
inferred exhumation are of the order of cm/yr and so not much slower than the
subduction process that carried the material to depth. Numerical simulations of the
complex deformation of the over-riding plate in the course of subduction indicate
that rather complex material trajectories can arise, with substantial variations in the
paths taken by points that are initially rather close (see, e.g., Gerya, Stockhert &
Perchuk, 2002).
Such numerical simulations require the coupling of the mechanical behaviour,
through the equations of motion and continuity conditions, with the heat ¬‚ow
equations. A very important role is played by the assumptions made about the
appropriate rheologies, particularly in the brittle regime where the distinction
between wet and dry rock has a marked in¬‚uence on behaviour. For a wet rock the
brittle strength decreases with depth, but for a dry rock the strength will increase.
Diffusion creep is expected to be signi¬cant in the subduction zone environment
and the interplay between temperature and grain size in (9.3.4) leads to only modest
variations in viscosity. An effective rheology can be found by combining creep and
quasi-brittle behaviour controlled by a yield stress σy. A ˜brittle™ viscosity ·b can
be de¬ned as

·b = 1 σy/( ™ II)1/2 with σy = (M1plith + M2)(1 ’ »f), (13.2.24)

where ™ II is the second invariant of the strain rate tensor 1 ™ ij ™ ij. The total effective
viscosity ·e is then determined by the relative size of the yield stress and the stress
associated with creep (Gerya et al., 2002)

2( ™ II)1/2 < σy,
·creep when
·e = (13.2.25)
2( ™ II)1/2
·b > σy.
A further complexity comes from the potential release of ¬‚uid from the
subducted slab with migration of a hydration front into the mantle wedge
introducing weak serpentenized mantle. Many of the necessary physical properties
are poorly constrained by conventional laboratory experiments, and in consequence
there is considerable latitude available in the choice of material parameters that can
lead to signi¬cant differences in the overall mechanical behaviour.
Continent“continent collisions can produce major mountain chains that are
susceptible to the in¬‚uence of erosion. The elevation of mountains can have the
effect of changing weather patterns, e.g., the rise of the Himalayas appears to
have caused a major modi¬cation of the monsoon system around 8 Ma ago. A
number of authors have demonstrated that erosion has the potential to modify the
deformation pattern and material trajectories in a collisional environment. For
13.3 The in¬‚uence of phase transitions 319

example, Pysklywec (2006) provides a number of examples of the way in which
surface erosion patterns can in¬‚uence the con¬guration of the lithosphere, and
indeed the nature of an orogeny, when buoyant continental material is brought into
the subduction environment.

13.3 The in¬‚uence of phase transitions
Evidence for the composition of the mantle comes from xenoliths, from slices of
mantle rocks emplaced at the surface in a few locations, from the composition
of mantle-derived melts and from cosmochemical arguments, supplemented by
geophysical results. The general consensus is that the mantle resembles a form
of peridotite with around 60 per cent olivine by volume. Most observations are
consistent with a single chemical composition through the bulk of the mantle, for
which the physical properties are modi¬ed through the progressive effect of high
pressure phase transitions over the depth range from 300“800 km (Figure 13.13).



P [GPa]

Figure 13.13. Depth-varying


proportions of different mineral
phases for a ˜pyrolite™ model of the
mantle. The phases represented are:
Depth [km]

(±) olivine,
(β) wadsleyite,
(γ) ringwoodite,
(opx) orthopyroxene,
(cpx) clinopyroxene,
(gt-mj) garnet-majorite,
(mw) magnesiow¨ stite,
(Me,Fe-pv) ferromagnesian silicate

perovskite, and

(Ca-pv) calcium silicate perovskite.
[Courtesy of C. Bina].


0 20 40 60 80 100
Volume Fraction [%]

The olivine component (±) transforms to wadsleyite (β) and then to ringwoodite
(γ-spinel) with increasing pressure, ¬nally ending up as a mixture of silicate
320 The In¬‚uence of Rheology: Asthenosphere to the Deep Mantle


phase A


phase B
p z

Figure 13.14. Where the transition from phase A to phase B does not take place
immediately there will be a gradient in the physical parameters linking the behaviour
associated with the two individual phases.

perovskite (pv) and magnesiow¨ stite (mw - also known as ferropericlase). The
remaining components of the peridotite are orthopyroxene (opx), clinopyroxene
(cpx) and garnet (gt). These components undergo a more gradual evolution
under the in¬‚uence of pressure as the pyroxenes dissolve in the garnet. The
resulting garnet-majorite solid solution (gt-mj) ultimately also transforms to silicate
perovskite, but this transition extends to greater pressure, and hence depth, than for
the olivine component.
Although a simple picture of a phase transition has the transformation taking
place at a line boundary in pressure“temperature space, depending on composition
some degree of coexistence of the two phases may be possible (Figure 13.14).
In this case there will not be a simple discontinuity in physical properties such
as seismic velocity across the phase transition, but rather a gradient zone whose
properties will depend on the way in which the transformation takes place. The
situation with regard to the depth dependence of physical properties can be further
complicated by several constituents of the mineral assemblage in the mantle having
transitions in a similar pressure interval.
When we wish to examine the behaviour of materials as a function of
temperature (T ) and pressure (p), the most convenient thermodynamic variable is
the Gibbs free-energy. The speci¬c Gibbs free-energy G satis¬es the ¬rst law of
thermodynamics in the form

dG(T, p) = ’S dT + v dp. (13.3.1)

where v is the speci¬c volume.
We will consider a simple univariant phase transformation for which there is no
mixed domain. Then, across the transition between the two phases A,B there can be
13.3 The in¬‚uence of phase transitions 321



β liquid


pv + mw

Figure 13.15. Schematic representation of the pressure“temperature relations for olivine
in the upper mantle, with ultimate transformation to perovskite and magnesiowustite (pv +
mw). The Clapeyron slope of the ± ’ β transition associated with the 410 km seismic
discontinuity is positive; whereas for the γ-spinel to perovskite transition the Clapeyron
slope is negative. The β ’ γ transition is not simple and as indicated takes place over a
range of pressure at a particular temperature. The dashed grey line indicates typical mantle
conditions; hotter or cooler trajectories will encounter the phase transitions at different

no net change in G. Thus, considering the situation at a particular p, T combination
on the boundary, we require
”GB,A = (vB ’ vA)dp ’ (SB ’ SA)dT = 0. (13.3.2)
The slope of the phase boundary, the Clapeyron slope γc, in pressure-temperature
space is given by
”S L
γc = = = = , (13.3.3)
”v T ”v T ”ρ
dT bdy

where L is the latent heat associated with the transition, L = T”S. For the upper
mantle minerals the higher pressure phase is always denser ”ρ > 0, and so ”v
is always negative. Latent heat may be absorbed or evolved depending on the
particular transition.
For most phase transitions, the transformation pressure increases with increasing
temperature, this means that the Clapeyron slope γc is positive and the reaction is
exothermic L < 0. The ± ’ β transition near 410 km depth has this character
(Figure 13.15). However, in the γ-spinel to perovskite transition associated with
the 660 km seismic discontinuity, the change in coordination of the silicon atoms
322 The In¬‚uence of Rheology: Asthenosphere to the Deep Mantle

from 4 to 6 is accompanied by an increase in the length of certain atomic bonds
with a consequence of an increase in entropy. The latent heat is then positive for
this endothermic reaction, and so the Clapeyron slope has to be negative (Figure
A further complication around 660 km depth is transformation of garnet to
perovskite that, depending on the temperature regime, can occur over an extended
depth interval. The combination of the γ-spinel and garnet transitions leads to
complex changes in the physical properties, so that, e.g., the seismic re¬‚ectivity of
the 660 km transition can be expected to be strongly temperature-dependent.
Impact of phase transitions on subduction
The material in a descending subducted slab interacts with the mineralogical
transitions, but because the thermal state of the slab differs from the surrounding
mantle, phase changes within the slab will occur at a different depth. The
displacement of the phase transition is governed by the Clapeyron equation
(13.3.3), under the pressure and temperature conditions within the slab. The
Clapeyron relation can be recast as a representation for the depth change
= , (13.3.4)
with the sign of the de¬‚ection controlled by the Clapeyron slope γc. The
exothermic ± ’ β transition, with positive Clapeyron slope, will be shallower
in the cool core of the slab compared with external conditions. In contrast, the
endothermic γ-spinel to perovskite transition would be expected to produce an
enhanced depth for the phase boundary in the colder slab. The 410 km discontinuity
tends to enhance the subduction potential since the presence of denser material at
shallower depth increases the downward forces, whilst the 660 km transition tends
to impede subduction. Careful seismological studies of slabs are consistent with
the predictions from the thermodynamic arguments, but the detailed shapes of the
de¬‚ected boundaries are dif¬cult to estimate in this heterogeneous environment.
The Clapeyron theory describes an equilibrium phase transition, but the
conditions within the slab are conducive to the preservation of metastable states.
Metastable material can be expected if the high-pressure phase is produced by
nucleation on grain boundaries with ultimate growth to consume the crystals of
the lower-pressure phase. Within the coldest slabs it is likely that a wedge of
metastable olivine can persist to depths approaching 660 km in a zone concentrated
towards the upper slab interface (Stein & Rubie, 1999). It has been suggested
that phase changes within this metastable zone may act as the trigger for deep
earthquakes via transformational faulting in which slip occurs along thin shear
zones where metastable olivine transforms to denser spinel. An alternative model
is that decomposition of hydrated minerals leads to the presence of water that can
cause reduction of effective stresses on faults in the slab. Deep earthquakes have
the characteristics of slip on faults that in some cases, such as the 1994 Mw 8.4
deep Bolivian earthquake, appear to cut through the whole slab.
13.4 The deeper mantle 323

Mineral physics studies in recent years have demonstrated that many nominally
anhydrous silicate phases can carry signi¬cant amounts of water to depth. As a
result the hydration potential adds another variable in the description of behaviour,
since even small amounts of free water lead to a major reduction of the effective
viscosity of materials.

13.4 The deeper mantle
Inferences on the long-term mantle viscosity have been made from two major
classes of observations: (a) patterns of sea-level change related to the adjustment
of the Earth to the effects of the ice loading and unloading during the glacial cycles
of the late Pleistocene glacial cycles, dominated by isostatic adjustment; (b) a suite
of surface geophysical observables linked to convective processes, such as long
wavelength free-air gravity anomalies, the horizontal divergence of tectonic plate
motions, and dynamic surface topography.
Models that satisfy both classes of data are characterised by a signi¬cant
increase in viscosity within the lower mantle relative to the upper mantle by
a factor of 30“100. The position of the increase has been a matter of some
debate, and although it is often modelled as a sharp change the data would be
entirely compatible with a transition spread over some hundreds of kilometres.
For example, there could be a ramp in viscosity from the mean upper mantle
of around 3—1020 Pa s to 1022 Pa s below 1000 km depth. Some models such
as that of Mitrovica & Forte (2004) invoke a thin layer of signi¬cantly reduced
viscosity just above the 660 km discontinuity and thereby reduce the absolute value
of viscosity required in the lower mantle. Such a zone would have the effect of
partially decoupling ¬‚ow patterns between the upper and lower mantle.
What might be happening to give such a signi¬cant change in viscosity within
the upper part of the lower mantle? Calculations of the properties of the mineral
assemblage for the mantle suggest that the proportion of garnet should diminish
with depth and drop to zero by 900“1000 km leaving just perovskites such as
(Mg,Fe)SiO3 and ferropericlase (Mg,Fe)O. A strong gradient in viscosity could
therefore arise from the diminishing in¬‚uence of the more deformable garnets.

13.4.1 Viscosity variations in the mantle and the geoid
A ¬‚uid body under rotation assumes the shape of an oblate spheroid as a result
of the balance between self-gravitation and the centrifugal effects of rotation. The
internal variations in density mean that the ¬gure of the Earth departs from a solid
of revolution. Density anomalies in the mantle contribute to the gravity ¬eld of
the Earth in two ways: a local mass anomaly (excess or de¬ciency) has a direct
effect on the gravitational ¬eld, but also causes ¬‚ow in the mantle that leads to
distortion of the surface and the core“mantle boundary as well as any internal
324 The In¬‚uence of Rheology: Asthenosphere to the Deep Mantle

Figure 13.16. Con¬guration of the long-wavelength components of the geoid, displayed
as an outer surface with an exaggerated scale of deformation from a spherical shell. The
inner spherical surface is present to provide a geographical reference. [Courtesy of M.

density interfaces. The ¬‚ow-induced distortions lead to additional mass anomalies
that can contribute to the gravitational ¬eld and its potential.
The surface of constant gravitational potential closest to mean sea level is termed
the geoid (Figure 13.16). Across the oceans the geoid lies close to the actual
surface. In the continents the geoid has no direct relation to topography, but
would be close to the water level cut across a continent to connect two oceans
directly. Because gravitational potential decays as R’1 with distance R from a
mass anomaly, the geoid is much more sensitive to the mass con¬guration at depth
than gravity with an R’2 dependence.
The topographic variation of the geoid from the rotational ellipsoid is of the order
of about 100 metres with a profound depression to the south of India (Figure 13.16).
The large-scale variations in gravitational potential • can be well represented
through a spherical harmonic expansion with the origin placed at the centre of mass

∞ l
GMe re l
[Cm cos(mφ) + Sm sin(mφ)] ,(13.4.1)
•= 1+ Pl (cos θ) l l
r r
l=2 m=0

where G is the gravitational constant and Me is the mass of the Earth. The l = 0
term is extracted from (13.4.1) as the central mass contribution, and the l = 1 term
13.4 The deeper mantle 325

vanishes with the choice of origin. The number of contributions to the expansion
increases rapidly with increasing degree l.
The oblate ellipsoid associated with rotation, with mean radius re and ellipticity
e, takes the form
(1 ’ e)2
≈ (1 + 1 e) + 2 e P2(cos θ),
= (13.4.2)
3 3
(1 ’ e)2 sin2 θ + cos2 θ

the higher-order spherical harmonic contributions decrease steadily but the P2 term
alone does not represent an ellipsoidal surface.
The shape of the geoid depends on the mass distribution within the Earth and the
in¬‚uence of self-gravitation and ¬‚ow that tend to deform the Earth™s surface and
the core“mantle boundary as well as any internal boundary surface. The deviations
of the geoid from the reference ellipsoid depends on instantaneous force balance
rather than ¬‚ow rates and hence are sensitive to the variations of viscosity rather
than the absolute values.

Dynamic compensation
Following Hager (1984) we consider a two-dimensional ¬‚ow system as illustrated
in Figure 13.17(a) consisting of a plane layer of unit depth with uniform viscosity
and density overlying an inviscid ¬‚uid half-space of greater density. A surface
density contrast ”ρm is placed at the midpoint of the layer with an amplitude

(a) (b)

Figure 13.17. Illustration of the components of the geoid anomaly due to a density contrast
in the form of cosine bell at the centre of a layer: (a) a layer with uniform viscosity ·, (b)
with the bottom half of the layer having a viscosity 30 times larger than that of the upper
half. In each case the total geoid anomaly (heavy solid line) is the sum of contributions
from the density contrast (light solid line), the dynamic deformation of the upper boundary
(long dashes) and the lower boundary (short dashes). In case (a) the total geoid anomaly is
negative for a positive density contrast, but in case (b) the corresponding geoid anomaly is
326 The In¬‚uence of Rheology: Asthenosphere to the Deep Mantle

speci¬ed by a cosine bell, the effect of this anomaly is to lead to an equal mass
displacement at both upper and lower boundaries.
At the top and bottom of the viscous layer the normal and shear tractions will
be continuous. The introduction of a sudden density contrast will lead to the
boundaries deforming until a steady state is reached. A sinking positive density
contrast will push material downwards and pull material from behind to create a
depression of the upper surface. The surface density contrast due to the depression
”ρt then acts to set up a secondary ¬‚ow tending to ¬ll in the topography in a way
that is precisely analogous to post-glacial rebound, with the same time constant.
In the steady state the rate at which the sinking heavy blob pulls material away
from the surface is balanced precisely by the rate at which ¬‚ow tends to ¬ll in the
depression. A similar behaviour occurs for the bulge at the base of the layer with
associated equivalent surface density contrast ”ρb. ”ρt and ”ρb will have the
opposite sign to ”ρm.
The amplitude of the surface depression depends on the size and depth of the
density contrast and the thickness of the layer, and is unaffected by the value of the
viscosity in the layer. However, the time needed to reach steady state is proportional
to the viscosity. The displacement of the density contrast during this interval is
similar to the ¬nal surface deformation.
The total geoid anomaly is the sum of contributions from ”ρm, ”ρt, and ”ρb
as indicated at the top of Figure 13.17(a). The surface density anomaly ”ρt can
contribute at short spatial wavelengths because of its proximity to the measuring
point, but the in¬‚uences of ”ρm and ”ρb will shift to longer wavelengths due to
the geometric decay of gravitational potential from a mass anomaly. In a uniform
viscosity ¬‚uid the total geoid effect from an anomaly with an increase in density
is negative, and such an effect does not explain the association of geoid highs with
regions where active subduction is occurring.
However, a modest change in the con¬guration produces the desired result,
Figure 13.17(b). We introduce now a much higher viscosity in the bottom half of
the layer, with a contrast of a factor of 30 with the upper half, For the same positive
density contrast much less deformation occurs at the upper boundary, but the effect
on the lower boundary is ampli¬ed compared with the uniform viscosity case. The
sinking of the density anomaly is slowed down and so the surface depression is
less, decreasing ”ρt; now more material needs to be moved at the lower boundary
to accommodate the motion and ”ρb is increased.
The total geoid anomaly is now positive, but only a fraction of that due to the
mass contrast ”ρm alone which is unchanged from the previous case. The increase
in the effect of ”ρt overcomes that from ”ρt, which is smaller than in the uniformly
viscous case.
These two models illustrate that the sign of the geoid anomaly for a given density
contrast is not immediately obvious, but depends on the viscosity variation with
depth. Both distance and viscous ¬‚ow act as ¬lters to determine the magnitude of
13.4 The deeper mantle 327

the total anomaly. The behaviour can be expected to be quite complex for features
such as the density contrasts associated with subduction zones.

Spectral representation
For spherical harmonic degree l and order m, the contribution to the geoid anomaly
δ•m can be found as an integral over the corresponding spherical harmonic load
contribution from density δρm in the form (Richards & Hager, 1984)
δ•m dr K(r, ·)δρm,
= (13.4.3)
l l
2l + 1 rc

where K(r, ·) is the geoid response kernel calculated from the deformation
produced by viscous ¬‚ow due to a surface density contrast at radius r, including
the boundary deformation contributions.
The issue then is how to ¬nd the requisite density contributions. A common
assumption is to use some class of scaling relation between seismic wavespeed
heterogeneity determined from seismic tomography and density variations. The
usual choice is based on a simple scaling from the shear wave heterogeneity
δρ δβ
= C(r) , (13.4.4)
ρ β
where C(r) is a scaling factor that is allowed to vary with radius. Such a relation
would be appropriate when the only effective in¬‚uence is temperature; however,
we can expect somewhat different behaviour when compositional variations come
into play. Consider the relative variation in shear wavespeed as a function of shear
modulus G, density ρ, temperature T and some compositional component X,
δβ 11 ‚G 1 ‚ρ 1 ‚G 1 ‚ρ
δT ’ 1 δT + 1 δX ’ 1
= 2 ρ ‚X δX. (13.4.5)
2G 2 ρ ‚T 2 G ‚X
β ‚T
We cannot expect a similar scaling of the in¬‚uences of temperature and
composition and so a simple relation as in (13.4.4) is likely to fail. There are
indications that such scaling works reasonably well for the longest-wavelength
features of heterogeneity where temperature may well dominate. However,
compositional effects can be expected to be localised with a strong effect for shorter
wavelength and also to be more prominent near the boundaries of the mantle.
A signi¬cant component of the geoid signal comes from subducted slabs and
once this has been removed using plausible slab models, based on tomographic
imaging and past plate con¬gurations, the remaining component is well described
using a two-layer mantle with a factor of 30 contrast in viscosity in the
neighbourhood of 1000 km deep. Such a change at depth is compatible with
the observations of glacial rebound, even though such data can be adequately
represented with hardly any change. The available data provide little support for


. 12
( 16)