. 4
( 16)


t0 V
t t
dm (giu— uig— ) dV (‚jσiju— ’ ui‚jσ— ).
= ’ +
dt dt
i i i ij
t0 V t0 V

The volume integral over the stress

dV(‚jσiju— ’ui‚jσ— ) = dV‚j(σiju— ’uiσ— ) ’ dV(σije— ’eijσ— ),
i ij i ij ij ij


σije— = cijklekle— = cklije— ekl = σ— ekl. (5.3.5)
ij ij ij kl
80 Linearised Elasticity and Viscoelasticity

Thus (5.3.3) can be written as
‚ui —
ui ’ ui i
‚t ‚t
V t0
t t
dm(giu— uig— ) dS(„iu— ’ ui„— ).
= ’ +
dt dt (5.3.6)
i i i i
t0 V t0 S

5.3.1 The reciprocal theorem
Take the time integration up to in¬nity and assume that u, u— are such that the
left-hand side of (5.3.6) vanishes, then
∞ ∞
dm giu— dS „iu— dm g— ui + dS „— u— ,
+ =
dt dt
i i i ii
’∞ V S ’∞ V S
and for a static ¬eld we can drop the time integral

dm giu— + dS „iu— = dm g— ui + dS „— u— . (5.3.8)
i i i ii

This is Betti™s reciprocal theorem.
In applications of Betti™s reciprocal theorem (5.3.8) we take the required ¬eld to
be u and choose an appropriate starred ¬eld for which the solution is known and
which has as close a character as we can achieve to the desired ¬eld.

(i) In¬‚uence of self-gravitation
Self gravity is switched on in a homogeneous sphere “ ¬nd the decrease in radius:
At the surface r = a, „ = 0. Let g be the acceleration due to self gravity, then the
interior body force per unit mass is g = ’gx/a.
Choose u— = x implying uniform σ— = 3κ δij , then g— = 0 because the ¬eld is
self-equilibrated. From (5.3.8)
· xρ dV = 3κ 4πa2 u(a),


and hence
’gρ 4π 5
a = 3κ 4πa2 u(a),
i.e., the decrease in radius is ’u(a) = ρga2 /(15κ). •
(ii) Loaded cylinder
A solid circular cylinder of height h and radius a is subjected to compressive loads L
applied via rigid plates in contact with its ends. Radial displacements over the plates are
opposed by a uniform frictional force „ (Figure 5.1). Find the reduction in height δh when
L is large enough to cause sliding everywhere over the ends.
5.3 Integral representations 81




. h Figure 5.1. Con¬guration of loaded
cylinder and force distribution


Let x3 be measured axially from the centroid, and choose
u— = (…x1 , …x2 , ’x3 )p— /E,
corresponding to a uniaxial pressure p— , so g— = 0. Now also g = 0 in V, „— = „ = 0
over the curved surfaces. Thus, denoting an end by A,

„ · u— dA = „— · u dA,

the load L = ’ „3 dA, and so
r„2πr dr = p— 1 πa2 δh.
2 HL ’… 2
E 0

h L 4…a„
δh = ’ .•
E πa 3h

5.3.2 The representation theorem
Consider a delta-function body force

gi = δikδ(x ’ ξ)δ(t ’ θ), (5.3.9)

acting in a volume V. The resulting displacement at x, t is speci¬ed by the Green™s
tensor with components

Gk(x, t; ξ, θ), (5.3.10)
82 Linearised Elasticity and Viscoelasticity

associated with the force in the k-direction at ξ, θ. Even when the properties of the
region vary Gk will depend on the time difference t ’ θ.
If we take gi as above and set
g— = δilδ(x ’ ζ)δ(σ ’ t), (5.3.11)

the displacement
u— = Gl(x, ’t; ζ, ’σ). (5.3.12)
i i

We assume that both Gk and Gl obey the same homogeneous boundary conditions
i i
on S, and then when we apply the reciprocal theorem (5.3.7) the surface integrals

dV δikδ(x ’ ξ)δ(t ’ θ)Gl(x, ’t; ζ, ’σ)
dt i
’∞ V

dV δikδilδ(x ’ ζ)δ(σ ’ t)Gk(x, t; ξ, θ),
= dt (5.3.13)
’∞ V
so that when ξ, ζ are in V
Gl (ξ, ’θ; ζ, ’σ) = Gk(ζ, σ; ξ, θ).
kξ lζ (5.3.14)
Equation (5.3.14) expresses Green™s tensor reciprocity.
Now set
g— = δikδ(x ’ ξ)δ(θ ’ t) (5.3.15)

in (5.3.6) and then
u— = Gk(x, ’t; ξ, ’θ), (5.3.16)
i i

and so
dVx{gi(x)Gk(x, ’t; ξ, ’θ) ’ ui(x)δ(x ’ ξ)δ(θ ’ t)}
0= dθ i
t0 V

dSx{„i(x)Gk(x, ’t; ξ, ’θ) ’ ui(x)Hk(x, ’t; ξ, ’θ)}
+ dθ i i
t0 S

dVx[‚tui(x)Gk ’ ui‚tGk]t=t0 ,
+ (5.3.17)
i i

where Hk is the traction associated with Gk.
i i
Now, using the reciprocity of the Green™s tensor we can obtain a representation
of the displacement as an integral over the boundaries, body forces and initial
dVx gi(x)Gi (ξ, t; x, s)
˜(ξ)uk(ξ, t) =
ξ ξ kξ
t0 V

