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)

2hl

νb δh

z

x

h

ν0

νb

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)

e

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

ν0

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

‚x

‚˜

∇4ψb ’ βRa = 0 in ’ 1 ¤ z ¤ ’ 1 + δ, 1 ’ δ ¤ z ¤ 1 . (13.1.9)

2 2 2 2

‚x

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

2

ψ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

2

effect of the normal viscous stress in the interior layer acts on the less viscous layer

as a compensating pressure term at z = 1 ’ δ,

2

‚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 ’δ

2

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 ,

z

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

2

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

1

(13.1.17)

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

1

(13.1.18)

‚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

3

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

1

ψ0 = 2 C(x ’ l ) + cos(amx){Am cosh(amz)+Bm sinh(amz)} ,(13.1.21)

m=1

where

π

am = (2m ’ 1) . (13.1.22)

2l

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

¯ ¯

1

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

δ3 δ

m=1

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

m

δ3βa4 d ’ 3 cosh2(amd)

m

Gm(δ) = (13.1.24)

2

3

3βa3 [cosh(a d) sinh(a d) ’ a d] ’

δ 2 cosh (amd)

m m m

m

The expression (13.1.23) can be well approximated by

z3 z

¯ ¯

2 2

1

ψ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

2

(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

x

ζ = 1 ‚zψb z,

ξ= d˜ ‚zψb(˜ ),

x x ¯ (13.1.27)

2

’l

the temperature equation can be transformed into a diffusion equation

2

‚˜ 1‚ ˜

= 4 ‚ζ2 . (13.1.28)

‚ξ

A solution which satis¬es the transformed boundary conditions

˜ = ’1 ζ = 0, ˜ = 0 for ζ ’ ∞,

at (13.1.29)

2

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

2

1

‚ C(2 + G(a1) l’x

2

3

Nuloc = ’ ˜(ξ(x)) = . (13.1.31)

4

‚z πδ (2l ’ x)1/2

The regular Nusselt number is the average of (13.1.31) over the interval ’l ¤ x ¤ l

1

C(2 + G(a1) 2

1

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)

1

H(a1) = 2 ;(13.1.34)

32βa3(sinh(a d) cosh(a d) ’ a d) + 3 cosh2(a d)

πδ 1 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,

1/3

’l+

[2 + G(a1)]l

2/3

1

C= 2 Ra dx ˜ = Ra . (13.1.35)

4π[1 + H(a1)]2δ

’l’

Hence the Nusselt number can be expressed as

1/3

[2 + G(a1)]2l2

1/3

1

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

values.

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

1

Island 4

2

Overriding Plate Subducting

Arc

3

Plate

6

5

7

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-

sion.

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

long

lat

long

lat depth

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.

410

660

-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

models.

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

ab

Sl

g

tin

uc

bd

Su

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

e

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

1

+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

s

the associated shearing stress

‚2f

· 2Vs· [δs cos(δs ’ θ) ’ sin δs cos θ]

σw = +f = ; (13.2.9)

rθ

δ2 ’ sin2 δs

‚θ2

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

rθ

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

10

5

20

2.5 20

10

5

Mantle Wedge

10

2.5

5

ab

Sl

2.5

g

tin

uc

bd

Su

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

wedge.

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

now

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θ

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

time

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

tage

¦

= . (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

a

the same distance by the ¬‚ow, La/U, where U is a characteristic velocity. Thus we

require

L2 U π2κH

a

∼ 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

e

UL

∼ π2

Pe = (13.2.20)

κH

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

1/3

5κr2δs

√

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

’¯

z

T (¯) = Ts + (Tw ’ Ts)erf

z , (13.2.22)

La

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

4

120

35.1 ± 0.9 Ma

100

3

Pressure [GPa]

e

coesit

artz

qu

Depth [km]

80

3.4 cm/yr

2

60

40

32.9 ± 0.9 Ma

1

1.6 cm/yr

20

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.

Hermann.]

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)

2

where ™ II is the second invariant of the strain rate tensor 1 ™ ij ™ ij. The total effective

2

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.

when

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

200

opx

cpx

±

+

P [GPa]

400

Figure 13.13. Depth-varying

β

gt-mj

proportions of different mineral

phases for a ˜pyrolite™ model of the

mantle. The phases represented are:

γ

Depth [km]

600

(±) olivine,

(β) wadsleyite,

(γ) ringwoodite,

(opx) orthopyroxene,

29.2

800

(cpx) clinopyroxene,

(gt-mj) garnet-majorite,

(mw) magnesiow¨ stite,

u

(Me,Fe-pv) ferromagnesian silicate

Ca-pv

mw

perovskite, and

(Mg,Fe)-pv

38.2

1000

(Ca-pv) calcium silicate perovskite.

[Courtesy of C. Bina].

47.5

1200

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

v/·

T

phase A

phases

A+B

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

u

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

T

±

β liquid

γ

p

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

pressures.

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

ρ2L

”S L

dp

γ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

13.15).

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

γc

dz

= , (13.3.4)

ρg

dT

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.

Sandiford].

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

m

[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/2

(1 ’ e)2

r(θ)

≈ (1 + 1 e) + 2 e P2(cos θ),

= (13.4.2)

3 3

(1 ’ e)2 sin2 θ + cos2 θ

re

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

negative.

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

l

contribution from density δρm in the form (Richards & Hager, 1984)

l

re

4πGre

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