. 11
( 16)


Lithosphere FTS


Figure 12.8. Schematic representation of the plate boundary force system, with both diver-
gent and convergent margins

The various forces can be grouped as:
(i) Forces associated with divergent plate boundaries, associated with the
gravitational effects of the relative elevation of the mid-oceanic ridge system
(ridge push FRP ).
(ii) Forces arising from frictional resistance to motion at transform fault systems
(transform fault resistance FTF “ not shown in Figure 12.8).
(iii) Forces linked to the complex processes at convergent plate boundaries: the
negative buoyancy (heaviness) of the cooler subducting plate (slab pull FSP );
12.3 Plate boundaries and force systems 273

resistance of the mantle to slab penetration (slab resistance FSR ); forces associated
with the impact of plate convergence (collisional resistance FCR ; and tensional
forces in the over-riding plate resulting from the presence of the oceanic trench
(trench suction FTS )
(iv) Forces related to the interaction of the plate and the underlying mantle (mantle
drag FMD ). With the complexity and thickness variations of the continental
lithosphere such drag effects are likely to be enhanced for portions of plates with
continental material with an additional drag term FCD compared with the oceanic
portions of the plates.
In the original analysis Forsyth & Uyeda (1975) estimated the relative
importance of the various forces (assumed to be the same for all plates) under the
assumption that the sum of the torques acting on each plate must be zero. This
leads to a set of linear equations in terms of the sizes of the forces, modulated by
the con¬guration of the lithospheric boundaries. Some contributions are simply
an integral along the length of the margin, but for the collisional resistance for
convergent margins and the resistance at transform faults the direction of relative
motion of the plates has to be taken into consideration. The inversion of the
over-determined system of equations (36 simultaneous equations in 8 unknowns)
can be accomplished using a singular value decomposition. The conclusion is that
slab pull and slab resistance are an order of magnitude larger than the other forces,
but that the net pull of downgoing slabs on the surface plates, FSP + FSR , is of the
same order as the rest of the forces. This means that although slab pull is signi¬cant
there can be a distinct contribution from ridge push.
A direct estimate of the ridge push per unit length of the ridge axis can be made
by looking at the gravitational effect produced by the elevation of the ridge crest.
Ritcher & McKenzie (1978) derive the simple form
+ 1 we ,
FRP = gwe(ρm ’ ρw) 3L (12.3.1)

in terms of the thickness L of the cooled plate and the elevation of the ridge axis
above the equilibrium depth we. For a typical plate thickness of 85 km with a 3 km
ridge elevation we obtain the estimate of 2—1012 N m’1 . The ridge push acts on
the entire plate and not just as a line force at the ridge axis.
The slab pull force arises from the fact that the descending plate is cooler and
denser than its surroundings. An estimate of the force FSP at depth z can be obtained
from models of the temperature distribution in the descending plate, so that, e.g.,
8g±th CpTMρ2 L2R ’π2z ’π2(dUM ’ L)
FSP (z) ≈ ’ exp ,
π 2cpρmR 2cpρmR
for a plate of thickness L, with mantle temperature TM, upper mantle thickness
dUM , ±th the coef¬cient of thermal expansion, Cp the speci¬c heat and R the
thermal Reynolds number introduced in (12.2.14) expressed in terms of the rate of
274 Lithospheric Deformation

slab sinking. As discussed in Section 13.3 the elevation of the olivine to spinel
phase transition within the slab, compared with the normal depth close to 410
km outside, provides an extra downward force because of the presence of denser
material at shallower depth. The total contribution to FSP is around 1013 N m’1 .
The descent of the plate is resisted by frictional forces acting on the sides of the
plate and by the phase transition near 660 km. Seismic tomography gives evidence
of impediments to penetration with stagnant slabs residing in the transition zone,
particularly in regions where there has been considerable trench retreat. Most plates
with a considerable length of subduction margin move with a similar speed which
suggests that there is a terminal velocity for the plate achieved by the near balance
of slab-pull and resistive forces.
The actual forces acting on any particular plate are likely to be more complex
than indicated in Figure 12.8, because negative buoyancy (heaviness) of the
subducting plate is age-dependent and the interaction of plates means that
convergent margins may be in retreat or advance with a consequent change of the
forces in the immediate neighbourhood of the descending slab. Modelling of the
expected stress distribution in a plate based on the varying force contributions from
its margins, taking into account age-dependence, can provide quite good agreement
with the estimated stress patterns with the continent (see, e.g, the work of Reynolds,
Coblenz & Hillis 2003, for the Indo-Australian Plate).

12.4 Measures of stress and strain
12.4.1 Stress measurements
A compilation of stress measurements from around the world is maintained
by the World Stress Map project of the Heidelberg Academy of Science and
Humanities. Nearly 16,000 reliable observations are included in their 2005 data
release (http://www-wsm.physik.uni-karlsruhe.de). The data base is built up
(i) earthquake focal mechanisms,
(ii) well-bore breakouts and drilling-induced features,
(iii) in situ stress measurements,
(iv) recent geological data such as fault slip analysis or alignments of volcanic
The primary indications of stress conditions in the lithosphere come from the
source mechanisms of well-located earthquakes. For larger events the global
centroid moment tensor catalogue (Section 11.5) provides an extensive suite of
information for the period since 1976. The ambiguity between the faulting
plane and the perpendicular auxiliary plane can only be resolved with additional
information such as, e.g., geological constraints. The interpretation of the focal
mechanism in terms of the principal stress directions requires the use of criteria for
brittle failure, such as the Anderson criteria introduced in Section 10.2.3.
12.4 Measures of stress and strain 275

66 77 88 99 110

40 40

32 32

24 24

Focal Mechanism
16 16
all depths

66 77 88 99 110

Figure 12.9. Stress conditions inferred from seismic focal mechanisms in central Asia.
[Courtesy of O. Heibach.]

In areas of distributed deformation such as the region from the Himalayan
front into Central Asia, the fault plane mechanisms for shallow events provide
a remarkable view into the changes in stress regime. There is thrust faulting
at the convergent margin that rapidly changes in the high Tibetan plateau to
normal faulting associated with east“west rifting with also strike-slip faulting. The
strike-slip events in particular show a dramatic rotation in direction around the
syntaxis in northern Burma, indicating a ¬‚ow of crustal material.
In regions with few seismic events, reliance has to be placed on information from
smaller events, and then it is important to have as good an azimuthal distribution of
stations as can be achieved to provide reliable results on source mechanism. As far
as possible determinations should be based on the main event in a sequence since
the mechanisms of aftershocks depend on the residual stress distribution rather than
the ambient stress ¬eld, and there can be signi¬cant variations between nearby
276 Lithospheric Deformation

Borehole breakout observations provide information on the orientation of the
stress ¬eld from the distortion of the hole. An initially circular borehole will tend
to be elongated in the direction of minimum horizontal stress by wall failure.
In-situ measurements are largely based on the use of the concept of stress relief
in which a volume of material is partially isolated from the surrounding rock and
the consequent change of shape is measured to determine the strain. The strain data
in conjunction with estimates of the elastic moduli can then be used to provide a
measure of the stress on the volume (see, e.g., Jaeger & Cook, 1979). A variety
of methods are used to provide the stress relief. In over-coring the measurement is
made of the deformation of a rock segment at the end of a borehole that is isolated
by coring the ¬nal section with an annular drill bit. Alternatively, a slit can be cut
in a borehole wall.
Hydraulic fracture methods rely on the creation and maintenance of fractures in
the borehole wall. Fluid is pumped into a sealed section of the borehole until the
wall fractures. The complete stress tensor can be determined from the pressure
information for the case when fracturing ¬rst occurs and the pressure needed to
hold the fractures open after failure. Unfortunately most deep observations are in
sedimentary basins, rather than the basement rocks of greater geophysical interest.

12.4.2 Strain measurements
Improvements in geodetic techniques, particularly though the use of space-based
methods, have revitalised studies of deformation at sub-seismic frequencies. With
a number of years of recording the use of very long baseline interferometry
(VLBI) has produced estimates of motion between the plates which are consistent
with geological estimates, but continue to improve in accuracy as further data is
On smaller spatial scales GNSS observations based on the analysis of signals
from constellations of satellites designed to aid navigation (e.g., GPS, GLONASS,
GALILEO) have become an increasingly important tool in the measurement of
earth deformation. With a minimum of four satellites observed simultaneously
the con¬guration of the station relative to the satellites can be determined with a
knowledge of the satellite orbits. The orbit information transmitted by the satellites
is suf¬cient for metre level accuracy, but precision observations at the centimetre
level or better require the orbital parameters to be improved using observations
from permanent ground stations and a set of corrections for other factors applied,
such as the presence of water vapour that can change the apparent distance to a
satellite. The GNSS observations can pick up coseismic signals, but are also well
suited to capturing deformation occurring at time scales greater than a month.
Such longer time scale phenomena are, e.g., the deformation associated with
viscoelastic relaxation following an earthquake, and also decadal estimates of strain
accumulation and plate motion, and their spatial variations.
GNSS observations indicate that the zone of deformation associated with the
12.4 Measures of stress and strain 277

boundary between the North American and Paci¬c plates, with its main expression
in the San Andreas fault system, is several hundred kilometres wide. This
helps to explain the discrepancy between the deformation estimated from the
cumulative slip of earthquakes and the predictions based on knowledge of the
velocities between the plates. Increasing density of geodetic information means
that seismological and geodetic information can be brought together to provide
more detailed information on the deformation processes whose most obvious
symptoms are earthquakes. There are indications that a number of different time
scales of post-seismic deformation may be present that could, e.g., represent rapid
afterslip on a time scale of a few weeks accompanied by slower viscoelastic
relaxation of mantle and lower crust layers.
In areas of active deformation direct observations of strain accumulation can
be obtained using appropriate instrumentation. The Plate Boundary Observatory
component of the EarthScope project in the western U.S.A. is using borehole
strainmeters that provide triaxial information and laser strainmeters for direct
measurements of strain in particular directions. A long-baseline laser strainmeter
measures the relative displacement between two end-piers separated by several
hundred metres. Laser light is directed through a straight evacuated pipe to a
re¬‚ector mounted on one pier that returns the light to a laser interferometer on
the other pier. Changes in optical path length are detected by the interferometer
in terms of the wavelength of light from a stabilized laser. The pipe is evacuated
to suppress the variations in the index of refraction of the air that would otherwise
affect the optical path-length measurement.
With high densities of permanent geodetic sites in many active source regions,
direct constraints are available on the mechanisms of earthquakes from the
deformation pattern. Models of the faulting process used in the analysis are tending
to become more sophisticated with multiple fault segments and realistic Earth
Multiple reoccupation of sites with short-term recording every few years
provides a means of extending deformation information with a limited number
of instruments. As an example, results from GPS studies across the complex
micro-plates of the Papua New Guinea Region have shown that the present-day
rates of change in relative position between stations are very consistent with the
average rates over the past few million years (Tregoning et al., 1998). In such an
active seismic region the patterns of observation can be disrupted on occasion by
the occurrence of a major seismic event between the reoccupation of sites.
The interferometric radar technique (InSAR) introduced in Section 11.5 provides
the opportunity for medium-term monitoring of deformation by the examination
of differences between repeat satellite passes across the same region. Such
information will be of particular value in detecting slow systematic vertical
movements, which are less easily detected using GNSS methods.
Holt et al. (2000) present an approach to extracting horizontal strain and rotation
278 Lithospheric Deformation

components through the determination of a smooth vector rotation function W
parameterising the horizontal velocity ¬eld

v(x) = re W(x) — x. (12.4.1)

Both the rate of strain elements and the rotation rate can be represented in terms
of the derivatives of W. Attention can be focussed on the horizontal components,
because even for rapid uplift (10 mm/yr), the contribution from vertical motion is
around 1%. An inverse problem can be posed with an objective function composed
of a sum of two parts, matching both direct observations of the rate of strain
tensor from geological and seismological information and velocity vectors from
geodetic measurements. The rate of strain ¬eld is required to satisfy St. Venant™s
compatibility relationships in the second derivatives of strain rate. Holt et al. (2000)
provide a number of applications to, e.g., western North America, Central Asia and
the Aegean.
An alternative approach for direct exploitation of geodetic information has been
developed by Spakman & Nyst (2002) using the combination of many observations
of differential velocities between pairs of stations. With observation of the
horizontal velocity vector v1 at site 1 and v2 at site 2, using the same reference
frame, the differential velocity
v1 ’ v2 = dx · + f12, (12.4.2)

where P12 is a suitable path between the sites and f12 is any contribution from faults
between the sites. In many regions it is a plausible assumption that faults are locked
at the surface, so that f12 = 0 . Then (12.4.2) forms the basis of a tomographic
inversion for the horizontal components of the velocity gradient tensor, and hence
the rate of strain and rate of rotation.
In Europe over 200 GPS observation sites provide more than 100,000
independent pairs of differential velocity, and the paths P12 can be chosen to
optimise areal coverage. The resulting deformation tomography can be cast in
terms of presentations of the principal components of rate of strain and rotation
patterns as in Figure 12.10. The resolution is very good and the standard deviations
in the model are generally much less than 5% of the imaged amplitudes. The results
in Figure 12.10 are obtained by assuming that all faults in Europe are locked at
the surface, and predominantly show the variation of the velocity gradient ¬eld
across Europe. Notable features at this scale are the extension in Fennoscandia,
the contraction in the Western Mediterranean and the distributed deformation in
the Italian“Aegean“Anatolian“Zagros plate convergence zones. A rotation rate of
2.5—10’8 /yr corresponds to 1.4—¦ /Myr; in the Aegean“Anatolian region the rotation
rates are beyond the limits of the contouring scale and attain values of the order of
1.0—10’7 /yr (5.7—¦ /Myr). Much more detail and systematic variations appear when
attention is concentrated on speci¬c regions.
12.4 Measures of stress and strain 279

Figure 12.10. Velocity rates for Europe determined by deformation tomography: (a) Prin-
cipal components of rate of strain, open vectors indicate extension, ¬lled vectors com-
pression; (b) rotation rates, dark shades represent counter-clockwise motion and patterning
indicates clockwise motion. [Courtesy of G. Tanasescu, W. Spakman, B.A.C. Ambrosius.]
280 Lithospheric Deformation

12.5 Glacial rebound
A convenient model for understanding the deformation of the Earth associated
with the loading and unloading of ice accumulation during the recent glacial
epoch is to use a model comprising a sequence of linear viscoelastic layers. The
correspondence principle (Section 5.5) can then be used to take results from perfect
elasticity and apply them to Laplace transformed variables.
As we have seen in the simple illustration in Section 7.4, the rapid removal of an
ice load allows the lithosphere to progressively return to its original con¬guration,
with maximum uplift in the zone which was originally most depressed. The uplift
of the Fennoscandian shield in the Baltic following the end of the last ice age has
produced raised beaches with current uplift rates reaching 9 mm/yr at the upper
end of the Gulf of Bothnia. The surface of the Earth deforms locally in response to
the changing ice“water distribution and signi¬cant effects can occur outside of the
zone directly affected by ice-loading.
One of the dif¬culties in exploiting observations of past sea level is that the ice
history is only imperfectly known and so aspects of the ice load have to be inferred
at the same time as information on the structure supporting the load. There can also
be signi¬cant water input from the melting of ice sheets at substantial distance, e.g.,
input from the Antarctic has an in¬‚uence in the northern hemisphere. The timing
of ice melting differs between locations and adds further complexity. The broad
distribution of the loading of major ice sheets makes a global approach using a
spherical harmonic representation very effective.
In Section 11.3 we have established the equations for elastic deformation
and shown how a system of six coupled differential equations in radius can be
established for the deformation of a self-gravitating Earth (11.3.57) in terms of the
coef¬cients of spherical harmonics. For the current long-term load problems we
can set the frequency ω = 0. The boundary conditions for surface loading are that
the radial stress σrr should balance the normal surface load, while the tangential
stress components σrθ and σrφ should vanish. In addition the load enters into
a continuity condition for gravitational potential gradient (see, e.g., Lambeck &
Johnson, 1998).
The surface load L can be expanded in spherical harmonics as
L(θ, φ) = LlmYl (θ, φ) = LlYl(θ, φ) (12.5.1)
l m=’l l

with associated gravitational potential •L(r, θ, φ) in r > re

re l
• (r, θ, φ) = Yl(θ, φ),
1 re l
= 4πre G LlYl(θ, φ). (12.5.2)
2l + 1 r
12.5 Glacial rebound 281

The new potential after deformation is

• = •L + •d, kl(r)•L(r),
where •d(r) = (12.5.3)

and kl is the Love number for gravitational potential.
The displacement vector u(r, θ, φ) is de¬ned in terms of a spherical harmonic
expansion by Love numbers h, as
•L(r) re l
u(r, θ, φ) = [hl(r)Yl(θ, φ)^r + r l(r)∇Yl (θ, φ)]
e (12.5.4)
g(r r
with ∇Yl the horizontal gradient of the spherical harmonic function, and g(r) the
radial distribution of gravity. The displacement Love numbers h, de¬ne the radial
and tangential deformation. The Love numbers are related to the variables U, V
and ¦ in (11.3.57) by
⎛ ⎞ ⎛ ⎞
hl(r) g(r)Ul(r)
⎝ l(r) ⎠ = 2l + 1 ⎝ g(r)Vl(r) ⎠ . (12.5.5)
1 + kl(r) ¦l
For linear viscoelastic strati¬cation, a Laplace transform as in (5.5.14), (5.5.19)
applied to the differential equations and boundary conditions reduces the problem
to an equivalent elastic problem, where the equations depend on the Laplace
transform variable p. In such a way a time-varying load can be introduced through
its Laplace transform and the deformation history described through the inversion
of the transformed results for displacement. Some care has to be taken in the
numerical inversion of the Laplace transform.
Realistic solutions of the viscoelastic equation systems require a high-resolution
of the coupled ocean“ice surface load. As ice melts load is transferred to the ocean
through increasing sea level; the ocean surface must remain an equipotential and
mass conservation is required between ice and ocean. The total change in the load
from oceans and ice can be expressed as
L(θ, φ, t) = ρI ”I(θ, φ, t) + ρW ”W(θ, φ, t). (12.5.6)
Here ”I is the change in effective ice height and ρI its density, ”W is the change in
water depth of density ρW. For an iceshelf grounded below sea level a correction
has to be made for the displacement of water by part of the ice load. The change in
water level has to be con¬ned to the area of the oceans and since coastline geometry
is a function of water depth the spatial constraint varies in time. It is useful to
introduce an ocean function O(θ, φ, t) equal to unity in the ocean and zero on
land. Then the change in water depth at time t relative to the initial conditions at
t0 can be expressed as
‚ζ(θ, φ, s)
”W(θ, φ, t) = ds O(θ, φ, s) , (12.5.7)
282 Lithospheric Deformation

where ζ(θ, φ, t) is the relative sea level function.
To use the Love number formalism outlined above, the surface loads need to be
expanded in terms of spherical harmonics as well as functions of time. Thus
L(θ, φ, t) = [ρI”Ilm + ρW”Wlm], (12.5.8)
l m=’l

where ”Ilm and ”Wlm are the spherical harmonic expansion coef¬cients. To
ensure an accurate representation of the time varying coastlines the expansion needs
to be carried to relatively high angular orders.
The primary information on rebound comes from estimates of past changes in
sea level inferred from the heights or depths of shoreline features at a known time
measured relative to present ocean surface (at time tp). Although a wide range
of observations can be exploited such as wave cut notches above current sea level
and submerged wavecut platforms, such information is not instantaneous. Other
proxies can only provide upper or lower bounds and good quality dating is essential
for good constraints on behaviour.
Sea level moves up and down relative to the crust, which is itself deforming due
to the changes in the loading from the combination of ice and water. The relative
change in sea level ζ from the present time can be represented as
”ζ(θ, φ, t) = ζ(θ, φ, t) ’ ζ(θ, φ, tp) (12.5.9)
The sea level change in (12.5.8) includes eustatic change ”ζe and the contribution
from the deformation of the crust and changes in the geoid as mass is displaced
”ζd. The eustatic component is the mean change in relative sea level over the
entirety of the oceans, and from mass conservation provides a measure of the
volume of ice on the continents,
1 dW00 ρI ”I00(t)
”ζ (θ, φ, t) = ≈’ ,
ds (12.5.10)
O00(s) ds ρW O00(t)

since the zero-degree harmonic corresponds to an average value. The deformational
components can be expressed in terms of the Love numbers as
•d(re )
”ζ (θ, φ, t) = ’ ur(re ) (12.5.11)

4πGre 1
= 1 + kl(re , t) ’ hl(re , t) — Ll(t)Yl(θ, φ),
g(re ) 2l + 1

where — denotes a convolution. To ensure mass conservation, a small contribution
dζ0 that is spatially uniform but time-dependent has to be added,
∞ l
”ζD Olm.
dζ0 ≈ ’ (12.5.12)
l=1 m=’l
12.6 Extension and convergence 283

The total relative sea-level change is then given by
”ζ(θ, φ, t) = ”ζe(θ, φ, t) + ”ζd(θ, φ, t) + dζ0(θ, φ, t). (12.5.13)
Lambeck & Johnson (1998) provide a number of examples of sea level behaviour
and describe the way in which inferences can be made about mantle structure even
though the ice history has major uncertainties. They demonstrate that there can be
considerable trade-off between different aspects of the response of the Earth, e.g.,
between upper mantle viscosity and apparent lithospheric thickness. Nevertheless
some very robust conclusions can be drawn, such as the requirement that the
average lower mantle has a higher viscosity by at least a factor of 20 than that
of the upper mantle (∼ 5—1020 Pa s). There is lateral variability in the viscosity of
the upper mantle and the apparent thickness of the lithosphere on broad horizontal
scales similar to those seen in seismic wavespeed and attenuation studies, but small
wavelength features cannot be resolved.
In most studies of glacial rebound a simple Maxwell body (Figure 1.3) has been
used for the viscoelastic constitutive equation. This is undoubtedly too simple a
rheology to provide a full description of the deformation, but can provide a good ¬t
to the limited available sea-level data. The estimates of viscosity are therefore best
regarded as an “effective” viscosity subject to the assumption of linearity.
Recent studies in and around Antarctica comparing analysis of sea-level data
and contemporary measures of uplift from direct space geodetic methods indicate
differences in the viscosity estimates for the millennia and year time scales. Such
results suggest that there may be a need to adopt a more complex rheology with
a non-linear dependence of strain-rate on stress, that yields a stress-dependent
effective viscosity. The ¬‚ow in the mantle needed to sustain the rebound has
to interact with the ¬‚ow patterns associated with mantle convection. Since the
rebound ¬‚ow rates are small the two effects can reasonably be superimposed (even
for a non-linear rheology).

12.6 Extension and convergence
The con¬guration of the lithosphere is modi¬ed through processes of extension
or compression, with the most obvious manifestations in the crust through the
in¬‚uence on surface geology. Failure in tension is signi¬cantly easier than
in compression (cf. Section 10.2.3) so that extension tends to impact on the
whole lithosphere in a localised region. Convergent processes tend to drive
mountain building with shallow zones of failure, but under certain circumstances
gravitational instability can develop in the lithospheric root.

12.6.1 Extension
Much work on tensional processes has concentrated on the development of rifts as
precursors to continental break-up. Many rifts appear to be linked to ancient zones
284 Lithospheric Deformation



Figure 12.11. Lithospheric detachment models: (a) lithospheric wedge; (b) Delamination
with ramps and ¬‚ats

of contrasting properties, such as major lithospheric lineaments or the edge of the
thick lithosphere beneath a craton as in the Tanzanian segments of the East African
Rifting can be accomplished by mechanical weakening of the lithosphere by
stretching, the intrusion of heat and interaction with the dynamical asthenosphere
beneath. Mechanical stretching accommodates the strain by the creation of faults
with large offsets in brittle layers and localised ductile deformation in weaker
regions. Depending on the speci¬c rheology of the lithosphere, strain may be
broadly distributed, or focussed on one large detachment fault. Rifting is frequently
associated with volcanism, and the intrusion of magmatic material, either passively
or actively, can be expected to enhance reduction of strength and localise strain.
Extensional strain can be taken up by magma injection into the lower lithosphere,
feeding dykes that penetrate the stronger lithospheric layers above. Small offset
faults will be formed above such zones of magma injection. Buoyant melt can
help to enhance stretching processes, whereas friction on fault surfaces will tend to
impede the thinning of the lithosphere.
Lister, Etheridge & Symonds (1991) have pointed out major differences in
continental margins that would have originally been contiguous before rifting
and the formation of the intervening ocean. The asymmetry of the margins can
be associated with their con¬guration relative to the detachment surface through
the lithosphere. If strength is concentrated near the surface, the situation before
break-up has a simple detachment surface with the development of a lithospheric
wedge structure (Figure 12.11a). However, multiple zones of strength can lead
to delamination with a set of ramps and ¬‚ats (Figure 12.11b). The subsequent
history will be signi¬cantly different for the upper and lower segments relative to
the detachment.
12.6 Extension and convergence 285

The mathematical description of lithospheric deformation is challenging because
of the strong temperature dependence of ductile rheologies, and the dif¬culty
of providing an adequate representation of the behaviour of brittle material. To
provide a continuum treatment, the brittle parts of the lithosphere can be assumed
to have a dense distribution of fractures on length scales smaller than the size of
lithospheric units so that a plastic rheology can be employed (see, e.g., Christensen,
1992). Such an approach cannot capture deformation on speci¬c faults or in narrow
shear zones, and so can only provide a generalised picture of structural evolution.
A simpli¬ed approach in which mechanical and thermal effects are decoupled
was introduced by McKenzie (1978) in a study of the evolution of sedimentary
basins. He considered the instantaneous stretching of the entire lithosphere
followed by thermal recovery. Subsequent work used a similar assumption
of pure shear with a uniformly strained lithosphere (with zero shear strain)
represented through a thin viscous plate. The non-linear effective viscosity
represents an average over the whole depth interval of the lithosphere. Such
models provide valuable insights into coupled thermal and mechanical effects in
deformation. For example, as the lithosphere thins in tension, conductive cooling
will increase the strength and thereby inhibit further extension (Houseman &
England, 1986). The changing con¬guration of the crust“mantle boundary and the
lithosphere“asthenosphere boundary during deformation bring gravitational forces
into play that themselves in¬‚uence the extensional process.
As we have seen in Section 12.2.4, the strength pro¬le of the continental
lithosphere is strongly dependent on water content, and can, in certain
circumstances, resemble a laminate of strong and weak components allowing
differential deformation between the different zones. For such a con¬guration the
use of a thin plate with an averaged viscosity will be somewhat misleading.
The treatment of the full complex rheology in deformation requires numerical
methods that can accommodate substantial distortion of the original material
con¬guration. This led Christensen (1992) to exploit an Eulerian formulation of
the coupled thermo-mechanical deformation in terms of ¬nite elements. In such a
treatment the material ¬‚ows through a ¬xed spatial grid and large deformation can
be readily followed, which is very useful in following the later stages of rifting.
Comparatively ¬ne sampling may be needed to provide a good representation
of material discontinuities in the Eulerian approach. If a Lagrangian approach
is used the grid is ¬xed to the material and as deformation proceeds the mesh
will be distorted, and remeshing is needed to keep the grid in a reasonable
con¬guration. Braun & Sambridge (1994) have demonstrated how such remeshing
can be accomplished dynamically to allow for the in¬‚uence of faulting.
Whatever representation is used for the numerical framework, there is a need to
suppress acoustic waves whose time scales are much shorter than the lithospheric
deformation. This can be done with either the assumption of incompressibility, or
the milder anelastic approximation ‚ρ/‚t = 0. In a two-dimensional problem the
286 Lithospheric Deformation

representation of ¬‚ow can be cast in terms of a stream function ψ (as in Chapter
7), but for three-dimensional simulations a full treatment in terms of velocities
is normally used. Density changes induced by extension or thermal effects are
unlikely to exceed 10%, and so the Boussinesq approximation (Section 7.5.2)
in which variations of density are included only in the gravitational forces can
generally be adequate. In the lithospheric case, because of the density contrasts
between crust and mantle, density changes associated with, e.g., the deformation
of the crust“mantle boundary must be included along with thermal effects.
A suitable rheological representation for ductile ¬‚ow is
™ ij = A exp ’ (E— + pLV — + μ— C)/RT σn’1σij, (12.6.1)
relating strain rate ™ to the stress-tensor σ. In (12.6.1), E— is the activation energy,
V — is the activation volume, μ— is the chemical activation potential associated
with changes in composition C, and pL is the lithostatic pressure. The power-law
dependence is introduced through the dependence on the second invariant of the
stress tensor σ = (σijσij)1/2 . The use of a chemical potential is a simpli¬ed
description of the variation of creep strength with composition. The next level of
complication would be to allow A and the index n to be functions of composition
C and so allow different creep behaviour in crust and mantle. The constant A in the
expression for the effective viscosity ·d can be eliminated by specifying a reference
viscosity ·ref for a reference temperature, pressure, and strain rate ™ ref :
·d = a exp ’ (E— + pLV — + μ— C)/nRT ™ (n’1)/n, (12.6.2)
a = ·ref exp ’ (E— + pref V — + μ— C)/nRT ™ ref . (12.6.3)
In the cooler, shallow parts of the crust we have brittle deformation that can be
approximated by a power law with a large stress exponent
™ ij = B(pL)σq’1σij, (12.6.4)
where, e.g., Christensen (1992) used q = 30. A brittle-plastic effective viscosity
·p can be extracted from (12.6.4), but again it is convenient to work in terms of
reference values and so de¬ne the yield stress invariant σy at which the plastic
strain equals ™ ref . The yield stress depends on lithostatic pressure as
σy = σy(0) + fpL, (12.6.5)
where the factor f represents the scaling between lithostatic pressure and pore-¬‚uid
pressure (cf. Section 10.2.3). The effective viscosity for the brittle (plastic)
deformation is then
·p = (σy ™ ref ) ™ ’(q’1)/q. (12.6.6)
The mechanical deformation is coupled to temperature through the advection of
12.6 Extension and convergence 287






Figure 12.12. Illustration of numerical modelling of lithospheric deformation in three di-
mensions through a three-dimensional perspective view of a surface at the base of the man-
tle lithosphere, with extension in one region leading to gravitational instability of dense
viscous mantle lithosphere just outside the extending zone. This model is used to explain
the occurrence of deep continental seismicity beneath the SE Carpathian mountains. [Cour-
tesy of P. Lorinczi and G. Houseman.]

heat and through heat generation processes. The temperature evolution equation
can be cast in the form, cf. (7.1.15),
ρCp + v · ∇T = ∇(k∇T ) ’ h(C) ’ [σ · ™ ’ ρ1g±th vzT ] (12.6.7)
where cp is the heat capacity, k the thermal conductivity, h(C) represents the rate of
volumetric heat production and ρ1 is the original density pro¬le. The two terms in
square brackets represent the contributions from frictional heating and the effect
of adiabatic compression. For consistency these two contributions must either
be neglected or included together; commonly if the Boussinesq approximation is
adopted they are neglected.
The dif¬culty with modelling extension processes is the choice of lithospheric
parameters, since the coupled mechanical“thermal behaviour can be quite sensitive
to the particular choice of parameters and approximations made (Christensen,
1992). The speci¬cation of boundary conditions plays an important role. There
is a considerable difference between the situation where a side boundary moves
with constant velocity or where it is ¬xed. Over a broad domain, extension in one
region will have the effect of throwing the neighbouring zones into compression
with consequent deformation (Figure 12.12).
288 Lithospheric Deformation

In a simple stretching scenario, the initial response is a slight updoming of
the crust-mantle boundary accompanied by a more pronounced upward shift in
the boundary between lithosphere and asthenosphere. Once this step has been
accomplished further extension can be comparatively rapid and extended crustal
thinning can be accomplished in a period of 10 Myr or so. Such a time interval is
quite short compared with thermal time scales and so there is a certain degree of
decoupling between the mechanical and thermal response.
McKenzie (1978) introduced a very simple model to capture the basic features
of sedimentary basin formation, using an instantaneous stretching of the entire
lithosphere followed by thermal recovery. When a segment of lithosphere of
thickness hL is stretched by a factor β the new thickness is hL/β, and to achieve
isostatic compensation hot asthenospheric material will rise to take the place of
the former thicker lithosphere. The thinned lithosphere will now have a steeper
thermal gradient than before and will be out of thermal equilibrium. The stretched
lithosphere will cool and thicken until it ¬nally reaches its original thickness and
temperature pro¬le. The immediate effect is a drop of the surface elevation by an
hL(ρM ’ ρL) + hc(ρL ’ ρc) 1
S= 1’ , (12.6.8)
ρM ’ ρw β
where hc is the initial thickness of the crust with density ρc, hL is the initial
thickness of the lithosphere with mantle density ρL, ρw is the density of water
¬lling the subsidence, and the asthenosphere has density ρM. Further subsidence
occurs as the lithosphere cools and an approximate solution of the equations of
thermal evolution is
4hLρL±th TM β π
S(t) = 1 ’ exp(’t/„s) ,
sin (12.6.9)
π2(ρL ’ ρw) π β
where the temperature of the asthenosphere is TM and the relaxation time „s =
h2/π2κH . As the surface continues to subside there is room for sediment to be
deposited. Since the sediment is denser than water the weight of sediment will also
make a small extra contribution to the subsidence.

12.6.2 Convergence
As can be seen from Figure 12.13, the majority of great earthquakes (magnitude
> 8) occur at convergent plate margins. Many of these earthquakes, like the 1960
Chilean event (Mw 9.6) and the 2004 Sumatra“Andaman event (Mw 9.3) have a
large thrust component associated with slip on the interface between the subducting
and over-riding plates. These thrust events are generally accompanied by tsunamis,
since the sea-¬‚oor is raised locally along a long fault segment. A few of the events,
e.g., the 1977 Sumbawa, Indonesia (Mw 8.1) and 2007 Kurile Islands (Mw 8.2)
events represent normal faulting in the plate as it enters the trench.
12.6 Extension and convergence 289

Certain segments of the convergent margins have had many great events in the
last century (e.g., Japan“Kuriles“Alaska, South America). The absence of such
events elsewhere does not imply absence of risk, merely that the build up of stress
has not yet reached critical levels. The 2004 Sumatra“Andaman segment had not
had an earthquake of this magnitude for at least 500 years. The last major event on
the Cascadia subduction zone (NW USA and Canada) occurred in January 1700,
and the timing can be pinned down quite precisely from the arrival times of tsunami
waves at sites in Japan.
The oceanic plate entering the subduction zone brings with it inherited structure,
due to age variations or the presence of seamounts. The result is that a typical
length of relatively homogeneous material entering the trench is of the order of
500 km and so can support a thrust earthquake of magnitude Mw up to around
8.7. Even larger events require multiple segments to fail through successive stress
transfer. The presence of obstacles in the seismogenic zone such as a seamount
being subducted beneath Shikoku in western Japan affects the slip pattern so that
the Nankai earthquake of 1944 has most energy release to the east of the subducting
seamount, but still a component to the west.
Large earthquakes along the continent“continent convergence zone in the
Himalayas and in the Burma subduction segment (cf. Figure 12.9) are infrequent,
but the population density in the region has grown rapidly over the last century. So
now the potential for human disaster is much enhanced.
A further class of great earthquakes at plate boundaries is dominantly strike-slip
events in transitional zones, e.g., Azores“Gibraltar or the Macquarie Ridge south

1952 1959 1906 1946 1949
2006 1965
1957 1963
1941 1920
1935 1950 1911

2005 1932 1906
1965 1996 2000
1935 1971
1938 19062007
1966 1974
1977 1995
1906 1985

1942 1929

Figure 12.13. Great earthquakes since 1905 (magnitude > 8). Indicative rupture areas are
shown for the largest events (dates in bold). The plate boundaries are represented with the
same convention as in Figure 12.7. The majority of great earthquakes occur at convergent
290 Lithospheric Deformation

of New Zealand. Most events away from speci¬c plate boundaries, e.g. in central
Asia, can be linked to the regional compressional regime. Nevertheless, the tectonic
settings of a few great earthquakes such as the 1998 event off Antarctica remain
unclear. The 1994 Bolivian event occurred at a depth of 590 km and was felt in
Toronto, Canada as a result of the radiation pattern from the source. This event
appears to have broken through the whole thickness of the deep plate.
In the convergent zones, almost all the energy and most of the slip is taken up
by the largest events and their major aftershocks. In between such events the faults
can be treated as locked so that stress (and accompanying strain) accumulates in
the vicinity of the plate boundaries. Because of the variations in the geometry of
convergent plate margins and a signi¬cant oblique component in convergence for
many subduction zones, stress concentrations can build up away from the main
seismogenic zones. The resulting faulting can have very different character from
the thrusts at the margin, e.g., the strike-slip faulting of the 1891 earthquake in the
Neo valley, central Japan, illustrated in Figure 10.13.
Present-day faulting patterns can be strongly in¬‚uenced by the prior
con¬guration of deformation. An example occurs in the Zagros mountains of Iran
where convergence is taken up by thrusting on listric faults, whose dip shallows
with depth, that originated in an extensional regime. Listric fault systems are found
in other thrust environments, as e.g., for the Chi-Chi earthquake in Taiwan whose
structural setting is illustrated in Figure 10.19.
Modelling of convergent processes can be carried out with similar tools to those
used for studying extension, subject to the appropriate boundary conditions. Rather
than undertake a full dynamic model of subduction, a number of authors have used
a kinematic prescription of the force distribution at the base of the lithosphere to
drive the convergence, as in Figure 12.14. In this way a wide range of conditions
can be explored through the variation of a limited number of parameters.
Not only do convergent processes produce increased elevation in mountain belts
that can impinge on weather patterns, but also the differential erosion arising from
prevailing wind and rainfall patterns can have an impact on the pro¬le of surface
relief and the stress/strain patterns in the crust. Such effects can be introduced into
the modelling by an appropriate speci¬cation of the surface conditions.
Figure 12.15 illustrates the deformation produced in the convergence of the
Paci¬c and Australian plates in South Island, New Zealand from the work of Batt
& Braun (1999). The main Alpine fault has a strong strike-slip component, but
the convergence leads also to signi¬cant vertical uplift and mountains of up to
4000 m elevation. High grade metamorphic material originating from lower crustal
conditions is frequently found at the highest elevations, and this tendency can be
explained by the ¬‚ow pattern associated with the convergence.
Batt & Braun (1999) have used a rheological model in which the material
behaves as a non-linear Maxwell viscoelastic body until the stress exceeds a critical
value so that brittle failure occurs. The viscosity is stress-dependent and thermally
12.6 Extension and convergence 291

Material removed
by erosion
Retro-shear Pro-shear

Upper crust
Imposed velocity singularity
Low-strength decollement

Net mass balance
assumed to be maintained
by ˜subduction™ or
distributed deformation
below the decollement

Figure 12.14. Kinematic model of convergence: the crust sits on a relatively strong mantle
with ¬nite ¬‚exural strength; deformation of the crust is driven by a velocity discontinuity
at the base leading to an internal wedge structure whose geometry is determined by the
internal strength of the crust and the strength of the decollement zone. [Courtesy of G.

activated with an Arrhenius form similar to (12.6.1) above, so that temperature can
feed back to in¬‚uence the mechanical properties of the model. The plastic failure
criterion is of the form

J2D + 12 T0p = 0, (12.6.10)

where J2D is the second invariant of the deviatoric part of the stress tensor, T0
the tensile strength and p the pressure, incorporating both lithostatic and dynamic
effects due to deformation. A uniform rate of internal heating is employed, with
a constant temperature condition at the base of the model. This approximation
neglects the cooling effect of the subducting material and is appropriate for young
orogens (<10 Myr duration). The calculations were carried out in a Lagrangian
frame using dynamic remeshing (Braun & Sambridge, 1994).
During the initial stages of deformation, two conjugate shear zones form, linked
at the velocity discontinuity at the base of the model, with dips close to 45—¦ .
The system evolves towards an assymmetric distribution of strain as the main
deformation is accommodated on the retro-side where the subducted material
descends. Material is rapidly advected through the other shear zone so that there is
292 Lithospheric Deformation

Secondary pro-shear zone
Retro-shear zone
Primary pro-shear zone
(a) Strain rate

(b) Thermal structure

Figure 12.15. Deformation patterns in convergence for a felsic crust, with total erosion
and a convergence rate of 10 mm/yr: (a) contour plot of the logarithm of the second in-
variant of the deviatoric part of the strain-rate tensor after 10 Myr of tectonic deformation
(higher rates are darker), with superimposed particle velocities indicated by the arrows.
The approximate location of the brittle“ductile transition is indicated by the dashed line;
(b) corresponding thermal structure, the initial thermal pro¬le is a gradient from 0—¦ C to
600—¦ C over the model depth of 30 km. [Courtesy of G. Batt.]

much less deformation on this side (Figure 12.15). The ¬‚ow pattern brings material
to the surface near the retro-shear zone. The rise of hotter material has the effect of
decreasing viscosity and so enhancing the concentration in the neighbourhood of
the retro-zone. The brittle-ductile transition is thereby brought closer to the surface.
Strain softening in the natural environment is likely to be more localised than in the
In the circumstances where the mantle lithosphere is more dense than the
underlying asthenosphere convergence will tend to thicken the lithosphere, with the
possibility of gravitational instablity leading to drips into the lower viscosity zone
beneath (cf. Figure 12.12). This Rayleigh“Taylor instability has been analysed by
Houseman & Molnar (1997) to see how the stress dependence of viscosity might
affect such convective instability for cold dense mountain roots. They employ both
a perturbation analysis for the initial growth of disturbances and ¬nite element
analysis for the case of ¬nite deformation. For a power-law rheology ( ™ ∝ σn)
the speed w at which the bottom of the dense layer descends grows with time as

n’1 ”ρ g 1/n
w= C h (tb ’ t) , (12.6.11)
n B
12.6 Extension and convergence 293

where ”ρ is the density contrast between the dense layer and the underlying
half space, h is the thickness of the layer and B is a measure of resistance to
deformation. C ∼1 is an empirical constant that depends on n, wavelength and
the gradient of the density distribution in the dense layer. For n > 1 the speed w
is small initially but increases signi¬cantly as t approaches tb at which time the
thickened portion of the dense layer detaches. The time tb is related to the initial
perturbation in the thickness of the layer zo by
n n
B n
tb = . (12.6.12)
”ρg C (n ’ 1)
This model may help to explain why many mountain belts that are built by crustal
thickening subsequently collapse with normal faulting and horizontal extension,
since the removal of dense mantle lithosphere would provide potential energy
to drive extension. The tongue of cool material projecting into the hotter
asthenosphere will produce lateral temperature gradients that drive convection. To
the sides of the descending material the lithosphere will thin. Once the dense
material drops off, warmer material will replace the lithospheric root and uplift
occurs. The process would be more effective for wet than dry materials and could
remove between 50 and 75 per cent of the mantle lithospheric material over a period
of about 10 Myr. The stress-dependent viscosity leads to a delay in the convective
thinning that would not occur for a Newtonian ¬‚uid, and could explain the time gap
between thrusting and normal faulting in the evolution of a number of current and
former mountain belts.
The In¬‚uence of Rheology: Asthenosphere
to the Deep Mantle

The contrast in behaviour between the lithosphere with stress transmission and the
asthenosphere with ¬‚ow linked to the mantle convection system represents one of
the distinctions that can be made between different parts of the Earth through the
nature of intrinsic rheology. However, we have to use indirect information to gain
information on the change with depth.
A further major contrast in rheological properties in the near surface arises from
the process of subduction where former oceanic lithosphere penetrates through
the mantle transition zone to interact with the phase boundary at 660 km and the
signi¬cant increase of mantle viscosity in the same neighbourhood.
A major control on rheology in Earth™s mantle is likely to be the in¬‚uence of
temperature, but this will be accentuated by the effect of other variables such as
grain size, which can at best be inferred indirectly.
Little is known of the rheology of materials under the pressure and temperature
conditions prevailing at the base of the mantle, but the complexities of the D
layer suggest that contrasts in physical properties comparable to those seen at the
surface may well prevail. The D layer has distinct, but somewhat variable, seismic
anisotropy that suggests a range of different deformation conditions.

13.1 Lithosphere and asthenosphere
As we have shown in Chapter 12, the de¬nitions of lithosphere are time-dependent
and also related to the phenomenon being investigated. The common feature is a
transition at depth from a region in which stress is transmitted and heat transport
is largely by conduction, to a new zone the asthenosphere where horizontal ¬‚ow is
expected in a mixed thermal regime.
This concept arose initially from interpretations of the dispersion of seismic
surface waves in terms of a zone with reduced shear wavespeeds at depth. This
low-velocity zone was associated with increased seismic attenuation that was
ascribed to the presence of a small amount of partial melting of the mantle
material. The low-velocity zone for S waves is well established in the oceanic
environment (Figure 13.1), but is less clear beneath the older parts of continents.

13.1 Lithosphere and asthenosphere 295

As the concept of plate tectonics developed, the asthenosphere was invoked to
accommodate the ¬‚ows associated with the movements of the plates. In this sense,
the lithosphere“asthenosphere boundary would represent the base of the zone that
moves coherently in plate motion. This is close to the concept of the ˜tectosphere™
introduced by Jordan (1975), who suggested major differences between oceanic
regions and the cratonic cores of continents based on the time differences in the
passage of the seismic phase ScS re¬‚ected from the core“mantle boundary. Jordan
(1978) added other evidence and geochemical arguments to suggest relatively deep
roots to the continents that would have a signi¬cant in¬‚uence on the way that plates
containing continents might move.
The distinction that we wish to make between the lithosphere and the
asthenosphere is dominantly rheological, based on a distinction between the time
scales for deformation. Such information is not directly accessible, and so we
must resort to a variety of proxy information to infer the presence of a change in
rheology. Most information comes from seismology, since this provides a means
of investigating present-day structure. The existence of a transition at depth is
corroborated by electrical conductivity studies that indicate the presence of a more
conducting layer. The logarithmic depth sensitivity of electrical sounding does not
place close bounds on the location of this layer.

Figure 13.1. Averages of shear wavespeeds in the Paci¬c Ocean as a function of plate age,
show a clear low velocity zone below a lithospheric lid. The averages are derived from a
model obtained using surface wave tomography for the fundamental and ¬rst four higher
modes. [Courtesy of K. Priestley.]

13.1.1 Seismic imaging
The zone of high seismic wavespeeds in the upper part of the mantle is often
interpreted as an expression of the lithosphere (cf Figure 12.1). Much of this
information comes from surface wave tomography exploiting frequencies in the
range from 0.006 to 0.20 Hz, so that the attainable vertical resolution will be of the
296 The In¬‚uence of Rheology: Asthenosphere to the Deep Mantle

order of 25 km. On this scale, the base of the high-wavespeed zone (the seismic
˜lid™) frequently appears quite sharp in oceanic regions,
The situation is much less clear for the continents, where the base of the
high-wavespeed region is not generally sharp. A variety of de¬nitions for the
˜seismic lithosphere™ have been used by different authors, e.g., based on the point
of strongest negative vertical gradient in velocity anomalies or where there is a drop
in absolute wavespeeds. Studies with good path coverage suggest that high seismic
wavespeed structures are mostly con¬ned above 220 km for SV waves, with some
faster zones at depth beneath Archaean regions. There are however indications of
some level of polarisation anisotropy so that faster SH wavespeed structures may
extend to greater depth than those for SV (Gung et al., 2004).




D ˜ 50 km
30S C ˜ 100 km
B 150 - 200 km

A > 200 km


100 km depth

130E 140E 150E 160E

Figure 13.2. Schematic representation of the variation of lithospheric thickness in central
to eastern Australia inferred from seismological studies, progressing from Proterozoic cra-
ton in the west to the eastern seaboard with Neogene volcanism. [Courtesy of M. Heintz
and S. Fishwick.]

Although it may be dif¬cult to get an absolute depth for the
lithosphere“asthenosphere transition, relative depths may be extracted by
using a consistent style of analysis for a tomographic model. Such results suggest
that lithospheric thinning may be accomplished as a sequence of steps related to
the history of continental assembly rather than a smooth gradient. Figure 13.2
shows the inferred progression of lithospheric thickness across the eastern half of
the Australian continent derived from interpretation of surface wave tomography.
The thickest lithosphere is associated with the oldest, cratonic material (early-
to mid-Proterozoic 2.0-1.8 Ga). The sequence of steps to thinner lithosphere
in the east appears to be related to the assembly of later belts of structure onto
the cratonic core. The thinnest lithosphere in the east has experienced Neogene
13.1 Lithosphere and asthenosphere 297

volcanism, after sea¬‚oor spreading that opened the Tasman Sea ceased around 80
The distinct steps in lithospheric structure appear to have been resistant to
thermal erosion over very long periods of time, given their clear association
with different age bands at the surface. For resistance to thermal erosion to
have persisted over such long intervals there needs to be some geochemical
distinctiveness, most likely the presence of larger amounts of harzburgite in the
older lithosphere. This mineral is both slightly less dense than other minerals and
very refractory. The lower density allows the thicker lithosphere to ¬‚oat on the
mantle beneath and still achieve isostasy.
The lithosphere“asthenosphere boundary between steps is likely to be somewhat
undulating and local thinning of the lithosphere may provide points where
kimberlitic volcanism can break through to the surface through a narrow pipe. Such
specimens as we have of upper mantle material come from the xenoliths brought to
the surface in the rapid ascent of the kimberlitic material, e.g., by removal of the
wall rocks in the conduit. The dominant colour is undoubtedly pale green, but the
extent to which the grain size and other characteristics of the xenolith material are
representative of their original state is unclear.

13.1.2 Seismic attenuation
Some of the most direct seismic evidence for the existence of a signi¬cant change
in properties at depth comes from seismic attenuation through the frequency
content of seismograms. As shown by Gudmundsson et al. (1994) seismic wave
propagation in ancient cratonic lithosphere shows very little loss with ef¬cient
transmission of high frequency energy for 2000 km or more. The P and S
wave arrivals on seismograms show similar frequencies for epicentral distances
to 16“18—¦ , and then rather dramatically the high frequencies in S are eliminated for
seismograms from sources at greater distances to the same stations. The loss of the
high frequencies is a direct result of enhanced seismic attenuation, and the simple
visual observations are reinforced by analysis of the spectral ratios of the P and S
waves. A region of low Q and reduced shear wavespeed has to lie beneath the high
wavespeeds of the Australian cratonic lithosphere.
The geometry of earthquakes is particularly convenient for such studies in
Australia, but it is clear that signi¬cant shear attenuation has to be sited beneath
the lithosphere for ancient continental regions. There are horizontal variations
in attenuation as well, and the lithosphere beneath younger domains tends to be
more attenuative than under cratons, but still has much less loss than for the
asthenosphere below.
The low shear wavespeeds and increased attenuation in the asthenosphere have
often been attributed to the presence of partial melt. However, geochemical
evidence suggests that, even where melt is present, it is unlikely to be much more
than 0.1%, which is insuf¬cient to make a noticeable decrease in shear wavespeeds.
298 The In¬‚uence of Rheology: Asthenosphere to the Deep Mantle

2 2
12 12
50 35 35
80 80 S O
110 110
Depth [km]



2 Myr Y - Yilgarn
12 Myr N - N.Aus
350 35 Myr Nishimura & Fishwick
S - SE.Aus
80 Myr
Forsyth, 1989 et al., 2005
O - Ocean (35 Myr)
110 Myr
0 500 1000 1500 4.2 4.4 4.6 4.8 0 500 1000 1500 4.2 4.4 4.6 4.8 5
β [km/s] β [km/s]
T [°C] T [°C]

Figure 13.3. Fitting of the representation of seismic wavespeeds, including attenuation,
from Section 11.1.2 to seismic velocity pro¬les from surface wave analysis. The ocean
pro¬les are taken from Nishimura & Forsyth (1989) and the continental pro¬les from a
study of the Australian region by Fishwick et al. (2005). The observational results are
shown by symbols and the ¬tted velocities and associated temperature pro¬les are indicated
by solid curves. A increase in grain size in the lower part of the asthenosphere is included
in the ¬tting. [Courtesy of U. Faul & I. Jackson].

Laboratory measurements, such as those shown in Figure 11.2, indicate that there
is a substantial reduction of shear wavespeed and dramatically increased shear
attenuation for temperatures 100-200—¦ below the solidus. The seismic properties
of the asthenosphere can therefore be explained by thermal effects without needing
to invoke signi¬cant melt. Figure 13.3 shows how the dependence of seismic
shear wavespeeds on depth, derived from surface wave observations, can be well
approximated using extrapolations of the results of laboratory experiments with
suitable temperature pro¬les. The ¬ts to the seismic velocities are achieved with
the modi¬ed Burgers model discussed in Section 11.1.2, and include a systematic
increase in grain size in the lower part of the asthenosphere.
The microstructural interpretation of the transition from lithosphere to
asthenosphere is in terms of a change of deformation regime from dominant
dislocation creep to a situation where diffusion creep, possibly enhanced by water,
is the most important process. The dislocation creep in the lithosphere is associated
with a power-law rheology, such as the n = 3 stress dependence in (9.3.9); whereas
diffusion creep is close to linear in stress.
The seismological aspects of the transition from the lithospheric domain to the
asthenosphere can certainly be accommodated within the modi¬ed Burgers model
(11.1.8)“(11.1.11), by a shift in the importance of the different components induced
by the increase in temperature with depth. Samples of upper mantle material from
xenoliths suggest that there may be a reduction in grain size accompanying the
change in thermal regime, enhanced by the presence of any traces of water; there
might also be an accompanying change in the activation energy.
13.1 Lithosphere and asthenosphere 299

13.1.3 Seismic anisotropy
Further evidence for the in¬‚uence of the asthenosphere comes from changes in
the pattern of seismic anisotropy with depth. The orientation of minerals under
strain ¬elds (as in Section 10.1) leads to anisotropy in elastic moduli depending
on the level of consistency in the crystallographic orientations. Major mantle
minerals, such as olivine and clinopyroxene, show strong anisotropy in both P and
S waves. Lattice preferred orientation (LPO) of crystals could arise from the active
deformation of the asthenospheric mantle that accommodates the absolute plate
motion, or be imposed during past deformation and subsequently ˜frozen™ in the
lithosphere during post-tectonic thermal relaxation.
Since olivine is the most abundant and deformable mineral in the upper mantle,
the direction of fastest S wavespeed is frequently taken as an indicator of the a-axis
([100]) of olivine, the dominant slip direction (Ismail & Mainprice, 1998).
There are two classes of rather different information that provide constraints on
seismic anisotropy: (i) the azimuthal variations in the dispersion of surface waves
can be interpreted in terms of the vertical distribution of angular variations in S
wave speeds, and (ii) the patterns of splitting in the times of arrival of orthogonally
polarised shear body waves as a function of the azimuth to the source can be
interpreted in terms of models of seismic anisotropy.

Surface wave anisotropy
For surface waves, Rayleigh waves provides information on the angular variations
in the wavespeeds of vertically polarised shear waves (SV) and Love waves provide
information on the variations for horizontally polarised shear waves (SH). The local
pattern of variation of the phase speed of the dispersed wavetrains as a function of
azimuth ‘ depends on sin(2‘), cos(2‘) and sin(4‘), cos(4‘). For Rayleigh waves
the dominant variation is with 2‘, whereas for Love waves the most important terms
are those in 4‘ that need much better path coverage to be resolved. As a result,
most studies have concentrated on Rayleigh waves. Debayle, Kennett & Priestley
(2005) have used Rayleigh wave dispersion on paths from 1000 km to 10 000 km
length across the globe, measured from waveform matching techniques, to produce
global maps of SV wave anisotropy as a function of depth (see Figure 13.4). Under
most continents and the oceans the anisotropy decreases in amplitude between 100
km and 200 km depth, with little change in the directions of the fastest SV wave
propagation. However, an exception is provided by the Australian continent which
has the highest plate velocity of any continent of nearly 70 mm/yr. Under Australia
the direction of fastest SV wave propagation shifts from roughly east“west at 100
km depth to more nearly north“south at 200 km depth. The shallow regime is
consistent with anisotropy acquired in past deformation, whereas the orientation
at 200 km is close to the direction of absolute plate motion. This suggests that
shear on the fast-moving plate has in¬‚uenced mineral orientation within a zone that
300 The In¬‚uence of Rheology: Asthenosphere to the Deep Mantle

2% peak to peak anisotropy
100 km

240 300 0 60 120 180
200 km

240 300 0 60 120 180

Figure 13.4. SV-wave heterogeneity and azimuthal anisotropy (black bars oriented along
the axis of fast propagation) at 100 km and 200 km depth obtained from the inversion
of 100779 Rayleigh waveforms (Debayle et al., 2005). Hot spot locations are indicated
with circles. The positive S wavespeed anomalies are emphasised, and negative anomalies
are shown with patterns. The length of the black bars is proportional to the maximum
amplitude of azimuthal anisotropy. The reference velocity for SV wave variation is PREM
(Dziewonski & Anderson, 1981).

still displays the fast shear wavespeeds that would normally be associated with the
The alignment of the fastest shear wave propagation with the direction of plate
13.1 Lithosphere and asthenosphere 301

movement is normally regarded as a characteristic of the asthenosphere, but the
example above indicates that the rate of strain may be important.

Body wave anisotropy
In the presence of anisotropy, time differentials arise between the arrivals of
quasi-shear waves. The ˜fast™ and ˜slow™ polarisations are orthogonal to each other
and the polarisation directions are characteristic of the anisotropy of the medium
(Figure 13.5).

Ss Sf


Sf Ss


Figure 13.5. Passage of an S wave from an isotropic zone into an anisotropic region leads
to the generation of two distinct quasi-S contributions (Sf , Ss ) with orthogonal polarisation
but different wavespeeds. As a result there is a separation in time, δt, between the arrivals
of the two different polarisations.

For seismic body waves, measurements of the time separation between different
S wave polarisations, shear wave splitting measurements, are commonly performed
on shear-wave phases with a core refraction such as SK(K)S or PK(K)S. Such phases
are generated with P to S conversion at the core“mantle boundary (CMB) and are
thus polarized along the radial direction as they enter the mantle, i.e., along the
great circle between source and receiver. These phases arrive at a station with a
nearly vertical incidence, and energy appearing on the tangential component has
to be acquired in passage through the mantle. Thus a difference in the shear wave
times measured at the Earth™s surface represents the vertically integrated effect of
anisotropy from the CMB to the surface, with no indication of the depth location of
the anisotropy source. Because conversion to SV waves occurs on the receiver side
at the CMB, no contamination arises from structure near the source.
The measurement of shear wave splitting at a seismological station yields two
parameters: •, the orientation of the polarization plane of the faster S wave
and δt, the delay between the arrival times of the fast and slow split waves.
Petrophysical studies suggest that anisotropy will mostly be located in the upper
mantle, between the 410 km olivine“spinel phase transition and the Moho. Some
smaller contributions from the D layer, the lower mantle and the crust may also be
302 The In¬‚uence of Rheology: Asthenosphere to the Deep Mantle

present. The orientation of the polarization plane of the fast S-wave, • is taken
as a proxy for the orientation of the [100] axis of olivine in the upper mantle.
The delay time δt will be a function of the intrinsic anisotropy, the thickness of
the anisotropic layer, the orientation of the ray path with respect to the elastic
matrix of the anisotropic medium, and the vertical coherence of the mantle fabric.
Where measured orientations appear to be correlated with the super¬cial geological
structures, the anisotropy should be shallow with a signi¬cant crustal component.
Otherwise an origin for the shear wave splitting needs to be sought at greater depth.

13.1.4 Asthenospheric ¬‚ow
In many areas it has proved dif¬cult to reconcile the observations of shear wave
splitting for body waves and azimuthal anisotropy from surface waves. Attempts


. 11
( 16)