dSx{„i(x)Gi (ξ, t; x, s) ’ ui(x)Hi (ξ, t; x, s)}
+ kξ kξ
t0 S

dVx[‚tui(x)Gi ’ ui‚tGi ]t=t0 .
+ (5.3.18)
k k
5.4 Elastic waves 83

˜(ξ) = 1,
ξ ξ in V,
= 0, otherwise .

5.4 Elastic waves
The equation of motion and the stress-strain relation for small disturbances in an
elastic medium, in the absence of body forces, can be combined to give
‚ ‚uk
cijkl =ρ 2 , (5.4.1)
‚xj ‚xl ‚t
as the governing differential equation for displacement u. Equation (5.4.1) controls
the spatial and temporal development of the displacement ¬eld and admits solutions
in the form of travelling waves. Consider then a plane wave travelling in the
anisotropic medium with frequency ω:
ui = Ui exp[iωpn · x ’ iωt]
= Ui exp[iωpnkxk ’ iωt].
n represents the direction of travel of the phase fronts and p the apparent slowness
(inverse of wave velocity) in that direction. On substituting this plane wave form
into (5.4.1) we obtain
ω2p2cijklnjnlUk = ’ω2ρUi, (5.4.3)
which constitutes an eigenvalue problem for the slowness p for waves travelling in
the direction n
[p2cijklnjnl ’ ρδik]Uk = 0. (5.4.4)
Equation (5.4.4) has three roots for p2 for each direction n, associated with
polarisations speci¬ed by the eigenvectors U(r)(n). The slownesses pr are
independent of frequency.
In a general anisotropic medium, the slownesses vary with direction and can give
quite complex slowness and wave surfaces as, e.g., those illustrated in Figure 5.2.

5.4.1 Isotropic media
For an isotropic medium, the eigenvalue equation for slowness takes the form
[p2(» + μ)nink + p2μδik ’ ρδik]Uk = 0, (5.4.5)
where we have used the isotropic form for the elastic modulus tensor cijkl.
(i) P waves
Consider a displacement ¬eld oriented along the propagation direction, i.e.,
uP = cn, (5.4.6)
84 Linearised Elasticity and Viscoelasticity

3 3




Figure 5.2. Wave surfaces for a transversely isotropic (hexagonal) medium formed from
the envelope of plane wavefronts for a unit time of propagation.

then the slowness equation is
[(» + 2μ)p2 ’ ρ]ni = 0, (5.4.7)
so that we have a wave disturbance with slowness
a2 = ρ/(» + 2μ), (5.4.8)
and associated phase velocity
± = [(» + 2μ)/ρ]1/2. (5.4.9)
This longitudinal wave solution is called the P wave
uP = APn exp[iω(an · x ’ t)]. (5.4.10)



Figure 5.3. Elastic waves in a uniform isotropic medium: P waves with longitudinal mo-
tion and S waves with transverse (shear) motion.
5.4 Elastic waves 85

(ii) S waves
There is an alternative type of elastic wave motion in which the displacement is
transverse to the direction of motion. Consider a displacement ¬eld
uS = cs, s · n = 0,
with (5.4.11)
for which the slowness equation reduces to
[μp2 ’ ρ]si = 0, (5.4.12)
so that we have a wave disturbance with slowness
b2 = ρ/μ (5.4.13)
and associated phase velocity
β = [μ/ρ]1/2. (5.4.14)
This form of solution holds for any direction orthogonal to the direction of motion,
i.e., we have a degenerate eigenvalue problem from which we can choose two
orthogonal S wave vectors. It is conventional to choose one vector in the vertical
plane (denoted SV) and the other purely horizontal (denoted SH). This choice
simpli¬es the analysis of wave propagation in horizontally strati¬ed media, and
so we represent the S wave ¬eld as
uS = {BvsV + BHsH} exp[iω(bn · x ’ t)], (5.4.15)
where sV lies in the vertical plane through n and sH in the horizontal plane
[sV · n = sH · n = 0].

5.4.2 Green™s tensor for isotropic media
The wave disturbance at an observation point r produced by a delta-function point
force in the direction d at the origin in an unbounded isotropic medium can be
represented as
δ(t ’ r/±) δ(t ’ r/β)
4πu(r) = (d · ^)^ + ^ — (d — ^)
rr r r
r±2 rβ2
+(3(d · ^)^ ’ d) 3 [H(t ’ r/±) ’ H(t ’ r/β)],
rr (5.4.16)
where ^ is a unit vector in the direction r and H(t) is the Heaviside step function.
This result was ¬rst derived by Stokes in 1851 and can be deduced by many
different techniques, e.g., Fourier transform methods or the superposition of
potential solutions.
The disturbance consists of two delta functions in time, one longitudinal
travelling with the faster P-wave velocity ± and the other with transverse
displacement and the slower shear wave velocity β (this can be represented in
terms of SV and SH components if required). Between these two wavefronts there
86 Linearised Elasticity and Viscoelasticity


r /± r /β

Figure 5.4. The displacement from a point force in an unbounded medium with distinct P
and S arrivals linked by a ˜near-¬eld™ ramp disturbance.

