. 10
( 16)


11.4.3 Comparison with observations
We have considered seismic wave phenomena in a spherical Earth, both the free
oscillations and the travel times for higher-frequency seismic phases. Even though
the dominant variation in physical properties is with depth, there are signi¬cant
three-dimensional effects. So how far do the simpli¬ed models represent the
observed behaviour?
In Figure 11.19 we show the travel times for the full set of picked phases for
a set of 104 events (83 earthquakes, 21 explosions) that have a well controlled
hypocentre. The travel times for the rich set of later phase readings (57 655 phases)
are corrected to surface focus and are compared with the predictions from the
AK135 model.
244 Seismology and Earth Structure




Time [min]




” [°]

20. 40. 60. 80. 100. 120. 140. 160. 180

Figure 11.19. Travel times for the AK135 model for surface focus superimposed on the
times of phases for a set of test events.

There is a good correspondence between the observed and calculated times, but
we note some spread about the theoretical curves. The differences arise in part
because of the limitations in picking the actual onset of phases on seismograms,
but a signi¬cant component comes from the departure of the real Earth from
the spherically symmetric model AK135. The differences in timing contain
information about the three-dimensional structure encountered in passage between
the source and receiver that can be exploited to delineate such structure (see Section
The long-period disturbances propagating from the 2001 Mw 8.4 event in Peru
across the globe are shown in Figure 11.20 for both the radial (R) and transverse
(T) components of motion at the stations of the GEOSCOPE network. The radial
component is taken along the great-circle between the hypocentre and the station
and the transverse component is orthogonal to the radial in the horizontal plane.
11.4 Probing the Earth 245

120. .RER

Delta [deg]



1200. 1800. 2400. 3000.

120. .RER

Delta [deg]



1200. 1800. 2400. 3000.
Time [s]

Figure 11.20. Radial (R) and transverse (T) component seismograms for the Mw 8.4
shallow focus event in Peru in June 2001 at GEOSCOPE stations, with a low pass ¬lter
to 0.025 Hz. The main seismic phases are indicated.

The main seismic phases are indicated in Figure 11.20; even for this very large
event with a complex source waveform, there is a good correspondence to the
predicted phase times from the AK135 model. The large late surface waves on both
R and T components are the travelling wave equivalents of the free oscillations. The
fundamental mode Love waves, LQ, on the transverse component are equivalent to
the 0Tl modes and arrive somewhat earlier in the seismograms than the Rayleigh
waves, LR, on the radial component. The fundamental mode Rayleigh waves
correspond to the 0Sl modes, and are preceded by overtones with faster propagation
The equivalence of the travelling wave and normal mode representations for
seismograms is illustrated by the synthetic seismograms presented in Figure
11.21 for the same stations. These seismograms are calculated by normal mode
summation from (11.3.37) for the PREM model, using the source mechanism
derived from analysis of waveforms for the Peru event (particularly the surface
waves). A summation over spheroidal modes is used for the radial component
246 Seismology and Earth Structure

120. •RER

Delta [deg]



1200. 1800. 2400. 3000.

120. •RER

Delta [deg]


•SCZ *

1200. 1800. 2400. 3000.
Time [s]

Figure 11.21. Simulation of the Peru event shown in Figure 11.20 using the sum of the ¬rst
15 normal mode branches for the PREM model, with again a low pass ¬lter to 0.025 Hz.
[Courtesy of K. Yoshizawa].

and over toroidal modes for the transverse component, ¬fteen mode branches are
included which is suf¬cient to provide a synthesis of the S body waves.
The representation of the source of the Peru event used in Figure 11.21 is not
exact, since it is derived from analysis of the observations; the point moment tensor
model does not include the full effects of the source in this very large event and
undoubtedly oversimpli¬es the seismic radiation characteristics. Even though the
synthetic seismograms in Figure 11.21 have been calculated with this simple point
source representation, there is a good general correspondence with the observed
seismograms in Figure 11.20. However, the energy in the body waves is somewhat
underestimated compared with the surface waves, and the balance between the S
and SS energy is not quite correct. The spherical PREM model provides a good
general prediction of the character and timing of both the Rayleigh and Love
wave trains. We see from Figure 11.20 that there are noticeable differences in the
nature of the observed seismograms for stations at similar distances. This arises
because of the cumulative effect of three-dimensional structure on the passage of
the surface waves, so that long oceanic paths impose a very different dispersion than
an equivalent continental path. The differences in the long-period seismograms
11.4 Probing the Earth 247

between the observations and the modal predictions form the basis of methods for
determining the large-scale variations in three-dimensional Earth structure.

11.4.4 Imaging three-dimensional structure
As we have seen, the use of a model with purely radial variations in properties
provides an excellent representation of the major behaviour of seismic phases. The
presence of three-dimensional structure manifests itself in a variety of ways, which
have been exploited to provide images of such structure; Kennett (2002) provides
a detailed discussion of the styles of tomographic inversion and their contrasting
The splitting of the frequencies of the free oscillations of the Earth with respect to
the angular order m is an important source of information for the lowest frequency
modes that have a sensitivity to density and seismic wavespeed. At intermediate
frequencies the commonest procedure comes through the ¬tting of portions of
seismograms using a perturbation development based on the normal mode theory
outlined above. Surface waves can be described through the summation of simple
mode branch contributions, but body waves need multiple branch contributions
with coupling between the coef¬cients.
The arrivals of seismic phases are exploited in a variety of ways. High frequency
information can be derived from the compilations of readings from seismic stations
across the globe. Careful reprocessing of such catalogues, including relocation
of events and association of arrivals, provides a major source of information for
P and S waves and many later phases (see, e.g., Engdahl, Buland & van der
Hilst, 1998). At somewhat lower frequencies a substantial data set has been built
up for both absolute and differential times through the use of cross-correlation
methods, including the use of synthetic seismograms. The way in which such lower
frequency waves interact with structure is not fully described by ray theory and a
number of schemes have been developed to represent the zone of interaction around
the ray path.
Many recent studies use a wide range of different styles of information to try to
achieve the maximum level of sampling of the Earth™s mantle (see, e.g., Masters
et al., 2000). However, such studies face the complication of assigning relative
weighting to different classes of information and also of combining information in
different frequency bands.
Three-dimensional models of the variation in seismic wavespeeds are normally
displayed as deviations from a radially-strati¬ed reference model. Many
studies of seismic tomography have concentrated attention on a single wavetype
(particularly S). Current S images derived from long-period seismic data (such
as free oscillations, waveforms of multiple S phases) provide good de¬nition of
heterogeneity with horizontal scales larger than 1000 km across the whole of the
mantle. The model of Megnin & Romanowicz (2000) based on just the use of
SH waves (Figure 11.22) illustrates the capabilities of waveform inversion on a
248 Seismology and Earth Structure

global scale. This model was derived using coupled normal modes to represent
long-period body wave phenomena, with an expansion of the perturbations of
wavespeed structure in terms of spherical harmonics to order 24 and a spline
representation in radius. We see in Figure 11.22 the strong variations in shear
wavespeed near the surface and the core“mantle boundary compared with the more
modest variations in the mid-mantle.
Somewhat higher resolution of three-dimensional structure can be achieved
using arrival times extracted from seismograms for both P and S waves, but at
the cost of less coverage of mantle structure. A prerequisite for good quality
imaging in tomographic inversion is multi-directional sampling through any

Figure 11.22. Maps showing the variation of shear wavespeed at depths of 200 km, 1300
km and 2700 km for the SH model SAW24B16 of Megnin & Romanowicz (2000) dis-
played as deviations from the PREM reference model (Dziewonski & Anderson, 1981).
The ellipses are scaled in proportion to the radius of the section.
11.4 Probing the Earth 249

zone. The con¬guration of earthquake sources and, mostly, continental seismic
stations restricts the sampling for body waves, unless multiply re¬‚ected body
waves can be exploited. In those regions where different studies have achieved
a comparable coverage, the major features of the tomographic images are in
reasonable correspondence, even though the details vary. The highest levels of
heterogeneity are found near the Earth™s surface and near the core“mantle boundary
(see Figure 11.22). More subtle features appear in the mid-mantle, including
relatively narrow zones of elevated wavespeed that are likely to be produced by
past subduction.
In the uppermost mantle, the ancient cores of continents stand out with fast
wavespeeds, while the mid-ocean ridge system and orogenic belts show slow
wavespeeds compared with the reference model. Below 400 km depth the
high-wavespeed anomalies are mostly associated with subduction zones; in some
regions they extend to around 1100 km depth, but in a few cases tabular fast
wavespeed structures seem to extend to 2000 km or deeper. The base of the mantle
shows long-wavelength regions of higher wavespeeds, most likely associated with
past subduction, and two major regions of slow wavespeed beneath the central
Paci¬c and southern Africa (see Figure 11.22) which may represent sites of
upwelling of hotter material.
The early success of seismic tomography came from the striking images of
large-scale three-dimensional structure and, later, of the details of subduction
zones. The interpretation of such images is based ¬rmly on the variations in
seismic wavespeed. Thermal processes can be expected to play a major role, but
chemical heterogeneity could also be important particularly in the regions with
strong variability at the top and bottom of the mantle.
Results for a single wavespeed are not suf¬cient to indicate the nature of the
observed anomalies. Recent developments in seismic imaging are therefore moving
towards ways of extracting multiple images in which different aspects of the
physical system are isolated. This may be from P and S images (preferably from
common data sources) or via the use of the bulk modulus, shear modulus and
density. Such multiple images of mantle structure encourage an interpretation in
terms of processes and mineral physics parameters, since the relative variation of
the different parameters adds additional information to the spatial patterns.
In principle, a signi¬cant increase in understanding of heterogeneity can be
achieved if both P and S information are exploited. The P wavespeed ± depends
on both the bulk modulus κ and the shear modulus μ as

± = [(κ + 4 μ)/ρ]1/2, (11.4.1)

where ρ is the density. We can thus isolate the dependence on the shear modulus μ
and the bulk modulus κ, by working with the S wavespeed

β = [μ/ρ]1/2, (11.4.2)
250 Seismology and Earth Structure

and the bulk-sound speed φ derived from both the P wavespeed ± and the S wave
speed β
φ = [±2 ’ 4 β2]1/2 = [κ/ρ]1/2, (11.4.3)

which isolates the bulk modulus κ. This style of parameterisation has been
employed in a number of studies (see Masters et al., 2000, for a comparison).
An unfortunate complication in the combined use of P and S information comes
from the uneven geographic distribution of data. Whereas S wave data are available
from globe circling paths; P wave information, dominantly from travel times, is
dictated by the location of seismic sources and receivers. There is some component
of P wave information in long period waveform data but this does not provide
strong constraints on mantle structure. It therefore remains dif¬cult to compare full
global coverage from S with information derived from P waves with a much more
limited geographic coverage. The inclusion of later phases and differential times
helps but the P travel times still provide the dominant information.
An alternative approach is to restrict attention to paths for which both P and S
information is available; the consequence is that sampling of the mantle is reduced
but the reliability of the images is high and direct assessment can be made of
the relative behaviour of either P and S wavespeed variations, or the variations in
bulk-sound speed φ and shear wavespeed β (see, e.g., Kennett & Gorbatov, 2004).
In Figure 11.23 we show results from the model of Kennett & Gorbatov (2004)
which was constructed with 2—¦ —2—¦ cells and 18 layers through the mantle, using
a joint inversion of P and S arrival-time data with common source“receiver pairs,
with light damping and a broad residual range designed to capture strong features
in the uppermost and lowermost mantle. A linearised inversion is ¬rst performed
for P and S separately, and then a joint inversion for φ and β is undertaken with
three-dimensional ray tracing.
For the main shield regions, such as western Australia, there is a strong shear
wavespeed anomaly (up to 6% or more) down to about 250 km accompanied by
a somewhat weaker bulk-sound speed signature (Figure 11.23a). In contrast, the
major orogenic belts from southern Europe to Iran (and also western N. America)
show rather slow S wavespeeds, accompanied by fast bulk-sound speeds. A similar
behaviour is evident for the Red Sea and E. African rifts. The anti-correlation
of bulk-sound and shear wavespeed is pronounced, a large component can be
thermal because of the very strong reduction of the shear modulus as the solidus is
reached, but volatiles may also be signi¬cant. In eastern Asia we see some portions
of subducted slab, e.g., the Ryukyu arc, standing out very clearly by their high
shear wavespeeds from the lower background; these also have some expression in
bulk-sound speed.
Only in a few places are there clear indications of slab-like behaviour extending
below 1300 km depth, most notably in southern Asia and beneath the Americas.
These two features have been linked to subduction in the past 80 Ma at the
northern edge of the Tethys Ocean and of the now-extinct Farallon plate beneath
11.4 Probing the Earth 251

Figure 11.23. Variation of bulk-sound speed δφ/φ and shear wavespeed δβ/β relative to
the AK135 reference model. (a) layer from 100-200 km depth, (b) layer from 1000-1200
km depth, (c) layer from 2000-2200 km depth.
252 Seismology and Earth Structure

the Americas. In Figure 11.23(b) we see striking fast shear anomalies extending
from Iran to Indonesia that are almost absent in the bulk-sound image. Away from
the major features, the two wavespeeds show comparable levels of variability on
intermediate scales, with a weak anti-correlation in the patterns that is compatible
with minor thermal ¬‚uctuations. The relatively narrow, slab-like structures in
southern Asia and the Americas become less coherent with depth and appear in
places to link with drip-like features in the lowermost mantle. Cross-sections of
these structures can tend to be misleading because of the in¬‚uence of oblique
cuts. Indeed, it is dif¬cult to follow the behaviour in three dimensions because
of the various factors that can in¬‚uence the amplitudes of the imaged wavespeed
The character of the heterogeneity patterns in the mantle changes with depth,
notably below 2000 km. In the interval around 2100 km, illustrated in Figure
11.23(c), the bulk-sound speed variations are subdued. There is, however,
signi¬cant variation in S wavespeed with striking anomalies, especially in Asia,
with a very different pattern from that seen at 1100 km depth.
The character of the heterogeneity regime in the mantle undergoes further change
as the core“mantle boundary is approached. The amplitude of bulk-sound variation,
which is very low as we have seen near 2100 km, increases with depth towards
the core“mantle boundary with a pattern of variation that increasingly becomes
uncorrelated with shear wavespeed. The amplitude of the shear heterogeneity also
increases rapidly with depth. Just above the core“mantle boundary, in the D zone,
a wide range of different pieces of information indicates the presence of extensive
but variable heterogeneity including variable seismic anisotropy and narrow zones
with very low seismic wavespeeds. A striking feature of the bulk-sound speed
and shear wavespeed distributions at the base of the mantle is the discordance in
the patterns of variation (Masters et al., 2000). Such behaviour is not compatible
with a simple thermal origin and suggests the presence of widespread chemical
heterogeneity at the base of the mantle.

Interpretation of seismic heterogeneity
The interpretation of the results of seismic tomography depends on understanding
the controls on seismic wavespeeds under the conditions prevailing in the Earth™s
mantle. This requires a strong input of information from mineral physics, both
experimental methods and ab initio calculations, as discussed in Sections 9.4, 9.5.
It is likely that a substantial component of the seismic behaviour is controlled by
temperature, but the in¬‚uence of composition via major element chemistry should
not be overlooked. A further complication comes from the dif¬culty of extracting
absolute velocity information from tomography, the in¬‚uences of the damping
and regularisation used in the inversion tend to lead to an underestimation of the
amplitude of the anomalies, even though the spatial pattern may be appropriate.
In particular, we should note that the requirements for tomographic interpretation
are somewhat different than for equations of state. In Chapter 6 we have discussed
11.5 Earthquakes and faulting 253

the development of material properties under hydrostatic pressure. The Eulerian
(Birch“Murnaghan) formulation of ¬nite strain is generally preferred because of
the weak dependence on the pressure derivatives of the bulk modulus. This Eulerian
approach gets us to a suitable reference state, but the information we now need
on perturbations from that state is best described by an incremental Lagrangian
treatment for non-hydrostatic strains, for which anelastic effects can be signi¬cant.
We have already noted in Section 11.1.2 the rapid changes in shear modulus with
elevated temperature as the solidus is approached. The consequence is a strongly
non-linear dependence of the elastic moduli with temperature. To achieve the same
size of change in seismic wavespeed will require a larger temperature contrast for
fast anomalies than for slow anomalies where the temperature derivative is larger.
Because of the strong in¬‚uence of shear wavespeed anomalies on P wavespeed
images (e.g., Kennett & Gorbatov, 2004), the ratio of the P and S anomalies
R±β = δ ln ±/δ ln β (11.4.4)
has little diagnostic value for the in¬‚uence of composition, whereas the equivalent
ratio for bulk-sound speed
Rφβ = δ ln φ/δ ln β (11.4.5)
is more suitable for recognising the competing effects of temperature and

11.5 Earthquakes and faulting
The description of an earthquake source as a propagating rupture across a fault
plane captures much of the behaviour with respect to seismological observations.
In reality there is a dynamic interaction between the propagating rupture front and
the surrounding medium. Because the rheological properties vary with depth, the
physical behaviour in a shallow event can differ signi¬cantly between the near
surface and depth (Scholz, 1990).
For smaller earthquakes, the point moment tensor description introduced in
Section 11.2 provides a good representation of the behaviour, and the radiation
patterns for P and S waves can be used to infer the mechanism for the event. For
shallow events the P waves overlap with the pP and SP arrivals re¬‚ected from
above the source to give a complex interference packet, whose character depends
on source depth. Waveform modelling procedures can then be used for the onset
of the body waves to extract the mechanism, source depth and time history of the
event (see, e.g., Stein & Wysession, 2003, Chapter 4).
For all events with magnitude greater than 5.5, a systematic procedure is applied
to extract the centroid moment tensor by matching long-period seismograms
from across the globe with modal synthetics for frequencies below 0.0125 Hz,
with some allowance for major three-dimensional structure. The catalogue at
http://www.globalcmt.org/CMTsearch.html covers the period from 1976 onwards.
254 Seismology and Earth Structure

Figure 11.24. ENVISAT interferogram of the Bam earthquake in Iran, 2003 December 26,
Mw 6.6; the InSAR analysis was carried out using DORIS software. Each interferometric
fringe represents 28.1 mm of line-of-sight deformation. The fault is predominantly strike-
slip, with a peak slip of around 2 m. The elastic deformation cannot be explained by slip
on a single planar fault (Funning et al., 2005). The areas of low coherence in the centre of
the image are the heavily damaged Iranian cities of Bam and Bavarat. ENVISAT SAR data
were made available through the European Space Agency (ESA). [Courtesy of J. Dawson.]

The number of high-quality broad-band stations across the globe has signi¬cantly
increased in recent years, allowing more events to be characterised.
The seismological observations do not, by themselves, resolve the difference
between the fault plane and the orthogonal auxiliary plane for a dislocation source.
However, the patterns of ground deformation for shallow events provide additional
information. As a consequence there has been a fusion of results from geodetic
and seismic techniques to understand earthquake behaviour. Measurements of
post-seismic displacements using space-geodetic methods such as GPS provide
important constraints on possible fault behaviour. For very large earthquakes
such as the Mw 9.3 Sumatra“Andaman event on 2004 December 26, perceptible
deformation can be measured thousands of kilometres from the event itself.
In certain circumstances, a very direct view of ground deformation induced
by an earthquake can be extracted from interferograms from satellite based
synthetic-aperture radar (InSAR). The phase difference between images before and
after the earthquake can be mapped into interference fringes that map the patterns
of displacement, with particular sensitivity to vertical deformation. In particular for
11.5 Earthquakes and faulting 255

near-surface faulting in arid areas, the InSAR approach provides vivid images of
ground deformation at ¬ne resolution (∼10 m) in areas 10“100 km wide, that can be
matched to models of fault behaviour (Figure 11.24). For the 2003 earthquake that
devastated the ancient town of Bam in Iran, InSAR images revealed the main fault
trace that had been missed in the original post-event survey because of its limited
surface expression (Funning et al., 2005).

8 1.8
Depth [km]

Slip [m]




Figure 11.25. The slip distribution for the 1995 Kobe earthquake, Japan, proposed by
Yoshida et al. (1996) related to a map view with contours for sediment thickness.

A major earthquake is not well described by a simple rupture, the displacement
on the fault is spatially varying and episodic in time. The combination of near-¬eld
observations from strong ground motion instrumentation and teleseismic results
from distant seismic stations can be used to constrain models of rupture behaviour,
described as the superposition of the effects of many sub-faults. Characteristically,
there are patches where relatively high moment release occurs, often referred to
as asperities, and there may also be barriers where rupture is impeded or even
stopped. Figure 11.25 illustrates the slip distribution inferred by Yoshida et al.
(1996) for the 1995 Kobe earthquake from a wide range of information. There is
bilateral slip on two segments with differing dip, and multiple concentrations of
moment release along the extended fault trace.
Estimates of the slip distribution can now be made quite quickly following a
major event. In such rapid analysis it is usual to adopt a single fault plane. More
¬‚exibly parameterised models, including multiple fault segments, tend to evolve as
the full range of information with geodetic constraints becomes available.
256 Seismology and Earth Structure





0 300

98 E
90 94

Figure 11.26. Image of radiated seismic energy as a function of position along the arc for
the 2004 Sumatra“Andaman earthquake obtained from stacking of high frequency energy
at the Japanese Hi-Net Array by Ishii et al. (2005). The changes in the intensity of the
radiated energy have a close correspondence to structural features in the slab.

For a truly great earthquakes such as the Mw 2004 Sumatra-Andaman event,
with a fault length exceeding 1300 km, the pattern of slip suggests a domino effect
whereby the stress induced by rupture on one fault segment triggers slip on the
next. The source time function includes at least two such stutters (Ammon et al.,
2005), which correspond to distinct geographic points in the earthquake rupture.
This great thrust event occurs where the Indo“Australian plate descends beneath the
Asian plate and the Andaman micro-plate. The variable energy release along the
Sumatra“Andaman arc is illustrated in Figure 11.26. It is likely that segmentation
of the slip is associated with the prior structure of the subduction zone, through
changes in the dip, strike and physical properties of the seismogenic part of the
slab interface.
Lithospheric Deformation

12.1 De¬nitions of the lithosphere
The concept of a lithosphere arises from the presence of a zone of strength near the
Earth™s surface that is capable of transmitting and resisting stress. A consequence
of the complex rheology of the Earth is that the apparent thickness of this zone
becomes smaller as the time scale of the processes being considered becomes
longer. This has lead to a variety of de¬nitions of lithosphere, based on particular
classes of observations, that are not necessarily mutually compatible.
We can think of the mechanical lithosphere as representing the outer part of the
Earth where stress can be transmitted over geological time scales. This is linked to,
but not identical with, the thermal lithosphere, which is that region where thermal
energy is largely transferred by heat conduction (Figure 12.1).

β „M
(a) (b) (c)
1280 C

boundary layer
100 100 100

thermal boundary layer

200 200 200
adiabatic gradient
z z z

Figure 12.1. Relation of different concepts of lithospheric thickness illustrated for a conti-
nental situation: (a) thermal concept through the transition from a conductive geotherm to
an adiabat, (b) seismic de¬nitions based on the properties of the shear wavespeed distribu-
tion with depth, and (c) the viscous relaxation time as a function of depth.

In seismological studies the lithosphere is commonly taken to be associated with
the region of elevated seismic wavespeeds (the ˜lid™), with a base assigned either in
the region of maximal vertical velocity gradient, or on entry into a zone of lowered
shear wavespeed. Estimates of the temperature distribution associated with seismic

258 Lithospheric Deformation

wavespeed (e.g., McKenzie, Jackson & Priestley, 2005) suggest that this occurs
in the same neighbourhood as the thermal transition (Figure 12.1). Although
the seismic lid may have a relatively sharp base beneath the oceans, there are
only a few observations suggesting very rapid velocity gradients in the continental
A useful viewpoint for understanding the issue of the dependence of the apparent
thickness of the lithosphere on the time scale of the phenomenon comes from the
Burgers model for a viscoelastic medium introduced in Section 4.6. Such a system
displays (i) instantaneous elasticity that can be associated with seismological
observations, (ii) intermediate term viscoelastic response that can be used to
describe the rebound after the removal of ice loading at the end of the ice age,
and (iii) long term viscous ¬‚ow. The time scale for the viscous process „M varies
substantially with depth as indicated in Figure 12.1c.

β · E
(a) (b) (c)

100 100

200 200

z z

Figure 12.2. Relation of models for lithospheric properties for different time scales rep-
resenting mechanical equivalents for different aspects of rheology: (a) seismic, (b) glacial
rebound, (c) effective elastic thickness for long term loads.

Much of the properties of lithospheric deformation can be understood by simple
mechanical models such as the ¬‚exure of an elastic plate to describe long term loads
or the simple viscous ¬‚uid model of glacial rebound in Section 7.4. As a result
simpli¬ed models have been widely used to extract proxy quantities to characterise
the more complex actual physical situation. The relation of such lithosphere models
for different time scales is indicated schematically in Figure 12.2.

12.2 Thermal and mechanical structure
12.2.1 Thermal conduction in the oceanic lithosphere
Oceanic lithosphere is formed at mid-ocean ridges, indicated by the double lines
in Figure 12.7, and then is carried away from the ridges by the movement of the
oceanic plates. Hot material from the mantle emerges at the axes of the mid-ocean
ridges and forms juvenile oceanic crust, further basaltic material solidi¬es onto
12.2 Thermal and mechanical structure 259

the base as the new oceanic lithosphere moves away from the axis. Hydrothermal
circulation is important through the rapidly cooled materials created just at the
ridge, and heat ¬‚ux measurements are rather scattered. Heat ¬‚ow stabilises for
older lithosphere as hydrothermal effects become less important. As the lithosphere
moves away from the mid-ocean ridge it cools and contracts with an increase of
depth that is roughly proportional to the square root of age out to 65 Ma or so, but
then deepens rather more slowly.
Much of this behaviour can be captured with relatively simple thermal models.
Consider a two-dimensional situation, independent of the y-coordinate, with
asthenospheric material ascending along the ridge axis at x = 0 with constant
temperature TM. With a plate transport speed vx in the x-direction, the temperature
is governed by (8.4.3)
‚2T ‚2T
‚T ‚T
ρCp + vx =k + 2. (12.2.1)
‚t ‚x ‚z
In a steady state ‚T/‚t = 0, and we can make a further simpli¬cation if we assume
that the conduction of heat in the horizontal direction is negligible compared with
the advective contribution from vx‚T/‚x and vertical conduction. Thus, neglecting
the ‚2T/‚x2 term in (12.2.1), we have
k ‚2T ‚2T ‚T
= κH 2 = vx . (12.2.2)
ρCp ‚z2 ‚z ‚x
We now introduce the spreading time ts so that x = vxts, and then we recover
the simple thermal conduction equation, in terms of ts
‚2T ‚T
κH 2 = . (12.2.3)
‚z ‚ts
The boundary conditions are T = TM at x = 0, and the top surface is maintained
at zero temperature (T = 0 at z = 0); the requisite solution is provided by (8.4.28),
so that
T (z, ts) = TM erf √ , (12.2.4)
with associated surface heat ¬‚ow
‚T kTM
= ’√
Q0(ts) = ’k . (12.2.5)
‚z z=0 πκts
Thus for this half-space cooling model, the surface heat ¬‚ux is proportional to

1/ ts. Isotherms, surfaces of constant temperature, will correspond to constant

argument of the error function zc = 2c κts, and hence the depth to a given
temperature increases as the square root of the lithospheric age. It is conventional
to take the base of the lithosphere as speci¬ed by such an isotherm, e.g., 1300—¦ C.
The cooling lithosphere shrinks and the depth of the seabed can be estimated by
assuming full isostatic compensation (Figure 12.3). We assume a compensation
260 Lithospheric Deformation




Figure 12.3. Oceanic depth pro¬le associated with lithospheric cooling with full isostatic

depth Hc in the mantle beneath the level of the deepest oceanic lithosphere and
then require the mass in any vertical column extending to Hc to be the same.
We take the water depth to be w0 at the ridge axis and the depth origin at the
sea¬‚oor at the axis. The vertical traction at the level Hc at the axis
σzz(x = 0, Hc) = gρww0 + g dz ρM, (12.2.6)
where ρw is the density of sea water and ρM is the density of the hot asthenosphere.
Off the ridge axis, for lithosphere thickness L and additional water depth w,
σzz(x, Hc) = gρw(w0 + w) + g dz ρ(z)
w+L Hc
= gρw(w0 + w) + g dz ρ(z) + g dz ρM. (12.2.7)
w w+L
For the tractions to be equal at different horizontal positions we require
ρww0 + ρMHc = ρw(w0 + w) + dz ρ(z) + ρM(Hc ’ w ’ L), (12.2.8)
and so the excess water depth is given by
w(ρM ’ ρw) = dz (ρ(z) ’ ρM). (12.2.9)
The density pro¬le ρ(z) is controlled by the temperature variation through the
effects of thermal contraction,
ρ(T ) = ρM[1 ’ ±th (T ’ TM)] (12.2.10)
in terms of the thermal expansion coef¬cient ±th , and so we can evaluate the integral
in (12.2.9) to give
w(ρM ’ ρw) = ρM dz ±th (TM ’ T (z)). (12.2.11)
12.2 Thermal and mechanical structure 261

With the temperature distribution (12.2.4), the excess water depth
ρM z

w(ts) = dz ±th TM 1 ’ erf , (12.2.12)
ρM ’ ρw 4κts

Once we are away from the axis the integral in (12.2.12) can be approximated by
its asymptotic value as L ’ ∞, (4κts/π)1/2, and so
2±th ρMTM κts 2
w(ts) ≈ . (12.2.13)
ρM ’ ρw π
We can improve on the simple half-space model by taking a plate model in
which the base of the rigid oceanic lithosphere derived from cooled asthenosphere
is de¬ned by an isotherm. In such a plate model the base of the lithosphere at
z = L is assumed to have the same temperature TM as at the vertical ridge axis.
We introduce the thermal Reynolds number R = vxL/2κ and can then ¬nd a series
solution to (12.2.1). In the steady state when ‚T/‚t = 0 we obtain a similar form
to (8.4.29)

z 2 nπz vxts
exp ’[(R2 + n2π2)1/2 ’ R]
T (z, ts) = TM + .
L nπ L L
The isotherms initially deepen as the square root of age, but ¬‚atten out for older
lithosphere as the steady-state thermal structure approaches a linear gradient in
depth. The excess water depth can again be estimated from (12.2.11) using the
isostatic assumption and has a similar age dependence to the isotherms.

12.2.2 Mechanical deformation
The lithosphere can be represented as a layered composite of differing rheological
properties underlain by viscoelastic material with reduced viscosity associated with
the asthenosphere. As noted in Section 12.1, the differing rheological attributes
manifest themselves through a variation in resistance to load depending on the time
scale of the imposed deformation (see Figure 12.1).
In the crust at shallow depth the dominant response is elastic, but the material
is brittle, with ¬nite strength, and failure occurs in localized zones. With
increasing depth we pass into the ductile regime where ¬‚ow is distributed with
a temperature dependent relation between stress and strain, such as (9.3.3). The
relative contributions of the different styles of rheological properties depend on a
range of factors such as mineral properties and chemistry, pore geometries, stress
conditions, and lithostatic and pore pressure. A useful summary of the nature of the
transition from the brittle to ductile and plastic regimes and associated experimental
evidence is provided by Evans & Kohlstedt (1995).
For deformation at seismic periods the dominant response comes from
262 Lithospheric Deformation

water or air

upper crust
lower crust


Figure 12.4. Model used to calculate crustal response to surface and internal loading.
Loads are imposed: (1) at the surface by imposing a layer of thickness s1 (x) and den-
sity ρu , (2) at the interface between upper and lower crust by moving the interface keeping
the crustal thickness constant, and (3) at the Moho by adding a layer with a thickness s3 (x)
and density ρm .

instantaneous elasticity; whereas on the scale of millennia, associated with loading
and unloading of ice in the glacial era, both the elastic and viscous components
come into play. For long term loads, such as the support of major topography with
time scales of many millions of years, the focus is again on elastic properties, but
now for the upper part of the lithosphere in the crust and, perhaps, the uppermost
A relatively simple model that can capture many aspects of crustal de¬‚ection
under long-term loads is provided by a thin sheet model with a laminated crust as
in Figure 12.4. The upper crust is underlain by a higher density lower crust and
loading is applied at the surface or Moho by the addition of layers with appropriate
density, and in the interior by moving the interface between the upper and lower
crust. The de¬‚ection w resulting from a load is assumed to be independent of
depth in the thin sheet approximation.
When a layer of thickness s1(x) is added to the surface the resulting de¬‚ection
w1 is given by

D∇4 + g(ρm’ρw) w1 = ’g(ρu ’ρw)s1, (12.2.15)

where D is the ¬‚exural rigidity of the plate, g the acceleration due to gravity, ρm
the density of the mantle, ρu the density of the upper crust and ρw is the density
of the ¬‚uid (water or air) overlying the crust. The hydrostatic restoring force
g(ρm ’ρw)w1 arises because mantle with density ρm has been replaced through
the de¬‚ection w1 by the upper ¬‚uid of density ρw. The ¬‚exural rigidity of the plate
D is related to the elastic thickness Te by
D= , (12.2.16)
12(1 ’ …2)
where E is Young™s modulus, and … is Poisson™s ratio. A simple derivation of
(12.2.16) is provided by Watts (2001) §3.4.
12.2 Thermal and mechanical structure 263

If we take the Fourier transform of (12.2.15) with respect to the horizontal
coordinate x we obtain
Dk4 + g(ρm’ρw) w1 = ’g(ρu ’ρw)¯1,
¯ (12.2.17)
where the transformed variables are indicated by overbars and k is the horizontal
wavenumber. The net surface topography e1(x) is the combination of the imposed
load and the de¬‚ection so that
e1(x) = s1(x) + w1(x). (12.2.18)
The resulting gravity anomaly ”g1 for wavenumber k is given by (McKenzie,
”g1 = 2πG{(ρu ’ρw)e1 + (ρl ’ρu)e’ktu w1 + (ρm’ρl)e’ktc w1}, (12.2.19)
¯ ¯ ¯ ¯
where G is the gravitational constant and, as shown in Figure 12.4, tu is the
thickness of the upper crust and tc that for the whole crust. From (12.2.17) and
(12.2.17) we have
Dk4 + g(ρm’ρu) w1 = ’g(ρu ’ρw)e1,
¯ ¯ (12.2.20)
and we can therefore express the gravity anomaly ”g1 in terms of the surface
topography e1 through an admittance Z1,
”g1 = Z1e1 = 2πG(ρu ’ρw)(1 + Z1)e1,
¯ ¯ ¯ (12.2.21)
where Z1 is given by
(ρl ’ρu)e’ktu + (ρm’ρl)e’ktc
Z1 = g . (12.2.22)
Dk4 + g(ρm’ρu)
When the load is at the interface between the upper and lower crust the de¬‚ection
w2(x) is determined by

Dk4 + g(ρm’ρw) w2 = ’g(ρl ’ρu)¯2,
¯ (12.2.23)
with a corresponding gravity anomaly
”g2 = 2πG{(ρu’ρw)w2 +(ρl’ρu)e’ktu (¯2+w2)+(ρm’ρl)e’ktc w2}.(12.2.24)

¯ ¯ ¯
The surface topography induced by this internal load e2 = w2 and the gravitational
¯ ¯
anomaly can again be expressed in terms of an admittance function Z2,
”g2 = Z2e2 = 2πG(ρu ’ρw)(1 + Z2)e2,
¯ ¯ ¯ (12.2.25)
Dk4 + g(ρm’ρw) ’ g(ρl ’ρu) ’ktu ρm’ρl ’ktc
Z2 = + . (12.2.26)
e e
g(ρu ’ρw) ρu ’ρw
For an internal load placed at the Moho
”g3 = Z3e3 = 2πG(ρu ’ρw)(1 + Z3)e3,
¯ ¯ ¯ (12.2.27)
264 Lithospheric Deformation

where now
Dk4 + g(ρl ’ρw ’ktc
ρl ’ρu ’ktu
Z3 = + .
e e (12.2.28)
ρu ’ρw g(ρu ’ρw)
The total surface topography will represent the sum of the effects of the three
different types of loading and the total gravitational anomaly will again be the sum
of the gravitational effects produced by the three classes of loads
3 3 3
e= ei, ”g = ”gi = Ziei.
¯ ¯ ¯ ¯ ¯ (12.2.29)
i=1 i=1 i=1
As in the work of Forsyth (1985), it is not unreasonable to assume that the surface
and internal loads are uncorrelated with each other; in which case
e1e— = e1e— = e3e— = 0,
¯ ¯2 ¯ ¯3 ¯ ¯2 (12.2.30)
where an asterisk indicates a complex conjugate and the angle brackets indicate
averages taken over a wavenumber band of width ”k centred on k.
The admittance Z relates the surface topography to the gravitational anomaly
”g = Ze + n
¯ ¯¯ (12.2.31)
where n represents that part of the gravity anomaly ”g that is not coherent with e
¯ ¯ ¯
such as that due to loads that have no surface topographic expression. Thus
¯ ¯—
”g e— i=1 Zi eiei
Z= = , (12.2.32)
ee— 3
¯¯ ¯ ¯—
i=1 eiei
no cross terms appear because we have assumed that the loads are uncorrelated
with each other (12.2.30).
The de¬‚ections due to a surface load can be determined directly from the solution
of (12.2.15). The general solution takes the form,
w = eζx[A cos(ζx) + B sin(ζx)] + e’ζx[C cos(ζx) + D sin(ζx)], (12.2.33)
where A, B, C and D are integration constants. The ¬‚exural parameter ζ controls
the amplitude and wavelength of the deformation
ζ= . (12.2.34)
If we impose a line load P at x = 0, the solution in x > 0 cannot contain growing
terms and, from symmetry dw/dx = 0 at x = 0, so
w = w0e’ζx[cos(ζx) + sin(ζx)], (12.2.35)
and from force balance w0 = Pζ/2g(ρm ’ ρw). The de¬‚ection vanishes when
cos(ζx) + sin(ζx) = 0 and thus when ζx = 3π/4, 7π/4, · · · . In addition to the
main downward de¬‚ection beneath the load, there is an upward peripheral bulge on
either side of the central depression.
12.2 Thermal and mechanical structure 265

Watts (2001) §3.5.2 provides a detailed discussion of loading applied at or near
a break using a formulation in terms of a semi-in¬nite beam for which
w = w0e’ζx cos(ζx), (12.2.36)
with the maximum de¬‚ection w0 = 2Pζ/[g(ρm’ρw)]. There is again a bulge, but
it is both narrower and of larger amplitude than for the continuous plate.
The thin-plate theory is based on the assumption that the thickness of the shell
is thin compared with the radius of curvature of any de¬‚ection, and that the stress
component σzz in the vertical direction can be neglected by comparison with the
other stress components.

12.2.3 Estimates of the elastic thickness of the lithosphere
The thin plate formulation of the de¬‚ection of the lithosphere has been used
by many authors to obtain estimates of the equivalent elastic thickness of the
lithosphere by exploiting gravity anomalies, either by spectral methods or using
forward modelling.
In the forward modelling approach the gravity anomaly due to a load, e.g. the
surface load due to the topography of an oceanic island or a mountain chain, and the
¬‚exural compensation arising from the de¬‚ection of an elastic plate are calculated.
The estimates of the gravity anomalies for different elastic plate thicknesses Te are
then compared with the observations and the value of Te that secures the best ¬t
between the observed and calculated gravity anomalies is selected.
The spectral methods are based on estimates of the admittance Z (12.2.32)
between the gravity anomalies and surface topography. In principle it should be
a relatively simple matter to determine the admittance through spectral division in
the spatial Fourier domain, e.g.,
”g(k)H— (k)
Z(k) = = (12.2.37)
H(k)H— (k)
where ”g(k) is the transform of the gravity anomaly and H(k) the transform of
the topography. However, in practice the situation is complicated by noise in the
data and limitations in the geometry of data collection that are most pronounced
in the distribution of ship tracks in the ocean. It is common practice to work with
auxiliary information such as the coherence:
C(k)C— (k)
with Eg(k) = ”g(k)”g— (k).
γ (k) = (12.2.38)
The coherence function derived from the observations is then compared directly
with model predictions. At short wavelengths the coherence approaches 0 because
the lithosphere appears strong irrespective of the value of Te. Such small features
would cause only minor ¬‚exure with little contribution to the gravity anomaly. At
long wavelengths the coherence normally approaches 1, because the lithosphere has
266 Lithospheric Deformation

¬‚exed under load. The presence of buried loads complicates the situation (Forsyth,
1985; McKenzie, 2003).
The two different approaches would be expected to give similar results for the
effective elastic thickness of the lithosphere Te and this is certainly true for the
oceans. For example forward modelling for the gravity anomalies associated with
the Hawaiian“Emperor seamount chain leads to an estimate of Te for the Paci¬c
Plate of 25 ± 9 km, whereas the spectral approach using free-air anomalies yields
values of 20“30 km. In particular, the various Te estimates provide the same
dependence on age of the lithosphere at the time of loading: Te increases from
2“6 km for young lithosphere to around 30 km for old lithosphere. Such values are
somewhat thinner than the estimates of lithospheric thickness in the oceans derived
from seismic observations, but re¬‚ect very different time scales for loading.
In the continents there are greater discrepancies between the two approaches
to estimating Te. The results of direct modelling are often smaller than those
obtained by spectral transfer function methods. Improved techniques for spectral
estimation (such as multi-taper techniques) have improved the spatial resolution
of Te. Forsyth (1985) introduced an allowance for the presence of internal loads
using the coherence with the Bouguer anomalies, where corrections are made for
the gravitational attraction of all material above sea level. However, McKenzie &
Fairhead (1997) have argued that the Bouguer approach provides upper bounds
on Te rather than the true value. They proposed instead a free-air admittance
method to argue that the effective elastic thickness of the continental lithosphere
was less than 25 km. Subsequently McKenzie (2003) produced a re¬ned coherence
method to allow for the presence of internal loads, particularly when there is no
associated surface topography. This approach is based on the model presented
above (see Figure 12.4), but intriguingly the effective elastic thickness Te derived
by this method is smaller than the crustal thickness.
The forward modelling approach is commonly applied to foreland basins and in
some cases such as, e.g., the Ganges basin has indicated Te signi¬cantly in excess
of the local crust thickness, 70 km compared with a crust of 40 km (Burov & Watts,
2006). Such a result suggests the presence of elastic strength in the mantle. The
spectral methods have generally been derived for old cratons, but results have also
been derived for orogenic belts and rifts. There is increasing evidence for relatively
rapid horizontal variations of Te on scales of a few hundred kilometres and also
of anisotropy in the gravitational response (e.g., Simons et al., 2000). Yet, the
theoretical development is based on a uniform isotropic plate, and so the values of
Te should be regarded as useful comparators between regions rather than having
direct physical signi¬cance.

12.2.4 Strength envelopes and failure criteria
The clearest manifestation of the ¬nite strength of the lithosphere comes from
the presence of shallow earthquakes. It is rare for a major event to break fresh
12.2 Thermal and mechanical structure 267

rock, so that the behaviour observed is dominantly that of a medium comprised
of units linked together by the strong friction at the boundaries of the blocks.
Experimental studies and computer simulations indicate that the friction laws are
dependent on prior deformation history and slip-rate, so that strength can drop
dramatically once relative motion starts. The friction is suf¬cient that for processes
with large horizontal scale, e.g., earth tides and glacial rebound, the crust behaves as
a single material. Yet the stress accumulating at places like convergent margins has
eventually to be released, as in the great 2004 Sumatra-Andaman earthquake where
a thrust front migrated over 1300 km along the plate boundary in a few minutes.
The passage of a large earthquake on a fault redistributes local concentrations
of stress that in turn are relieved by aftershocks, but can also change the stress
distribution in the neighbourhood of another fault system, rendering that zone
closer to the stress conditions required for failure.
The changing conditions from brittle failure at shallow depth to ductile ¬‚ow
at depth can be summarised though a yield strength envelope that attempts to
provide a single strength pro¬le through the lithosphere, based on the dominant
constitutive law at each depth. As we have seen in Chapter 9, the details
of the stress“strain (rate) relations depend strongly on material properties and
temperatures, so that a single strength envelope for a particular regime is an
oversimpli¬cation. Further, estimates of strength envelopes are based on data
from experimental rock mechanics at much higher strain rates (10-4 “10-6 s-1 ), than
those in the Earth which rarely exceed 10-13 s-1 and where slower process such as
dynamic recrystallisation might be important.
For the upper portion of the lithosphere, brittle deformation by frictional sliding
on pre-existing fracture planes of suitable orientation can be described through
the use of a Mohr“Coulomb fracture criterion (cf. Section 10.2.3), in these
circumstances frequently known as Byerlee™s law. The critical stress at which
frictional sliding takes place σsm is then related to the normal stress σn and the
ratio »f between the pore ¬‚uid pressure pf and the lithospheric pressure plith ,
»f = pf/plith , by
σsm = (M1σn + M2)(1 ’ »f), (12.2.39)
where M1 and M2 are empirical constants, equivalent to μ and c in (10.2.6). The
relation (12.2.39) incorporates the shift in effective stress state due to pore ¬‚uid
pressure. For materials in compression, Brace & Kohlstedt (1980) suggest the use
M1 = 0.85, M2 = 0 MPa, σn < 200 MPa,
M1 = 0.6, M2 = 60 MPa, σn > 200 MPa.
In tension, the differential stress needed to produce failure is signi¬cantly reduced
(see Figure 10.10).
When material is fractured or porous, the controlling factor for brittle strength is
pore ¬‚uid pressure pf. In the upper crust a hydrostatic pore pressure gradient with
268 Lithospheric Deformation

a pore pressure coef¬cient »f = 0.4 works quite well. In the deep KTB borehole in
Germany, hydrostatic pore pressure gradients persist to 9 km depth in metamorphic
basement rock.
In the ductile regime, the creep can be described by the grain size dependent
relation (9.3.3), that can be recast in a form where the stress σ is related to the
strain-rate ™ as
n m
1 d0 TM
σ=B exp ’g , (12.2.41)
™G d T
for grain size d and temperature T . Since the exponent n is generally greater than
1, the stress dependence on ( ™ )’n is non-linear and highly temperature dependent.

Oceanic lithosphere
For the oceanic lithosphere the dominant mineralogy is olivine, and away from
the ridge crest, the brittle regime extends beneath the oceanic crust. At greater
depths the behaviour is ductile and the overall stress envelope is determined by
the intersection of the linear yield stress relation from brittle failure and the
nonlinear ductile law. At the brittle“ductile transition equivalent differential stress
is predicted for the two different regimes. The resulting yield strength envelope,
indicating the stress that can be supported by the lithosphere before it yields, has a
con¬guration similar to a sail (Figure 12.5). The area inside the strength envelope
can be regarded as a measure of the integrated strength of the lithosphere.

Differential Stress [MPa]
-1000 0 1000 -1000 0 1000
0 0
Tension Compression Tension Compression

20 20
Depth [km]

Depth [km]

40 40

60 60
16 Ma 100 Ma OCEANS

80 80

Figure 12.5. Yield strength envelopes for oceanic lithosphere, based on a strain rate of
10-14 s-1 and a half-space cooling model for the oceanic geotherm.

The typical strength envelope behaviour for young lithosphere (16 Ma) and older
lithosphere (100 Ma) is shown in Figure 12.5 using a strain rate of 10-14 s-1 and a
half-space model for the geotherm in the oceanic lithosphere (Section 12.2.1). For
ages less than 100 Ma there would be little change using a cooling plate model.
12.2 Thermal and mechanical structure 269

The brittle“ductile transition corresponds approximately to the 450—¦ isotherm in
the oceanic lithosphere, and hence roughly to the estimates of Te that show good
correlation with the age of the oceanic lithosphere.
The majority of oceanic earthquakes are relatively shallow, but a few occur in
the oceanic mantle, and are con¬ned at depths above the 650—¦ C isotherm. The
maximum depth of earthquakes is thus somewhat deeper, but generally consistent
with the apparent elastic thickness of the lithosphere Te. Deeper events at the
outer-rise of subduction zone trenches tend to indicate a compressional regime
whilst the shallower events indicate tension. The separation in mechanism is not so
distinct for the much rarer events within the oceanic plates.

Continental lithosphere
The details of the yield strength behaviour depend strongly on the particular
mineralogy of the highly heterogeneous continental crust. A typical feature is a
signi¬cant reduction in strength below 20 km depth, because the onset of ductility
in quartz-rich rocks occurs at a shallower depth than for olivine under the same
thermal gradients. There is the possibility of an increase in strength in the lower
crust, and again beneath the Moho particularly when the uppermost mantle material
is dry (Figure 12.6).
Evidence for a strong uppermost mantle has been adduced from occasional
earthquakes below the Moho, but there could well be a residual in¬‚uence of the
assigned depth (33 km) used by international agencies for earthquakes for which a
reliable depth cannot be found. Reinterpretation of possible mantle events at 70“85
km beneath Tibet, in light of improved constraints on crustal structure, suggests
that they lie in the lower crust (Jackson, 2002), and may be associated with the
transition from metastable granulite to eclogite in the presence of water.
Even very small amounts of water can have a dramatic effect on the strength
distribution in a material, as can be seen from Figure 12.6. The conventional
viewpoint has favoured a scenario with a wet lower crust and a dry uppermost
mantle, in which case a considerable portion of the total strength of the lithosphere
would reside in the mantle. This viewpoint has been challenged by Jackson (2002)
building on the observations of deep lower crustal earthquakes and estimates of the
effective elastic thickness Te for the continents as less than that of the seismogenic
layer. Jackson suggests that the normal continental scenario is both wet lower crust
and uppermost mantle so that nearly all strength would reside in the seismogenic
upper crust. In some regions, the lower crust could be dry and sustain earthquakes
as in the Indian Shield material in the doubled crust under Tibet. Certainly most of
the infrequent earthquakes in shield areas initiate above 15 km depth; though there
is some evidence for lower crustal events in the Fennoscandian region associated
with glacial unloading.
The concept that strength resides solely in the upper crust has been disputed by
Burov & Watts (2006) who favour the retention of a strength contribution from
the uppermost mantle. Simulations of lithospheric dynamics with weak or strong
270 Lithospheric Deformation

Differential Stress [MPa]
00 1000 2000 3000


Temperature [ C]
Depth [km]



Figure 12.6. Yield strength envelope, in compression, for the continental lithosphere as a
function of depth at a strain rate of 10’15 s’1 and a surface heat ¬‚ow of 60 mW m’2 , based
on experimental results for dry materials, and those containing small amounts of water. Dry
conditions are indicated with grey lines and wet with black. The upper crust is taken as a
wet quarzite, the lower crust as either a dry diabase (D) or undried granulite (W) and the
upper mantle as wet or dry olivine. A wide range of strength conditions can be produced
by different combinations of wet and dry elements. When both lower crust and mantle are
wet, strength resides in the seismogenic upper crust. In contrast, with a dry lower crust,
deep crustal earthquakes would be possible. If the whole crust is wet then most strength
resides in the uppermost mantle. [Courtesy of J. Jackson.]

uppermost mantle (but similar effective Te) show very signi¬cant differences in
the stability of mantle structures. A weak uppermost mantle tends quite rapidly
(<5 Myr) towards the development of drip structures sinking into the mantle,
and the collapse of surface topography. Strength in the uppermost mantle allows
topography induced by orogeny to be sustained for 100 Myr or more.
A major dif¬culty with resolving such fundamental issues is that most of the
available information that bears on strength comes from regions that are clearly
abnormal, because of the occurrence of sizeable earthquakes. The widespread
development of seismic anisotropy in the upper part of the mantle lithosphere that
shows preferred directions linked to past tectonic history provides evidence for
12.3 Plate boundaries and force systems 271

preservation of structures over substantial time periods, which would be dif¬cult to
sustain without some strength in this zone.

12.3 Plate boundaries and force systems
12.3.1 Nature of plate boundaries
The lithosphere is suf¬ciently resilient to transmit stress over long distances and
although there are some regions of internal deformation marked by sporadic
earthquakes, the main features of plate tectonics are dominated by convergent and
divergent plate boundaries (see Figure 12.7). The divergent boundaries mark the
mid-ocean ridges where upwelling of hot material leads to the creation of new
ocean crust and are shown by a double line in Figure 12.7. The rates of divergence
vary from less than 10 mm/yr to more than 90 mm/yr, and in some cases there is
some asymmetry in spreading rate on the two sides of the divergent boundary. The
relative motion of the relatively rigid plates has to be accommodated by jogs in the
divergent margins associated with transform faults. The traces of the transforms
form small circles about the pole of rotation that describes the motion of the two
separating plates. Most of the oceanic earthquakes occur on or near the transform

Figure 12.7. The earthquake distribution across the globe in relation to plate boundaries.
Earthquakes below 150 km deep associated with subduction are shown in slightly darker
tone. Divergent margins are indicated with a double line and convergent plate boundaries
with a heavier line weight.

The addition of more oceanic material has to be accommodated by the
re-assimilation of former oceanic lithosphere into the mantle at subduction zones,
which are the major form of convergent margins. The convergent margin segments
around the globe are indicated with heavier lines in Figure 12.7, and it is
272 Lithospheric Deformation

immediately apparent that they are accompanied by the majority of the world™s
earthquakes. The consumption of oceanic lithosphere in front of a continental
component on the plate can lead to the scenario where, as under Tibet, an attempt
is being made to force relatively buoyant continental lithosphere to depth.
The plate boundaries are not static but evolve over time, and only a limited
number of possible triple junctions between plates are stable. When a plate
overrides a speading centre, as appears to have happened for the extinct Farallon
plate in North America, the resulting motions need still to be taken up and
the modern San Andreas transform fault system links subduction in Mexico and
The modern-day con¬guration of plates does not have each plate with convergent
and divergent boundaries. The African plate has no convergent margins, and the
Philippine plate no current ridge system. The balance of forces associated with
the motion of the plates will thus depend strongly on the geometry of the plate

12.3.2 Plate boundary forces
Forsyth & Uyeda (1975) provided a useful description of the nature of the forces
acting at the boundary of the lithospheric plates (Figure 12.8). As we have just seen
the oceanic lithosphere is suf¬ciently strong to transmit stress over long distances
and hence for the entire plate to behave as a single entity.


. 10
( 16)