is a ˜near-¬eld™ disturbance with intermediate character which decays more rapidly
with distance r (Figure 5.4).
The speci¬c form of the Green™s tensor Gk(x, t; 0, 0) for an observation point
with direction cosines γi is thus
Gk(x, t) = (3γiγj ’ δij) 3 [H(t ’ r/±) ’ H(t ’ r/β)]
1 1
+γiγj 2 δ(t ’ r/±) + (γiγj ’ δij) 2 δ(t ’ r/β). (5.4.17)
r± rβ
The radiation pattern for the far-¬eld terms which have a r’1 decay with distance
r away from the source takes a fairly simple form: the P waves depends on the
cosine of the inclination to the force direction and the S wave pattern on the sine of
the inclination (Figure 5.5).
More realistic sources can be obtained by using couples and dipoles whose
associated displacements can be derived by differentiation of the Green™s tensor.
The sense of differentiation is that of the displacement in the couple or dipole. The
double couple in Figure 5.5(c) is the equivalent force system for an in¬nitesimal
shear displacement on a surface aligned with one of the arms of the couple, and the
resulting radiation patterns are used in the interpretation of faulting mechanisms
from seismological observations.

5.4.3 Interfaces
The solution presented in the previous section is for elastic waves in an unbounded
medium, and more complicated behaviour results once material interfaces or a free
surface is introduced.
At a free surface, the requirement is that the traction should vanish. The surface
acts as a mirror for horizontally polarised S waves (SH); an incident plane SH wave
is re¬‚ected without change of amplitude and the SH component of the Green™s
5.4 Elastic waves 87



1 1 1
2 2 2

3 3 3


1 1 1
2 2 2

3 3 3


1 1 1
2 2 2

3 3 3

Figure 5.5. Far-¬eld radiation patterns for P and S waves from (a) a single point force, (b)
a vertical dipole, (c) a double couple.

function within the half-space can be found by introducing a ¬ctitious mirrored
source lying above the free surface. The situation is more complex for vertically
polarised waves. An incident plane P or SV wave will be re¬‚ected from the free
surface with an equal inclination to the normal to the surface, in addition there will
be conversion between P and SV waves with the inclination of the converted wave
determined by Snell™s law
sin i sin j
= , (5.4.18)
±0 β0
where i is the inclination to the normal of the P wave with surface wavespeed ±0,
and j is the inclination of the SV wave with surface wavespeed β0. In addition
there is a special form of wave that satis¬es the vanishing traction condition called
a Rayleigh wave. This consists of coupled exponentially decaying P and S waves
travelling with a horizontal velocity about 0.9 of the shear wave velocity. A surface
source in a half-space excites the Rayleigh wave in addition to the P and SV body
88 Linearised Elasticity and Viscoelasticity

waves; an elegant treatment of the development of the Green™s functions for this
case is given by Hudson (1981).
At an interface between two different media the boundary conditions are that
(i) the displacement u is continuous,
(ii) the traction „(n) associated with the normal n to the interface is continuous.
These boundary conditions can only be satis¬ed by the coupling of P and S waves
at the interface. For a plane wave incident on a horizontal interface the P and SV
wave ¬elds are coupled, but SH is independent. Such incident plane waves are
both re¬‚ected and transmitted through the interface with inclinations to the vertical
dictated by Snell™s law, which may be viewed as a condition for the continuity of
phase along the interface.
A major application of linearised elastic wave theory is in seismology, and
a comprehensive development of seismic wave theory with applications to the
interpretation of observed seismograms can be found in Kennett (2001, 2002), see
also Chapter 11.

5.5 Linear viscoelasticity
As in (4.6.1) we assume that the stress tensor σ depends linearly on the history of
strain et so that
σij(x, t) = Cijpq{et }, (5.5.1)

with ij, pq symmetry associated with the symmetries of the stress and strain tensor.

Creep and relaxation
Consider imposing a step-function stress cycle on a wire (i.e., a one-dimensional
σ(t) = σ0H(t), (5.5.2)
then the resulting strain for t > 0 can be written as
e(t) = [1 + ψ(t)], ψ(0) = 0, (5.5.3)
where ED is the dynamic modulus and σ0/ED is an instantaneous elastic response.
The monotonically increasing function ψ(t) is known as the creep function. If ψ(t)
approaches a constant as t ’ ∞, then we have long-term solid behaviour, whereas
for materials like pitch whose ultimate behaviour is that of a (viscous) ¬‚uid ψ(t)
increases without limit.
We can extend this concept to a more general stress history, by considering the
superposition of many small steps and then taking the limit as the step interval tends
to zero, to obtain
e(t) = σ(t) + ds σ(s)ψ(t ’ s) ,
™ (5.5.4)
5.5 Linear viscoelasticity 89

linear viscoelastic liquid
E 0 e(t)
σ0 viscoelastic solid

elastic solid

Figure 5.6. Strain evolution as a function of time.

a convolution integral over stress rate. If the creep function is differentiable for
t ≥ 0, we may use the properties of the convolution to rewrite the strain as

e(t) = σ(t) + ds σ(s)ψ(t ’ s) . (5.5.5)

An alternative viewpoint is to consider a strain cycle for which

e(t) = e0H(t), (5.5.6)

with associated stress for t > 0

σ(t) = EDe0[1 ’ φ(t)], φ(0) = 0. (5.5.7)

φ(t) is termed the relaxation function and is a monotonically increasing function
with a monotonically decreasing slope. If the long-term behaviour φ(∞) = 0 we
have a liquid.


σ (t)

partially relaxing (solid)

completely relaxing (liquid)


Figure 5.7. Stress evolution as a function of time.
90 Linearised Elasticity and Viscoelasticity

For an arbitrary differentiable strain history we can represent the stress behaviour
σ(t) = ED e(t) ’ ds e(s)φ(t ’ s) .
™ (5.5.8)

Once again, if the relaxation function is differentiable for t ≥ 0, we can reorganise
the convolution to give

σ(t) = ED e(t) ’ ds e(s)φ(t ’ s) . (5.5.9)

Since (5.5.4), (5.5.8) have to represent the same mechanical behaviour we can
establish the relation between the creep and relaxation functions. Consider the
application of a constant unit stress, then
e(t) = E’1[1 + ψ(t)], (5.5.10)

and so from (5.5.8) for t > 0

1 = 1 + ψ(t) ’ dsψ(s)φ(t ’ s). (5.5.11)
Simplifying (5.5.11), we require

ψ(t) = ds ψ(s)φ(t ’ s), t > 0. (5.5.12)
Similarly, we ¬nd

φ(t) = ds φ(s)ψ(t ’ s), t > 0. (5.5.13)
These two coupled integral equations have their simplest solution in the Laplace
transform domain. We de¬ne, e.g.,

¯ dt ψ(t)e’pt,
ψ(p) = (5.5.14)
and then (5.5.5) and (5.5.9) can be written as
e(p) = E’1[1 + pψ(p)]σ(p),
¯ ¯ (5.5.15)
σ(p) = ED[1 ’ pφ(p)]e(p).
¯ ¯ (5.5.16)
Since these relations are equivalent we require
¯ ¯
[1 + pψ(p)][1 ’ pφ(p)] = 1, (5.5.17)
which connects the Laplace transforms of ψ(t), φ(t). We can now de¬ne a
transform modulus
¯ ¯
E(p) = ED[1 ’ pφ(p)] = ED[1 + pψ(p)]’1, (5.5.18)
5.6 Viscoelastic behaviour 91

in terms of which the one-dimensional constitutive relation for linear viscoelasticity
can be written as
σ(p) = E(p)e(p).
¯ ¯ (5.5.19)
The formal analogy between the Laplace transform domain result (5.5.19) and the
equivalent elastic relation provides the formal basis for solving many viscoelastic
problems via the correspondence principle.

5.6 Viscoelastic behaviour
The relaxation function for the standard linear solid model introduced in Section
4.6 is
φ(t) = [1 ’ („R/„C)][1 ’ e’t/„R ], (5.6.1)
and the corresponding creep function is
ψ(t) = [(„C/„R) ’ 1][1 ’ e’t/„C ]. (5.6.2)
Although the standard linear solid model goes a long way towards providing a
description of a viscoelastic material it does not provide an adequate description
of any real material. More satisfactory representations are provided by joining
standard linear solid elements in series or parallel.
In the series case, the stress in each element is the same whilst the strain is the
sum of the strains in each element. The creep function is then a sum of terms of the
form (5.6.2)
f(„n)[1 ’ e’t/„n ],
ψ(t) = (5.6.3)

where f(„n) may be regarded as the contribution of a creep process with a
characteristic time „n.
For a set of parallel elements, the stress is cumulative and the relaxation function
can be built from a set of terms of the form (5.6.1)
g(„n)[1 ’ e’t/„n ].
φ(t) = (5.6.4)

In the most general case we can consider a spectrum of relaxation and creep
The creep function for the Burgers material takes the form
ψ(t) = ψU + ”ψ[1 ’ e’t/„] + , (5.6.5)
where ψU is the reciprocal of the relevant instantaneous modulus, ”ψ describes the
extent of viscoelastic relaxation with a relaxation time „C and ·M is the viscosity
of the Newtonian viscous element. The equivalent Maxwell relaxation time for
the in¬‚uence of viscosity is „M = ·MψU. Once again more general behaviour
92 Linearised Elasticity and Viscoelasticity

can be produced by a sum or spectrum of creep and relaxation processes, and such
a combination of Burgers elements can be used as a versatile representation of
deformation in the mantle of the Earth.

Isotropic linear viscoelasticity
For an isotropic medium the tensor function Cijpq in (5.1.5) can be represented as
an isotropic tensor of the form
Cijpq = Lδijδpq + M(δipδjq + δiqδjp), (5.6.6)
where L and M are scalar functions, and so
σij(x, t) = δijL{et } + 2M{et }, (5.6.7)
kk ij

and we can separate the dependence of the dilatation component from the deviatoric
stress and strain.
By analogy with the equations of isotropic elasticity we can write the
stress“strain relation as
σij(x, t) = »0δijekk(x, t) + 2μ0eij(x, t)
t t
™ ™
+ ds δijekk(x, s)R»(t ’ s) + ds 2eij(x, s)Rμ(t ’ s), (5.6.8)
0 0
in terms of relaxation functions R», Rμ and a representation analogous to (5.5.9)
for the one-dimensional case.
An equivalent representation to (5.6.8) expressing the strain as a function of the
stress history can be made in terms of creep functions C», Cμ. These are related
to the relaxation functions via coupled integral equations similar to (5.5.13). The
correspondence principle in the Laplace transform domain can be extended to these
three-dimensional problems and provides a means of solving many viscoelastic

5.7 Damping of harmonic oscillations
We consider a harmonic strain cycle of the form
eij = Eije’iωt, (5.7.1)
and then the corresponding stress is
σijeiωt = »0δijEkk + 2μ0Eij
™ ™
ds [δijEkkR»(t ’ s) + 2EijRμ(t ’ s)]eiω(t’s),
+ (5.7.2)
where t = ’T is the time at which the disturbance starts. Now let T ’ ∞ and
introduce the Fourier transforms of the time derivatives of the relaxation functions
∞ ∞
™ ™
dt Rμ(t)eiωt.
»1(ω) = dt R»(t)e , μ1(ω) = (5.7.3)
0 0
5.7 Damping of harmonic oscillations 93

We can express the stress as
σij = {[»0 + »1(ω)]δijEkk + 2[μ0 + μ1(ω)]Eij}e’iωt, (5.7.4)
with a stress“strain relation in the same form as for an elastic medium but with
frequency-dependent complex moduli
» = »0 + »1(ω), μ = μ0 + μ1(ω).
¯ (5.7.5)
We may therefore employ the descriptions of P and S waves developed in section
5.4, but have to allow for the frequency dependent velocities and consequent energy
The rate of energy dissipation per unit volume due to deformation of the medium
is from (4.2.11)
ρ U = σij ≈ σijeij,
™ (5.7.6)
Dt ‚xi
in the linear approximation. The average work dissipated per unit time and volume
is thus
¯^ ^ ij
¯ — iωt
= ’ 1 ω Im[κ””— + 2μEij E— ],
W= 2 Re[iωσijEije ] ¯ (5.7.7)
where the dilatation ” = Ekk, and the deviatoric strain Eij = Eij ’ 1 ”δij. The
complex bulk modulus
¯ ¯ 3¯
κ = » + 2 μ = κ0 + κ1(ω). (5.7.8)
A convenient measure of the rate of energy dissipation is the loss factor Q’1(ω)
which may be de¬ned as
Q’1 = , (5.7.9)
where ”W(ω) is the energy loss in a cycle at frequency ω and W0(ω) is the
˜elastic™ energy stored in the oscillation, i.e., the sum of the strain and kinetic
energy associated with just the instantaneous elastic moduli. For purely deviatoric
^ ^ ij
’ 1 Im μ(ω)EijE—
¯ ’Im μ1(ω)
Qμ (ω) = = . (5.7.10)
1 ^ ^— μ0
Similarly for purely dilatational disturbances
’Im κ1(ω)
Q’1(ω) = ; (5.7.11)
normally loss in pure dilatation is much less signi¬cant than loss in shear and
Q’1 Q’1.
κ μ
In general the loss factor for dilatation or shear can be represented as
’Im m1(ω)
Q’1(ω) = , (5.7.12)
94 Linearised Elasticity and Viscoelasticity

where m is the appropriate complex modulus. Thus, for the standard linear solid
„’1 ’ iω
m(ω) = m2 (5.7.13)
„’1 ’ iω
ω2 + „’1„’1 ’ iω(„’1 ’ „’1)
= m2 , (5.7.14)
ω2 + „’2
and the instantaneous modulus is m2, and thus the loss factor
„’1 ’ „’1
Q’1(ω) R C
=ω . (5.7.15)
ω2 + „’2
By superposing a spectrum of relaxation times, a wide range of frequency
behaviour can be simulated.
The dissipation via Q’1(ω) is accompanied by a frequency dependent
correction to the real part of the complex modulus and hence to the apparent
ω2 + „’1„’1
Re m(ω) = m2 . (5.7.16)
ω2 + „’2
This property is a consequence of a causal dissipation mechanism. Since the

relaxation contributions depend only on the past history of the strain Rm(t)
vanishes for t < 0, so that the transform m1(ω) must be analytic in the upper
half-plane. In consequence the real and imaginary parts of m1(ω) are Hilbert
transforms of each other,
∞ ∞
ω Q’1(ω )
1 Im m1(ω ) 2m0
Re m1(ω) = P =’ ,
dω P dω
ω 2 ’ ω2
π ω ’ω π
’∞ ’∞
where P denotes the Cauchy principal value of the integral.
Continua under Pressure

The linearised development of elasticity in the previous chapter can be applied
about a stress-free state and has therefore direct application to the shallower parts
of the Earth. As we shall see, we can make an incremental treatment about
the stress-state at depth, which will be dominantly hydrostatic, and derive the
incremental elastic moduli for small disturbances such as those created by the
passage of seismic waves.
It is dif¬cult in experimental studies to achieve the pressures and temperatures
appropriate to a direct study of the Earth. For the lower mantle, in particular,
recourse needs to be made to extrapolation from available conditions. A convenient
summary of the property of Earth materials at depth is provided by the equation of
state, which is commonly expressed as a relation between the pressure, the speci¬c
volume and the temperature. This simpli¬ed approach ignores non-hydrostatic
stresses and any non-isotropic response to the assumed hydrostatic pressure. Finite
strain representations are used in the Birch“Murnaghan equations to relate the
conditions at depths to the properties at zero pressure.
Long-term response to large non-hydrostatic stresses leads to irreversible
(plastic) deformation, and so over geological time materials are expected to have
yielded to such stresses, leaving just the hydrostatic load at depth. Non-hydrostatic
stresses are, however, important at shallow depth where they approach the
magnitude of the pressure and are responsible for the occurrence of earthquakes.
Although materials within the Earth are undoubtedly anisotropic on a
microscopic scale, the properties of the polycrystalline aggregates averaged over a
typical seismic wavelength are close to isotropic except towards the top and bottom
of the mantle.

6.1 Effect of radial strati¬cation
A convenient simpli¬cation when considering the properties of the Earth at depth
is to concentrate on the variations as a function of radius and to regard any
three-dimensional variation as a perturbation around this state. This enables the
use of a simple equation of state in terms of the hydrostatic pressure as a function

96 Continua under Pressure

of radius, which depends directly on knowledge of the density distribution within
the Earth.

6.1.1 Hydrostatic pressure
The predominant stress-¬eld within the Earth is imposed by the hydrostatic
pressure induced by gravitational effects. The gravitational force is directed
radially inwards and so g(r) = ’g^r, with the acceleration due to gravity g given
ds s2ρ(s),
g= 2 (6.1.1)
r 0

in terms of the radial distribution of density ρ(r). For hydrostatic balance, the
pressure is determined by
+ ρg = 0, i.e., p(r) = ds ρ(s)g(s), (6.1.2)
‚r r

with the integration taken from the free surface at re , where p(re ) = 0, down to
radius r.

16 400


12 300

[Mg/m3], [m/s2]


8 200

4 100

0 0
0 1000 2000 3000 4000 5000 6000
Depth [km]

Figure 6.1. Pressure p, gravitational acceleration g for the density distribution associated
with the AK135M model.
6.1 Effect of radial strati¬cation 97

The behaviour of the pressure and the gravitational acceleration within the Earth
is illustrated in Figure 6.1, using the density distribution ρ for the AK135 model
(Figure 1.6) modi¬ed to remove the low density zone in the upper mantle. Both the
pressure and the gravitational acceleration are continuous across the discontinuities
in density associated with the phase transitions in the upper mantle, the core“mantle
boundary and the boundary between the inner and outer core.
The gravitational acceleration grows through the core as we move from the centre
of the Earth outwards and reaches its maximum value at the core“mantle boundary,
because the mantle material is lighter and the increased volume associated with
spherical shells closer to the surface of the Earth does not compensate for the
lowered density. The secondary maximum at the base of the transition zone again
re¬‚ects the step down to lower densities in the outer layers.
The pressure increases steadily with depth with slight changes in slope at each of
the discontinuities in density. The pressure tends to a constant value at the centre of
the Earth of about 370 GPa; the reduction in slope for the inner core is associated
with the steady decrease of g towards zero as the centre of the Earth is approached.

6.1.2 Thermodynamic relations
Most equations of state can be regarded as statements of the conservation of energy.
In terms of the total energy U of a system, the ¬rst law of thermodynamics can be
written as
dU(S, V) = T dS ’ p dV, (6.1.3)
as a function of temperature T , entropy S, pressure p and volume V. Because
we want to consider the thermal state of a material it is most convenient from a
theoretical viewpoint to work in terms of temperature T and volume V, since these
can be related directly to the reference state, whereas stresses are most naturally
introduced in the deformed state (Chapter 3). The ¬rst law can then be written in
terms of the Helmholtz free energy F = U ’ T S:
dF(T, V) = ’S dT ’ p dV. (6.1.4)
In experiments on earth materials, it is usually the pressure rather than the volume
which is adjusted, and we have noted in section 4.2 that propagation is adiabatic
rather than isothermal. From (6.1.4)
p= . (6.1.5)
‚V T
To provide a key to the mineral physics literature we will use the notation K rather
than κ for the bulk modulus and G for the shear modulus instead of μ. The bulk
modulus (incompressibility) K is de¬ned as
‚p ‚p
KT = ’V , KS = ’V . (6.1.6)
‚V ‚V
98 Continua under Pressure

KT applies to isothermal conditions in slow laboratory experiments and KS to the
adiabatic conditions in transient wave propagation.
A related quantity is the ˜seismic parameter™
KS ‚p
¦= = . (6.1.7)
ρ ‚ρ S

For an isotropic material, the bulk-sound speed

(±2 ’ 4 β2),
φ= ¦= (6.1.8)

where ± is the P wavespeed (5.4.9) (also designated vp) and β is the S wavespeed
(5.4.14) (also designated vs).
Regarding pressure p as a function of V, T , S
‚p ‚p ‚p ‚T
= ’ , (6.1.9)
‚V ‚V ‚T ‚V

and so the difference between the isothermal and adiabatic moduli is
‚p ‚T
K T ’ KS = V . (6.1.10)
‚T ‚V

The volume expansion coef¬cient ±th
1 ‚V ‚ ln V
±th = = ; (6.1.11)
V ‚T ‚T
p p

and so we can use the chain rule for partial derivatives to write
‚p 1 ‚V ‚p
=’ V = ±th KT . (6.1.12)
‚T V ‚T ‚V
V p T

The adiabatic derivative of temperature with respect to speci¬c volume can be
written as
‚T T
= ’γth , (6.1.13)
‚V V

where γth is the thermodynamic Gr¨ neisen parameter. The speci¬c heats at
constant volume CV and constant pressure Cp are given by
‚S ‚S
CV = T , Cp = T , (6.1.14)
‚T ‚T
V p

where S is the speci¬c entropy.
We are therefore able to express the Gr¨ neisen parameter in the form
±th VKS ±th VKT
γth = = . (6.1.15)
6.1 Effect of radial strati¬cation 99

Note also that
‚T T
= γth . (6.1.16)
‚p KS

In terms of these new variables we can express the difference between the
isothermal and adiabatic moduli as
KT ’ KS = ’±th KT Tγth , (6.1.17)
and so we ¬nd that
= 1 + γth ±th T. (6.1.18)
The behaviour of the bulk and shear moduli in the Earth with depth, inferred
from the properties of seismic waves, is illustrated in Figure 6.2 for the AK135

16 16

14 14

12 12

10 10

[1011 Pa]
8 8
[1011 Pa]

6 6

4 4

G p
2 2

0 0
0 1000 2000 3000 4000 5000 6000
Depth [km]

Figure 6.2. The behaviour of the bulk modulus K and shear modulus G as a function of
depth for the AK135M model, compared to the variation in pressure p.

The effect of compression with depth plays a major role in determining the
physical properties, but the increase in temperature with depth also plays a
signi¬cant role. Seismic wave propagation will be close to isentropic conditions
and hence measurements are of the adiabatic modulus KS. For both the bulk of the
100 Continua under Pressure

lower mantle and the core, where the composition is expected to be homogeneous,
the bulk modulus K is approximately proportional to pressure p.

6.2 Finite strain deformation
As we have seen in Section 4.3, the elastic properties of a material can be directly
related to a local Lagrangian treatment about a suitable reference state. However,
until recently, the appropriate conditions for much of the Earth™s mantle were not
accessible to experimental studies. In particular, even where measurements were
carried out at high pressure, the temperature was close to ambient.
Estimates of the material behaviour in the Earth have therefore been made by
using a systematic extrapolation from the low-pressure regime of most experiments
to likely Earth conditions. A suitable equation of state for mineral systems can be
derived by examining the in¬‚uence of large ¬nite strain, suf¬cient to change the
volume of a unit cell from its zero-pressure value to that at depth. The aim of
such studies is to try to understand the distribution of physical properties revealed
from seismological work, as in Figure 6.2, in terms of the appropriate mineral
assemblages and temperature regimes.
We recall from Section 2.2.5 the de¬nition of the Green strain in the Lagrangian
‚ui ‚uj ‚uk ‚uk
E = 1 (FT F ’ I), Eij = 1
2 ‚ξ + ‚ξ + ‚ξ ‚ξ , (6.2.1)
j i i j

and the Cauchy strain in the Eulerian frame
‚ui ‚uj ‚uk ‚uk
e = 1 (I ’ (FFT )’1), 1
eij = + ’ . (6.2.2)
2 2 ‚xj ‚xi ‚xi ‚xj
The volume transformation under ¬nite strain can be expressed either in terms
of the Green strain in a Lagrangian viewpoint,
= J2 = det[FT F] = det[I + 2E], (6.2.3)
or alternatively in terms of the Cauchy strain in an Eulerian perspective,
= J’2 = det[FT F]’1 = det[2e ’ I]. (6.2.4)
Consider hydrostatic pressure on a cubic or isotropic material. Then in the
Lagrangian formulation the Green strain tensor takes the form
det[I + 2E] = (1 + 2·)3,
Eij = ·δij, (6.2.5)
and so, from (6.2.3),
2/3 2/3
V ρ0
1 1
·= ’1 = ’1 ; (6.2.6)
2 2
V0 ρ
6.2 Finite strain deformation 101

this scalar · varies from 0 to ’ 1 as the pressure varies from 0 to ∞.
Alternatively, in the Eulerian formulation, the Cauchy strain is given by
det[2e ’ I] = (2 ’ 1)3,
eij = δij, (6.2.7)
and so
2/3 2/3
V0 ρ
’1 ’1
= ’1 = ’1 . (6.2.8)
2 2
V ρ0

In terms of the strain variables · and , the volume change
V0 ρ
= (1 + 2·)’3/2 = (1 ’ 2 )3/2.
= (6.2.9)
V ρ0
In an elastic material the strain energy W can be equated with the
speci¬c Helmholtz free energy F/ρ (we neglect thermal effects in either slow
near-isothermal deformations or the nearly isentropic incremental deformation
associated with the passage of seismic waves). The Helmholtz free energy can
then be expressed as a function of strain, e.g., through its dependence on the Green
strain E. With a generalisation of (6.1.4) for speci¬c quantities
ρ dW = ’ρS dT + ρ0σij dAij, (6.2.10)
in terms of the displacement gradient A = F ’ I (5.1.1). The Cauchy stress tensor
σij is thus
ρ ‚F
σij = ; (6.2.11)
ρ0 ‚Aij
which is equivalent to our previous expression (4.3.10)
det F σij = Fim Fjn, (6.2.12)
The effective elastic moduli Cijkl are de¬ned by
ρ ‚2W
Cijkl = = ’ σijδkl, (6.2.13)
‚Akl ρ0 ‚Aij‚Akl
since ‚ρ/‚Akl = ’ρδkl. The differentiation may be carried out either isothermally
or adiabatically to give the appropriate moduli. In terms of the dependence on
strain, the effective moduli in the case of hydrostatic pre-stress take the form
FjnFlq ’ pδkl,
det F Cijkl = FimFkp (6.2.14)
δkl = δikδjl + δilδjk ’ δijδkl. (6.2.15)
102 Continua under Pressure

6.3 Expansion of Helmholtz free energy and equations of state
A commonly used approach to handling ¬nite strain effects is to expand the
Helmholtz free energy F as a Taylor series in the elastic strain from the
zero-pressure state. This yields a power series expansion in terms of a suitable
measure of ¬nite strain.
The level of strain required to attain the required density at depth in the Earth
is then employed to determine the pressure and the other desired properties of the
medium in the strained, i.e. pressured, state. For the Lagrangian viewpoint we
make a power-series expansion in the strain · from (6.2.6),

F = V0 (6.3.1)
Whereas for the Eulerian viewpoint (Birch“Murnaghan) the expansion is in terms
of the strain from (6.2.8),

Bi i.
F = V0 (6.3.2)
The isothermal thermodynamic relations
‚F ‚p ‚K
p=’ , KT = ’V , KT = (6.3.3)
‚V ‚V ‚T
enable us to identify the coef¬cients in these series expansions.
For the Lagrangian expansion the coef¬cients in the expansion (6.3.1) for the
Helmholtz free energy F = V0 ∞ Ai·i are

A0 = A1 = 0, (6.3.4)
A2 = 9 K0, (6.3.5)
A3 = ’ 9 K0K0, (6.3.6)

A4 = 9 K0[K0K0 + K0(K0 + 1) ’ 1 ], (6.3.7)
2 9
where the subscript 0 indicates evaluation at zero pressure and ambient temperature
T0. The Lagrangian equation of state is then
’1/3 ∞
‚F ‚· V
(i + 1)Ai+1·i,
= ’1
pL = ’ (6.3.8)
‚· ‚V V0
‚· V
= 1 (1 + 2·)’1/2 = 1
. (6.3.9)
3 3
‚V V0
Thus the pressure pL is represented as
· ’ 3 K0·2 + · · · .
pL = ’3K0 (6.3.10)
6.3 Expansion of Helmholtz free energy and equations of state 103

and the isothermal modulus can be found in terms of · from (6.1.6).
For the Eulerian viewpoint the coef¬cients in the expansion (6.3.2) for the
Helmholtz free energy F = V0 ∞ Bi i can be expressed as

B0 = B1 = 0, (6.3.11)
B2 = A2 = 9 K0, (6.3.12)
B3 = ’ 9 K0(K0 ’ 4), (6.3.13)

B4 = 9 K0[K0K0 + K0(K0 ’ 7) ’ 143
9 ], (6.3.14)

and the equation of state for the pressure takes the form
V 2
’ 3 (K0 ’ 4)
pE = ’3K0 + ··· . (6.3.15)
When we reinstate the full dependence on by expressing the volume change in
terms of the strain, we ¬nd

pE = ’(1 ’ 2 )5/2 C1 ’ C2 2/2 + C3 3/6 , (6.3.16)

C1 = 2 B2 = 3K0, C2 = ’2B3, C3 = 6B4. (6.3.17)

The corresponding expression for the isothermal bulk modulus is

KT = 1 (1 ’ 2 )5/2 C1 + (C2 ’ 7C1) + (C3 ’ 9C2) 2 3
’ 11C3 (6.3.18)

= (1 ’ 2 )5/2 K0 + K1 + K2 2
+ ··· . (6.3.19)

For moderate changes in volume 1.0 > V/V0 > 0.9 the predictions of the
Lagrangian and Eulerian equations of state are essentially equivalent, but for larger
strains the pressure estimate pE increases more rapidly than pL.
Following Birch (1952) the preferred formulation for studies of the Earth is that
based on the Eulerian viewpoint because the in¬‚uence of the third-order term in
strain is much reduced. We note that the third-order coef¬cient B3 depends on the
pressure derivative of the bulk modulus through K0 ’ 4. For most minerals the
pressure derivative K0 lies in the range 3 to 6 and is frequently close to 4 and so
B3, C2 are small, with only a limited in¬‚uence from uncertainties in measurement.
The third order coef¬cient in the equivalent Lagrangian expression depends solely
on K0, and has a stronger in¬‚uence.
The third-order truncation of the Birch“Murnaghan equation of state has the
advantage that no knowledge is required of second pressure derivatives of physical
properties that are more dif¬cult to measure than the ¬rst derivatives, and hence are
rather poorly known. The fourth-order representation is nevertheless more ¬‚exible
because of the larger number of parameters.
104 Continua under Pressure

The isothermal seismic parameter
¦T = (6.3.20)
(1 ’ 2 ) C1(1 ’ 7 ) + C2( ’ 9 2/2) + C3( 2
’ 11 3/6) ;
the grouping of the terms emphasises that the coef¬cients of the powers of are
not all independent.
With the aid of an Eulerian representation comparable to (6.2.14) we can develop
representations for the wavespeeds of both P and S waves in the form
±2 = v2 = 2
(1 ’ 2 ) L0 + L1 + L2 + ··· , (6.3.21)
β2 = v2 = 2
(1 ’ 2 ) G0 + G1 + G2 + ··· , (6.3.22)
L0 = K0 + 4 G0, L1 = K1 + 4 G1, etc. (6.3.23)
3 3

The ¬rst and second pressure derivatives of the shear and bulk moduli at zero
pressure are then given by
G0 = 1 (5G1 ’ G2)/K0,

K0 = 1 (5L1 ’ L2)/K0 ’ 4 G0,
3 3
G0 = 9 [5G3 + 3(K0 ’ 4)G2 ’ 5(3K0 ’ 0

K0 = 1 [5L3 + 3(K0 ’ 4)L2 ’ 5(3K0 ’ 5)L1]/K2 ’ 4 G0.
9 3
The coef¬cients for the contributions from the bulk modulus can be related to
quantities which can be measured, or estimated, from mineral physics experiments.
However, the corresponding results for the shear modulus are rather poorly
known, it is only recently that measurements of shear properties under conditions
appropriate to the Earth™s transition zone have been accomplished and the
derivatives are weakly constrained. Ab initio quantum mechanical calculations of
mechanical properties are also easier to perform for the bulk modulus, because
allowance does not need to made for anelastic thermal effects (see Chapter 9).
The ¬nite strain theory does not, by itself, provide a relation between the bulk
and shear moduli. Current seismological models are consistent with a dependence
of the shear modulus G on pressure p and bulk modulus K of the form
G = aK ’ bp. (6.3.25)
This behaviour is ilustrated in Figure 6.3 for the AK135 model, by the linear
segments associated with different parts of the mantle. The slopes are quite
similar, except for the zone just below the 660 km discontinuity and the highly
heterogeneous D region at the base of the mantle.
6.4 Incremental stress and strain 105








. 4
( 16)