. 6
( 16)


We also require from the continuity condition that
w = 0. (7.6.19)
The three components of the momentum equation are
1 ‚p
’2 ©v = ’ +ν 2,
ρ ‚x ‚z
1 ‚p
2 ©u = ’ + ν 2,
ρ ‚y ‚z
0= , (7.6.20)
and hence
‚p ‚p0 ‚p ‚p0
= = 0, = = 0. (7.6.21)
‚x ‚x ‚y ‚y
The equations to be solved for the ¬‚ow are
’2©v = ν 2 ,
’2©(u0 ’ u) = ν 2 , (7.6.22)
with boundary conditions
u = v = 0, z = 0,
u ’ u0, v ’ 0, z ’ ∞. (7.6.23)
The solution to this equation system is
u = u0 [1 ’ e’z/” cos(z/”)],
v = u0 e’z/” sin(z/”), (7.6.24)
” = (ν/©)1/2. (7.6.25)
We notice that the distance over which the solution approaches the interior solution
is of order ”. In other words, ” gives the thickness of the boundary layer; this layer
is called the Ekman layer.
Comparing ” with the expression for the Ekman number Ek con¬rms that the
boundary layer thickness corresponds to a local Ekman number of unity.
Continuum Equations and Boundary

We summarise here the main equations for continuum behaviour that need to
be employed in a broad range of applications, by drawing on the development
in earlier chapters. We start with the conservation laws for mass, momentum
and energy expressed in differential form, that need to be supplemented by the
appropriate constitutive laws to express the rheological state of the medium. We
then consider the boundary conditions that prevail at the surface of a continuum or
an interface between two different continua. The normal components of velocity
and the stress tensor are required to be continuous across a general interface.
The tangential components of velocity are also continuous for solid“solid and
¬‚uid“¬‚uid boundaries and also for a solid“viscous ¬‚uid interface under a “no-slip”
condition. Heat ¬‚ux is continuous across an interface, but because of variations
in thermal conductivity temperature gradients can have a jump. At a phase
boundary, additional thermodynamic constraints need to be applied to describe the
equilibrium scenario along the interface.
Hitherto, we have concentrated on the way in which a continuum responds
to deformation or imposed stress, but we also need to take into account
electromagnetic phenomena. The iron-rich core of the Earth is a conducting ¬‚uid
where a complex interaction of motion and electromagnetic effects leads to dynamo
action and creates the Earth™s internal magnetic ¬eld. We therefore provide a
brief development of the topic of continuum electrodynamics and show how the
continuum equations need to be modi¬ed to accommodate magnetic effects.
In many situations we have to understand the interaction of the temperature and
¬‚ow ¬elds, and so provide an introduction to the solution of diffusion problems for
temperature and the corresponding heat ¬‚ow.

8.1 Conservation equations
In earlier chapters we have established a number of results based on conservation
properties, working from integral forms to derive the appropriate local equations.
These results based on the conservation of mass, momentum and energy are
brought together here and are recast into explicit differential representations of the
conservation equations.

132 Continuum Equations and Boundary Conditions

8.1.1 Conservation of mass
The velocity ¬eld v(x, t) is related to the density ¬eld in a moving continuum by the
continuity equation (2.5.6) that represents a local formulation of the conservation
of mass,
+ ∇ · (ρv) = 0; (8.1.1)
this relation is equivalent to the Lagrangian form (2.5.7),

+ ρ∇ · v = 0; (8.1.2)
so that incompressibility (Dρ/Dt ≡ 0) requires the divergence of the velocity ¬eld
to vanish (∇ · v = 0).

8.1.2 Conservation of momentum
The acceleration of a continuum produced by the various forces acting upon it is
given by (3.2.4) in conjunction with (2.4.5), in terms of the velocity v(x, t),
D ‚v
v= + (v · ∇x )v = ∇x σ + g, (8.1.3)
Dt ‚t
i.e., in component terms
D ‚ ‚
vj = + vk vj = σij + gj. (8.1.4)
Dt ‚t ‚xk ‚xi
The relation between the stress and deformation is determined by the relevant
constitutive equation.
For a viscous ¬‚uid, from Section 4.5, the stress tensor is given by
‚vi ‚vj 2 ‚vk ‚vk
σij = ’pδij + · + ’ δij + ζδij , (8.1.5)
‚xj ‚xi 3 ‚xk ‚xk
where we have allowed for the possibility of bulk viscosity (ζ).
For linearised isotropic elasticity, from Section 5.2, we ¬nd that
‚vi ‚vj
σij = (K ’ 2 G)δij + 2G + , (8.1.6)
3 ‚xk ‚xj ‚xi
in terms of the displacement u, the bulk modulus K and shear modulus G ≡ μ.
For linearised viscoelasticity we also need to include terms that depend on the
history of deformation.
We can write (8.1.3) in a form that more directly re¬‚ects the conservation of
momentum in the current frame,
‚ ‚
(ρvj) = ’ Πij + ρgj, (8.1.7)
‚t ‚xi
8.1 Conservation equations 133

where Πij is the momentum ¬‚ux density tensor
Πij = ρvivj ’ σij; (8.1.8)
as can be veri¬ed by direct differentiation.

8.1.3 Conservation of energy
In the discussion of energy balance in Section 4.2 we have established the relation
ρ U = h ’ ∇x q + σ.∇x v, (8.1.9)
where h is the rate of internal heat production and the heat ¬‚ux vector q is
determined by the temperature gradient
q = ’k∇x T, (8.1.10)
where k is the thermal conductivity. When we include the kinetic energy
contribution as in (4.2.9) we have
(U + 1 v2) = h ’ ∇x q + ∇x · [σ · v] + g · v,
ρ (8.1.11)
where g is the external force ¬eld (normally gravitational). We express g in terms
of a potential •
g = ’∇x •. (8.1.12)
When g is a function of position rather than time
g · v = ’v · ∇x • = ’ •, (8.1.13)
and then we can include the potential energy of the gravitational ¬eld • as an
additional energy contribution. For ¬‚ow in the mantle and core, self-gravitational
effects will be small so the time-independent approximation in (8.1.13) is well
The material derivative of the speci¬c energy
D 12 ‚
( 2 v + U + •) = ρ ( 1 v2 + U + •) + ρv · ∇x ( 1 v2 + U + •), (8.1.14)
‚t 2 2
and the partial derivative of the total energy
‚ ‚ ‚ρ 1 2
[ρ( 1 v2 + U + •)] = ρ ( 1 v2 + U + •) + ( v + U + •); (8.1.15)
‚t 2 ‚t 2 ‚t 2
further, from the continuity condition (2.5.8),
= ’∇x · (ρv). (8.1.16)
134 Continuum Equations and Boundary Conditions

We can combine the expressions (8.1.11) and (8.1.14)“(8.1.16) to produce a
representation of the rate of change of energy in the spatial domain
‚1 2
( ρv + ρU + ρ•) = h ’ ∇x Q, (8.1.17)
‚t 2
in terms of the energy ¬‚ux density
Q = ( 1 ρv2 + ρU)v ’ v · σ ’ k∇x T. (8.1.18)

Equation (8.1.17) is a direct differential representation of the conservation of

8.2 Interface conditions
The surface of the Earth represents a substantial density contrast between rock and
air, and for many purposes the presence of the air can be disregarded and the surface
treated as if it were a free surface in vacuum. Thus the surface traction „(nf) must
vanish, where nf is the normal to the free surface. If the surface can be regarded
as locally horizontal, the free surface condition translates to a requirement on the
3-components of the stress tensor
σ13 = σ23 = σ33 = 0. (8.2.1)
The corresponding condition for a static ¬‚uid would be the vanishing of pressure.
However, the density contrasts between the ocean and atmosphere are smaller than
for rock and so it is preferable to use the full ¬‚uid“¬‚uid boundary conditions.
For a stationary rigid boundary SB of a ¬‚uid, with normal nB we require
v · nB = 0. (8.2.2)
At an internal interface between two solids or ¬‚uids we require continuity of the
traction vector acting on the surface (3.2.11)
[„j]+ = ni[σij]+ = 0, (8.2.3)
’ ’

where n is the interface normal. For the isotropic pressure distribution in static
¬‚uids, (8.2.3) requires continuity of the pressure.
At welded contact between two solids or a ¬‚uid“¬‚uid interface, the velocity must
be continuous so that
[v]+ = 0. (8.2.4)

At a ¬‚uid“solid boundary, the normal component of the velocity will be
continuous, so we need
[v · n]+ = 0. (8.2.5)

For a viscous ¬‚uid in contact with a solid boundary, the no-slip boundary condition
requires the velocity to be continuous as in (8.2.4). For a boundary between
8.3 Continuum electrodynamics 135

an inviscid ¬‚uid and a solid slip can occur, and there is no requirement that the
tangential velocities need to be the same.
Where sediments with slow elastic wavespeeds lie close to the surface, coupling
between air and ground can take place and the ¬‚uid“solid boundary conditions need
to be used.
Across a material interface we require the heat ¬‚ux to be continuous and so
k(1)∇T (1) ’ k(2)∇T (2) · n = 0 (8.2.6)
so that there will, in general, be a discontinuity in the temperature gradient n · ∇T
across the boundary, since the thermal conductivity will be different on the two
At a phase transition in a material, the physical properties will normally
differ on the two sides of the phase boundary and so we need to impose the
appropriate continuity conditions as derived above on the deformation and force
¬elds. The different phases have to be in equilibrium at the phase interface and
as a result the position of the phase boundary will be a function of the conditions
to which it is subjected. Thermodynamic arguments require that there be no net
change in the Gibbs free energy across the phase transition. The slope of the
pressure“temperature (pb ’ Tb) relation for the mutual phase boundary is then
given by the Clausius“Clapeyron equation
= , (8.2.7)
Tb ”v
where L is the latent heat of the phase transition and ”v is the change in speci¬c
volume across the phase interface. The speci¬c entropy change is ”Sb = L/Tb.

8.3 Continuum electrodynamics
We have so far concentrated on the effects of deformation, but alongside
mechanical effects we have the possibility of electromagnetic in¬‚uences. The
electrical conductivity of the mantle modulates how external magnetic disturbances
interact with the Earth, and the character of the magnetic ¬eld of the geodynamo
as seen at the surface. If a conducting ¬‚uid moves in a magnetic ¬eld, electric
¬elds are induced with resulting ¬‚ow of electric current. The magnetic ¬eld exerts
forces on the electric currents that can alter the ¬‚uid ¬‚ow; the currents themselves
also alter the magnetic ¬eld. Such intertwining of ¬‚uid ¬‚ow and electromagnetic
phenomena within the Earth™s core creates dynamo action to produce a sustained
internal magnetic ¬eld.

8.3.1 Maxwell™s equations
The electromagnetic ¬eld can be described by the electric vector E and the
magnetic induction B, and the interaction of the ¬eld with matter requires the
136 Continuum Equations and Boundary Conditions

speci¬cation of a further set of vectors: the electric displacement D, the magnetic
vector H and the electric current density j. The space and time derivatives of the
¬ve vectors in the current state are related by
‚D ‚B
∇— H ’ = j, ∇— E + = 0,
‚t ‚t (8.3.1)
∇ · D = ρ, ∇ · B = 0,
where ρ is the density of electric charge. Charge conservation requires
∇·j+ = 0. (8.3.2)
In free space
D= B = μ0H,
0E, (8.3.3)
where 0 is the permittivity and μ0 the magnetic permeability of free space, such
that 0μ0 = 1/c2, where c is the speed of light.

8.3.2 Electromagnetic constitutive equations
The levels of electromagnetic ¬elds within the Earth are such that the constitutive
equations are linear relations between ¬elds. The frequencies of interest normally
lie well below those of crystal vibrations so that the dielectric constant and the
magnetic permeability μ are essentially constant and show no noticeable frequency
For a conductor, the current j is related to the electric vector E by
j = σE,
˜ (8.3.4)
where σ is the speci¬c conductivity; note that we have written σ rather than the
˜ ˜
more conventional σ to avoid con¬‚ict with the notation we have used for stress.
Inside a material the electric displacement D is related to the electric vector E
through the dielectric constant ,
D= E. (8.3.5)

In a similar way the magnetic induction B is derived from the magnetic vector H
by the action of the magnetic permeability μ,
B = μ0μH. (8.3.6)
Except for a few magnetic minerals, μ ≈ 1 for Earth materials. The temperature
coef¬cient of the dielectric coef¬cient can be signi¬cant, and there will be distinct
changes in physical properties associated with the in¬‚uence of phase transitions.
We can therefore write the Maxwell equation for the curl of the magnetic vector
H as
‚D ‚E
∇— H = j + = σE + 0 ,
˜ (8.3.7)
‚t ‚t
8.3 Continuum electrodynamics 137

where we have assumed that the rate of change of geological materials is slow
compared with the frequency of electromagnetic phenomena. For good electrical

conductors such as metals, the second term 0 E on the right-hand side of (8.3.7)
is negligible, and this will be the regime that prevails in the iron-rich core of the
Earth. Nevertheless, for less good conductors such as the silicate minerals in the
mantle this contribution can be signi¬cant at high frequencies ω > σ/ 0 .
For crystalline materials, the relations (8.3.4)“(8.3.6) are replaced by tensor

jk = σklEl, Dk = 0 klEl, Hk = μ0μklBl.
˜ (8.3.8)

The tensors σ, , μ will be invariant under the symmetry group of the crystal. From
thermodynamic arguments, the conductivity tensor σ, the dielectric tensor , and
the magnetic permeability tensor μ are all symmetric.
When a dielectric is deformed, the dielectric tensor has strain dependence
= ij + aijklekl, (8.3.9)

in terms of the linearised strain ekl. For certain classes of crystals deformation
is accompanied by the appearance of an electric ¬eld proportional to the
accompanying stress σ. In this case, the electric induction takes the form

Di = 0( ijEj + γi,klσkl), (8.3.10)

where the piezoelectric tensor has the symmetries

γi,kl = γi,lk, (8.3.11)

because of the symmetry of the stress tensor σ. An analogous relation to
(8.3.11) describes piezomagnetism. In each case the electromagnetic ¬elds have
to be determined along with the deformation. In general the ˜piezo™ coef¬cients
are rather small and it requires very considerable stress to produce signi¬cant
electromagnetic ¬elds. Nevertheless, such contributions have been invoked to
explain electromagnetic phenomena accompanying earthquakes.

8.3.3 Electromagnetic continuity conditions
At a material discontinuity where the conductivity and dielectric properties change,
there is the possibility of the build up of a surface charge density ρ and a surface
current ^ and discontinuities in the electromagnetic ¬eld quantities. By an analysis
similar to that in Section 3.2.2 using small volumes and surfaces crossing the
interface between media (1) and (2), with local unit normal n12, we can derive
the continuity conditions from the Maxwell equations (8.3.1).
The two equations for the divergence of the B, D ¬elds require the application
of the divergence theorem to a small volume crossing the interface. The normal
138 Continuum Equations and Boundary Conditions

component of the magnetic induction B is continuous across such a material
n12 · (B(2) ’ B(1)) = 0, (8.3.12)
but, as in Section 3.2.2, this means that the tangential components are not
constrained. In the presence of a layer of surface charge density ρ, the normal
component of the electric displacement D has a discontinuity equal to ρ,
n12 · (D(2) ’ D(1)) = ρ.
^ (8.3.13)
For E, H, whose curl enters the Maxwell equations, Stokes™ theorem is applied
to a surface crossing the interface. The tangential component of the electric vector
E is continuous across the material discontinuity
n12 — (E(2) ’ E(1)) = 0; (8.3.14)
whereas the tangential component of the magnetic vector H has a discontinuity
equal to the surface current density ^
n12 — (H(2) ’ H(1)) = ^
j. (8.3.15)

8.3.4 Energy equation for the electromagnetic ¬eld
From the Maxwell equations (8.3.1),
‚D ‚B
E · (∇ — H) ’ H · (∇ — E) = E · j + E · +H· , (8.3.16)
‚t ‚t
and from a standard vector identity, the left-hand side can be recognised as
∇ · (E — H),
’∇ · (E — H) = E · (∇ — H) ’ H · (∇ — E). (8.3.17)
Hence we can write (8.3.16) in the form
‚D ‚B
E· +H· + j · E + ∇ · (E — H) = 0. (8.3.18)
‚t ‚t
Consider then an arbitrary volume V with surface S and integrate (8.3.18)
throughout the volume
‚D ‚B
E· +H· dV + j · E dV + (E — H) · n dS = 0.
^ (8.3.19)
‚t ‚t
where n is the unit outward normal to S, and we have used the divergence to convert
the last term into a surface integral. The relation (8.3.19) is a direct consequence of
the Maxwell equations, and does not depend on the form of the material equations.
If we assume the material properties do not vary in time we can write (8.3.19) as
‚ 1
(E · D + B · H) dV + j · E dV + (E — H) · n dS = 0,
^ (8.3.20)
‚t 2
8.3 Continuum electrodynamics 139

which we can interpret as the energy conservation equation
+Q+ S · n dS = 0,
^ (8.3.21)
‚t S
W= V 2 (E · D + B · H) dV is the stored energy in the magnetic ¬eld,
Q= · E dV is the dissipation due to Joule heating, and
S = E — H is the Poynting vector.
S can be interpreted as the energy ¬‚ux per unit area, normal to E and H, associated
with the ¬eld.
Although ¬‚ow in the mantle might modify material properties it will do so on
a scale so much longer than the frequency of any electromagnetic phenomena that
the time derivatives of the conductivity σ and dielectric constant can be ignored.
For the Earth™s core ef¬cient stirring ensures a sustained radial pro¬le of variation
that does not vary signi¬cantly in time.

8.3.5 Electromagnetic disturbances
In a dielectric material described by (8.3.4)“(8.3.6), we can reduce the Maxwell
equations (8.3.1) to a suitable form to describe the propagation of electromagnetic
disturbances in a region without charge density. We can generate two equations for

∇ — H from the equation for the curl of the electric vector E,
1 ‚H
∇— ∇— E = ’∇ — , (8.3.22)
μ0μ ‚t
and the time derivative of the equation for the curl of the magnetic vector H,
‚j ‚2D
∇— = + 2. (8.3.23)
‚t ‚t ‚t
We can eliminate ∇ — ‚tH between (8.3.22) and (8.3.23) to give
1 ‚
’∇ — ∇— E = (σE) + 2 (
˜ E). (8.3.24)
μ0μ ‚t ‚t
Equation (8.3.24) can be simpli¬ed with the aid of the two vector identities
∇ — (∇ — v) = ∇(∇ · v) ’ ∇2v. (8.3.25)
∇ — (uv) = u∇ — v + ∇u — v,
Thus, we ¬nd
‚E 1
’∇(∇ · E) + ∇ E = μ0μσ + 0μ0 μ 2 + μ0μ∇ — E. (8.3.26)
‚t ‚t μ0μ
We use ∇ · D = 0 in association with the constitutive relation D = E, and then

∇( E) = ∇·E+ 0∇ · E = 0. (8.3.27)
0 0
140 Continuum Equations and Boundary Conditions

On combining (8.3.26) and (8.3.27) we have
∇ E ’ μ0μσ ’ 0μ0 μ =
’∇(log μ) — (∇ — E) ’ ∇(∇(log ) · E). (8.3.28)
For a homogeneous dielectric the material gradients vanish and (8.3.28) takes the
∇ E ’ μ0μσ ’ 0μ0 μ 2 = 0,
˜ (8.3.29)
‚t ‚t
with a comparable form for the magnetic vector H.
In the absence of conductivity, σ = 0, the equation (8.3.29) reduces to the wave
μ ‚2E
∇2E = , (8.3.30)
c2 ‚t2
= 1/c2.
where we have introduced the velocity of light in free space c through 0μ0
In the homogeneous dielectric we have
n2 ‚2E
∇ E= 2 , (8.3.31)
c ‚t2

with refractive index n = μ, so that the velocity of electromagnetic disturbances

in the material is c/ μ = c/n. The transverse electromagnetic waves have
energy ¬‚ow in the direction of the Poynting vector E — H.
For a disturbance with frequency ω represented via an complex exponential time
dependence exp(’iωt) we can write the equation (8.3.29) for the electric vector in
a conductive medium in the form
∇2E + ω ’ iμ0μσω E,
˜ (8.3.32)
which can be viewed as a wave equation for a complex refractive index,
corresponding to attenuation as the wave propagates. The effect can be readily
illustrated by considering a plane wave travelling in the x direction for which we
will have a y-component of the electric ¬eld Ey and a z-component of the magnetic
¬eld Hz,
n n
^ ^
Ey = exp iω x ’ iωt Ey, Hz = exp iω x ’ iωt Hz, (8.3.33)
c c
From (8.3.32) we have
ω2n2 μ ω2
’ 2 = ’ 2 ’ iωμμ0σ, ˜ (8.3.34)
c c
and thus the complex refractive index for an isotropic medium takes the form
n2 = μ + i . (8.3.35)

8.3 Continuum electrodynamics 141

We set n = n + in and then
n2 = n 2 ’ n 2
+ 2in n = μ + i . (8.3.36)

which can be solved for n , n . The propagating electric ¬eld takes the form
ω ω ^
Ey = exp ’ n x exp i n x ’ iωt Ey, (8.3.37)
c c
so that the plane wave decays as it travels in the x-direction. The amplitude falls to
e’1 after the wave has advanced a distance
d= , (8.3.38)
known as the penetration depth or skin depth, which is often a small fraction of the
wavelength in a vacuum. The relation of the electric and magnetic components is
cμμ0 ^ μμ0
^ ^ ^
Ey = ZHz = Hz = Hz, (8.3.39)
n 0 + iσ/ω
where Z is the intrinsic impedance of the medium. In the decaying transverse wave
the ratio of the amplitudes of Ey and Hz is everywhere the same |Z|, but a phase
difference φ is introduced between the electric and magnetic vectors, where
tan φ = n /n . (8.3.40)
For a very good conductor such as a metal, the conduction term dominates and
1/2 1/2
iσμ 1+i σμ
˜ ˜
n + in ≈ , (8.3.41)
0ω 0ω
with penetration depth
d= . (8.3.42)
The varying depth of penetration of electromagnetic disturbances with frequency
is exploited in the analysis of magnetotelluric signals to determine the conductivity
structure of the Earth. Such signals may be man-made or exploit the induced
electric and magnetic ¬elds produced by the varying magnetic ¬elds in the
ionosphere. Long period waves (∼104 s period) can penetrate into the transition
zone in the upper mantle, which appears as an increase in electrical conductivity.

8.3.6 Magnetic ¬‚uid dynamics
When a conducting ¬‚uid moves in a magnetic ¬eld, an electric ¬eld is induced and
electric currents ¬‚ow. The magnetic ¬eld exerts forces on the currents that can alter
the ¬‚ow, and the presence of the currents modi¬es the magnetic ¬eld. The result is
a complex interaction between electromagnetic effects and ¬‚uid dynamics.
142 Continuum Equations and Boundary Conditions

We will assume that the magnetic permeability of the conducting ¬‚uid is very
close to unity (μ ≈ 1), and work in terms of the E and B ¬elds. From Maxwell™s
equations (8.3.1)
∇ · B = 0, ∇— E = ’ ; (8.3.43)
the current density in a moving conductor is
j = σ(E + v — B).
˜ (8.3.44)
so that
∇ — B = σE + σv — B.
˜ ˜ (8.3.45)
The displacement current ‚D/‚t is negligible for phenomena whose time scale
is long compared with the time for passage of electromagnetic waves across the
region of interest. The neglect of the displacement current acts as a ¬lter to remove
electromagnetic waves from the situation in much the same way as we used the
incompressibility condition in Chapter 7 to suppress acoustic waves from ¬‚uid ¬‚ow.
We can use (8.3.45) to represent the electric ¬eld in terms of the magnetic ¬eld as
E= ∇ — B ’ v — B, (8.3.46)
and then from (8.3.43) and (8.3.44),
1 ‚B
∇— ∇— B =’ + ∇ — (v — B). (8.3.47)
σμ0 ‚t
In a homogeneous conductor with constant conductivity, we can simplify (8.3.47)
‚B 1
= ∇ — (v — B) + (8.3.48)
‚t σμ0
since ∇ · B = 0. The quantity » = (σμ0)’1 is known as the magnetic diffusivity of
the ¬‚uid.
The force acting on the magnetic ¬‚uid can be expressed as
F = μ0j — B = (∇ — B) — B = (B · ∇)B ’ 1 ∇B2. (8.3.49)

For a magnetic continuum the continuity equation for conservation of mass is
unchanged, but in the equation of motion we need to include the electromagnetic
force (8.3.49). The Navier“Stokes equation for the magnetic ¬‚uid thus has the form
ρ + (v · ∇)v = ∇σ + ρg + (∇ — B) — B
= ∇[σ ’ 1 B2] + μ0(B · ∇)B ’ ρ∇•, (8.3.50)

where the stress tensor σ for the viscous ¬‚uid is related to the velocity ¬eld through
(8.1.5), and ρ is the density of the ¬‚uid.
8.3 Continuum electrodynamics 143

As in Section 8.1, we can recast (8.3.51) into a direct expression of the
conservation of momentum
‚ ‚m
(ρvj) = ’ Π + ρgj (8.3.51)
‚xi ij
where the momentum ¬‚ux density tensor Πm now takes the form

Πm = ρvivj ’ σij ’ μ0(BiBj ’ 1 B2δij), (8.3.52)
ij 2

and the additional term relative to (8.1.8) is the Maxwell stress tensor for the
magnetic ¬eld
The rate of change of internal energy for the magnetic ¬‚uid has now to include
the in¬‚uence of Joule heating from current dissipation (j2/σ). As a result (8.1.9) is
replaced by
D 1
(∇ — B)2.
ρ U = h + ∇q + σ.∇v + (8.3.53)
Dt ˜0
The representation of the conservation of energy for the magnetic ¬‚uid has to
include the magnetic energy 1 B2/μ2 and also the electromagnetic energy ¬‚ux from
the Poynting vector E — B/μ0. Thus we ¬nd that
‚ 1B
= h ’ ∇ · Qm.
2 ρv + ρU + ρ• + (8.3.54)
‚t μ0
where the energy ¬‚ux density vector Qm is
Qm = ( 1 ρv2 + ρU)v ’ v · σ ’ k∇T
1 1
+ 2 B — (v — B) + (∇ — B) — B. (8.3.55)
μ0 ˜0
As in our discussion of non-magnetic ¬‚uids in Chapter 7 we can usually make the
approximation that the ¬‚uid is incompressible, and thereby both suppress acoustic
waves and simplify the equation system.

Frozen ¬‚ux
In the limiting case of very high conductivity, the magnetic diffusivity » ’ 0 and
(8.3.48) simpli¬es to
= ∇ — (v — B) = (B · ∇)v ’ (v · ∇)B ’ (∇ · v)B. (8.3.56)
After rearrangement of the terms in (8.3.56) we see that
D ‚B
B= + (v · ∇)B = (B · ∇)v ’ (∇ · v)B. (8.3.57)
Dt ‚t
With the aid of the continuity equation (8.1.1) we obtain
DB 1 DB B Dρ B
= ’2 = · ∇ v. (8.3.58)
Dt ρ ρ Dt ρ Dt ρ
144 Continuum Equations and Boundary Conditions

The rate of change of the quantity B/ρ satis¬es the same equation as that for a
material line element dx (2.4.8). As a result, if these vectors are initially in the
same direction they will remain parallel, with a constant ratio of length. Thus the
line of force B moves with the ¬‚uid particles that lie on it. We can regard the B-lines
as ˜frozen™ in the ¬‚uid, so that the topological structure of the ¬eld cannot change
with time.
Hence as any closed ¬‚uid contour moves over time it will not cut any B-lines
and so the ¬‚ux of B passing through the contour remains constant. This property
of frozen ¬‚ux for a highly conducting ¬‚uid has been invoked for the Earth™s core to
try to resolve ambiguities in the relation of the magnetic ¬eld and the velocity ¬eld
that are related to the requirement ∇ · B = 0 (see Chapter 15).

Consider the in¬‚uence of small wave disturbances in a constant magnetic ¬eld, when
dissipative effects from viscosity, thermal conductivity and electrical resistance can be
We assume high-frequency disturbances so that an adiabatic state is maintained and
introduce a small perturbation to the magnetic ¬eld, density and pressure
B = B0 + B , ρ = ρ0 + ρ , p = p0 + p , (8.3.e1)
where the subscript 0 indicates the constant equilibrium value. From the adiabatic con-
dition p = c2 ρ , where c0 is the sound speed. To a ¬rst-order approximation the
pertubations satisfy the magnetic equations
∇ · B = 0, = ∇ — (v — B0 ), (8.3.e2)
the continuity condition
+ ρ0 ∇ · v = 0, (8.3.e3)
and the Navier“Stokes equation
= ’c2 ∇ · p + B0 — (∇ — B ).
ρ0 (8.3.e4)
We now seek solutions for the additional disturbances in the form of plane waves with
slowness vector p and frequency ω so that, e.g.,
B ∝ exp[iω(p · x ’ t)],
the corresponding wavevector k = ωp. From (8.3.e2)“(8.3.e3) we ¬nd
’B = p — (v — B0 ) = v(p · B0 ) ’ B0 (p · v), (8.3.e5)
ρ = ρ0 p · v, (8.3.e6)
’ρ0 v + c2 ρ p = ’B0 — (p — B ). (8.3.e7)

The ¬rst equation (8.3.e5) requires that the vector B is perpendicular to the slowness
8.4 Diffusion and heat ¬‚ow 145

vector p, i.e., p · B = 0. We can eliminate ρ from (8.3.e7) with the aid of (8.3.e6) to
ρ0 c2 (p · v)p ’ v = ’B0 — (p — B ) = ’p(B0 · B ) ’ B (p · B0 ). (8.3.e8)

The properties of the waves are governed by the two coupled equations (8.3.e5), (8.3.e8)
for B and v.
Consider a disturbance transverse to the plane containing the slowness vector p and the
ambient magnetic ¬eld B0 , such that B · B0 = 0, p · v = 0 , then
’ρv = v(p · B0 )2 ;
B = ’v(p · B0 ), (8.3.e9)
√ √
so that p · B0 = ρ0 and the wavespeed u = 1/|p| = |B |/ ρ0 , where B is the
component of B parallel to p. The velocity

v = ’B / ρ0 , (8.3.e10)
and the dispersion relation for these transverse disturbances is

ω = ’B0 · k/ ρ0 , (8.3.e11)
with group velocity

= ’B0 / ρ0 , (8.3.e12)
which does not depend on the wavevector k.
A disturbance with v and B lying in the plane of p and B0 will have associated density
perturbations ρ = ρ(p · v). By eliminating B between (8.3.e5) and (8.3.e8) we ¬nd
ρ0 c2 p2 ’ 1 (p · v) = p2 (p · v)B2 ’ (p · B0 )(B0 · v) , (8.3.e13)
0 0

which leads to a quartic in the slowness p.
For very weak ambient magnetic ¬elds we have one class of waves for which c2 p2 ≈ 1
so that the magnetohydrodynamic waves reduce to ordinary sound waves. In the second

v ≈ ’B / ρ0 , (8.3.e14)
so that the situation is similar to (8.3.e10), but now with v, B in the plane of p and B0 .
In an incompressible ¬‚uid, for which ρ = 0, we have a single kind of wave propagation
with two independent directions of polarisation; both v and B are perpendicular to p

v = ’B / ρ0 , (8.3.e15)
such Alfv´ n waves are exact solutions of the equations and do not depend on the
assumption of small perturbations. •

8.4 Diffusion and heat ¬‚ow
The mechanical forms of deformation we have introduced include material ¬‚ow
and the transmission of wave phenomena, but in many circumstances diffusive
phenomena are important, as in the way that temperature effects spread through
a medium.
We have encountered the in¬‚uence of heat effects in our discussion of energy
146 Continuum Equations and Boundary Conditions

balance in Section 4.2, where we demonstrated that the local form of the
conservation of energy can be written as
ρ U = h ’ ∇x · q + σ · ∇x v, (8.4.1)
for internal energy density U, heat production h, heat ¬‚ow vector q, velocity ¬eld
v and stress tensor σ. The heat ¬‚ux q is related to the temperature gradient by
q = ’k∇x T, (8.4.2)
where k is the thermal conductivity.
The internal energy density U can be related to temperature through a thermal
capacity C and from (7.1.15), in the absence of deviatoric strain, the equation for
the temperature ¬eld can be written as
D ‚
ρ (CT ) = ρ (CT ) + v · ∇x (CT ) = h + ∇x · [k∇x T ] . (8.4.3)
Dt ‚t
The particular thermal capacity C depends on the thermodynamic state; normally
we would use Cp the thermal capacity at constant pressure. In equilibrium,
‚T/‚t = 0 and so
ρv · ∇x (CpT ) ’ ∇x · [k∇x T ] = h. (8.4.4)
In a material at rest with constant properties Cp, k the equation for temperature
simpli¬es substantially to

T = h + κH ∇2T . (8.4.5)
A comprehensive set of analytical solutions to heat ¬‚ow problems based on (8.4.5)
are to be found in Carslaw & Jaeger (1959).

8.4.1 Equilibrium heat ¬‚ow
Consider the simple case where the temperature T and the thermal conductivity k
depend only on depth z. Then the conduction of heat in steady state is governed by
d d
k(z) T (z) = 0, (8.4.6)
dz dz
which can be readily integrated to yield

T (z) = a + b . (8.4.7)
The constants a, b have to be determined by the boundary conditions in the speci¬c
situation. Suppose we have a layer whose boundaries are at ¬xed temperatures
T (zu) = Tu, T (zl) = Tl, then
T (z) = Tu + (Tl ’ Tu) . (8.4.8)
8.4 Diffusion and heat ¬‚ow 147

In the special case where k is constant, (8.4.8) reduces to a linear gradient in
temperature, independent of k:
z ’ zu
T (z) = Tu + (Tl ’ Tu) . (8.4.9)
zl ’ z u
When heat sources are evenly distributed through a region so that
d d
k(z) T (z) = ’ho, (8.4.10)
dz dz
the solution (8.4.7) is augmented by the in¬‚uence of the heat generation
z z

T (z) = a + b ’ ho ;
dζ (8.4.11)
k(ζ) k(ζ)
once again a, b are to be found from the imposed boundary conditions. For a
homogeneous layer the additional temperature contribution is quadratic in depth.
The exponentially decaying heat production h = hre’z/L represents an
approximation to the likely radiogenic heat production in the Earth. For this
situation the temperature distribution in z > 0 satis¬es
d d
= hre’z/L.
k(z) T (z) (8.4.12)
dz dz
The heat ¬‚ow
T (z) = hrLe’z/L + qm,
q(z) = ’k(z) (8.4.13)
if the heat ¬‚ow at great depth is qm. Thus the surface heat ¬‚ow
q(0) = qm + hrL, (8.4.14)
and is linear in the surface heat production rate.

For a homogeneous material the equilibrium temperature satis¬es
k∇2 T = ’h.
Cylindrical and spherical problems can be tacked with the appropriate form for the
Laplacian operator ∇2 .
Cylindrical symmetry: for steady radial ¬‚ow in the presence of heat production hc , the
temperature distribution is controlled by
1‚ ‚T
k = ’hc ,
‚ ‚
where is the cylindrical radial variable. The general solution of this differential
equation is
T ( ) = A + B ln ’ hc .
148 Continuum Equations and Boundary Conditions

For a solid cylinder the term in ln is not allowed, but needs to be included for problems
with hollow cylinders.
Spherical symmetry: the radial temperature distribution is controlled by
1‚ ‚T
k = ’hs ,
r2 ‚r ‚r
with solution
B hs 2
T (r) = A + + r;
r 6k
once again the singular term has to be excluded if the domain includes the origin r = 0.•

A simple model of the motion of a ¬‚uid through a porous rock is provided by the
scenario where there are no heat sources, but the advective term is retained. For a
homogeneous region with porosity g and a ¬‚uid with vertical transport speed vz,
the temperature has to satisfy
γ = 2, (8.4.15)
dz dz
where γ = ρCpvzg/k = vzg/κH . The temperature solution is
T (z) = c ’ d (8.4.16)
with c, d determined by the boundary conditions. The departure of the temperature
distribution (8.4.16) from the linear solution for conductive heat transport in a
uniform medium can be used to infer the value of vz.

8.4.2 Time-varying problems
We continue to consider problems that depend just on the depth coordinate z, but
now allow for the evolution of the temperature ¬eld with time.
Consider the response to an impulsive source of heat h(z, t) = δ(z)δ(t) in a
uniform medium, the Green™s function for the temperature has to satisfy
’ κ 2 = δ(z)δ(t), (8.4.17)
‚t ‚t
where for notational convenience we have written κ for the thermal diffusivity κH .
We can ¬nd a solution by means of a Fourier transform. We express T (z, t) as a
combination of frequency and wavenumber components as
∞ ∞
1 ¯
dω T (q, ω)ei(qz’ωt),
TG(z, t) = dq (8.4.18)
2π ’∞ ’∞

and then
¯ ¯
’iωT (q, ω) + κq2T (q, ω) = 1/2π. (8.4.19)
8.4 Diffusion and heat ¬‚ow 149

We now have to evaluate the multiple integrals over vertical wavenumber q and
frequency ω
∞ ∞
TG(z, t) = .
dq dω 2 (8.4.20)
(2π)2 κq ’ iω
’∞ ’∞
Performing the frequency integral ¬rst using contour integration we can extract the
contribution from the polar residue at ω = ’iκq2, when t > 0, to give

1 2 t)
dq ei(qz’κq
TG(z, t) = H(t), (8.4.21)
2π ’∞

where H(t) is the Heaviside step function. We can rewrite (8.4.21) as

1 ’z2 /4κt 2
dq e’κt(q+iz/2κt) ,
TG(z, t) = H(t)
2π ’∞
1 2 /4κt
= H(t), (8.4.22)
with the aid of a change of variables. Subsequent to the impulse of heat the
temperature spreads out with a Gaussian pro¬le whose scale length depends on
the inverse square root of the time t’1/2.
We can ¬nd a comparable solution for the situation where the initial temperature
pro¬le T0(z) is speci¬ed at t = 0,

T (z, t) = dζ T0(ζ)TG(z ’ ζ, t). (8.4.23)
An important special case is provided by
T0(z) = T0 H(z), (8.4.24)
which is equivalent to the sudden juxtaposition of two uniform half-spaces with
different temperatures. The resulting temperature distribution

∞ z/ 4κt
T0 2
2 /4κt 2
dζ e’(z’ζ) dξ e’ξ ,
= 1 T0 √
T (z, t) = √ 2 π
4πκt 0 0
= 1 T0 1 + erf √ , (8.4.25)
where the error function erf(x) is de¬ned as
2 2
dξ e’ξ ,
erf(x) = √ (8.4.26)
π 0
and erf(’x) = ’erf(x).
The situation of the cooling of a half-space initially at temperature T0 with a
requirement that the surface temperature vanishes is equivalent in z > 0 to the
T0(z) = T0 sgn(z), (8.4.27)
150 Continuum Equations and Boundary Conditions

and then
z ’z
T (z, t) = 1 T0 1 + erf ’ 1 T0 1 + erf
√ √
2 2
4κt 4κt

= T0 erf . (8.4.28)

The relative penetration of temperature ¬‚uctuations imposed at the surface can be assessed
by the analysis of a harmonic ¬eld.
Consider a temperature ¬eld with a frequency dependence exp(’iωt), then
d2 T
κ + iωT = 0
with solutions of the form
1 1
ω ω
2 2
T (z, t) = exp ’ z exp i z ’ iωt .
2κ 2κ
The temperature behaviour is in the form of a diffusive wave, decaying as it penetrates
into the medium (cf. (8.3.39)). The rate of decay with depth is proportional to the square
root of frequency (ω1/2 ). This means that for a comparable reduction of the amplitude
of temperature variation relative to the surface, the penetration of the yearly cycle will

be 365 greater than the diurnal variation. The diurnal cycle can be perceptible through
a masonry wall about 20 cm thick, whereas in a 4 m deep cellar the yearly variations
are nearly out of phase with the surface variations and serve to minimise the overall
temperature ¬‚uctuations. •

For a situation where boundary conditions are imposed on a layer the solution
may be found in terms of a Fourier series in z rather than a Fourier transform. Thus,
for example, for the situation where a layer of thickness L has an initial temperature
Ta and boundary conditions T = 0 at z = 0, T = Ta at z = L, a solution for the
subsequent temperature distribution can be found as

z 2 nπz
T (z, t) = Ta + exp ’ ,
sin (8.4.29)
L nπ L

with heat ¬‚ow Q0(t) at z = 0 of

Q0(t) = 1+2 exp ’ . (8.4.30)
Part II
From the Atomic Scale to the Continuum

The continuum equations we have encountered in previous chapters make use
of the device of a constitutive equation representing the relation between stress
and strain at the macroscopic level. Such constitutive equations are therefore
the intermediaries by which properties at the atomic level are translated to their
implications at much larger length scales. Behaviour at the atomic scale depends
on the properties of the constituent crystals, e.g., elastic moduli, and on departures
of the material from ideality. Grain boundaries, vacancies, dislocations and other
forms of defects play their role in controlling transport properties in the crystals.
Such transport phenomena in turn control the way in which composite materials
respond to externally imposed forces.
A wide variety of semi-empirical constitutive relation have been devised to
represent speci¬c aspects of material behaviour. The variety of functional forms
represents the balance between different classes of deformation mechanisms at the
atomic level. The net result is a non-linear relation between stress and strain (or

9.1 Transport properties and material defects
9.1.1 Grains and crystal defects
The polycrystalline aggregates that make up the Earth™s silicate mantle are
composed of multiple minerals, each with ¬nite grain size within which there may
be further classes of interruption of the regular crystal lattices.
The distribution of grain size with depth is therefore a very important, but
imperfectly known, property of the mantle. It is tempting to assume that the
millimetre to centimetre size crystals of ma¬c materials brought to the surface in
xenoliths are representative of the depths from which they have been detached.
Kimberlite and lamprophyre pipes show every evidence of very rapid eruption so
that material has been projected to the surface very rapidly, but interaction with the
eruptive process cannot be excluded.
Some exposures of former upper mantle rocks in mountain belts suggest that
the original crystal sizes were rather large. Nevertheless, the common view is

154 From the Atomic Scale to the Continuum

that grain sizes in the continental lithosphere are unlikely to be larger than a few
centimetres and would be expected to decrease in size in the asthenospheric regime.
Faul & Jackson (2005) have included a grain size dependence in a model of the
behaviour of shear moduli and attenuation with temperature for olivine systems;
¬ts to velocity pro¬les for both continental and oceanic domains suggest the need
for a signi¬cant increase in grain size in the lower part of the asthenosphere. In the
transition zone and below there are no direct indicators of grain size. However,
there are suggestions that there should be a grain size reduction at the phase
transitions in silicate minerals that lead to closer packed structures, particularly
just below the 660 km discontinuity due to the spinel to perovskite phase change.
Increasing evidence suggests the presence of signi¬cant amounts of water at depth,
and variations in water content could well affect grain size.
Within the crystal grains, the lattice structures are not perfect; there may be point,
line, surface or volume defects that disturb locally the regular arrangements of the
atoms. Minor departures from the theoretical composition (non-stoichiometry) will
lead to point defects, either vacancies due to the absence of one atomic species,
or interstitials where there too many atoms of one variety. A common feature
of the silicate minerals is the substitution of small amounts of different atoms.
The presence of some iron atoms in a magnesio-silicate will lead to some lattice
strain around the iron sites, that may need to be accommodated by other defects.
Even more signi¬cant is the case when elements with different valency such as
aluminium are introduced into the lattice since the charge distribution will also be
affected (see, e.g., Li et al., 2005, for a discussion of substitution in perovskite and
the effect on material properties).
The most signi¬cant class of line defects is the presence of dislocations that
offset parts of the regular crystal structure. Such dislocations can be described by
two vectors: the line direction l, which is the tangent to the defect at any point, and
the Burgers vector b, which describes the closure failure in a circuit made around
the line defect. The line defect moves in the glide plane with normal n = l — b
(Figure 9.1a). In an edge dislocation (Figure 9.1b) an extra half plane of atoms is
inserted into the lattice, so that the Burgers vector b is perpendicular to the line
direction l. For a screw dislocation (Figure 9.1c), which may be thought of as
a shear of part of the lattice, b and l are parallel so that glide can occur on any
plane including l. The most general dislocation will have both both edge and screw
components, as in Figure 9.1a.
Dislocation lines can end at the surface of a crystal and at grain boundaries, but
never inside a crystal. Dislocations must therefore form closed loops or branch into
other dislocations (for which the net Burgers vector must vanish). The result is a
complex net of dislocations permeating the crystal that will be affected by crystal
The importance of dislocations is that they provide an energetically favourable
way for a crystal lattice to take up deformation. Rather than the high energy
9.1 Transport properties and material defects 155



b b
(b) (c)

Figure 9.1. (a) The geometry of a dislocation line defect is described by the Burgers vector
b, line direction l and the glide plane with normal n. The Burgers vector describes the
displacement produced by the passage of the dislocation. (b) Con¬guration of an edge dis-
location, where the Burgers vector is perpendicular to the line direction. (c) Con¬guration
of a screw dislocation, where the Burgers vector is parallel to the line direction.

requirement to force one atom past another, a dislocation can progress by the
breaking and re-formation of atomic bonds. Thus, e.g., an edge dislocation can
move in the slip direction de¬ned by the Burgers vector b by the transfer of the
additional half plane of atoms under relatively small applied stress.
Estimates of the stress and strain ¬eld around a dislocation can be obtained by
using continuum approximations. A simple model of a screw dislocation in an
isotropic medium is provided by the distortion of a cylindrical ring surrounding the
line vector l. For a uniform shear strain along l, taken as the z-axis, the associated
stress ¬eld for displacement b at radius r from the dislocation line is
σzθ = σθz = , (9.1.1)
where G is the shear modulus. All other stress tensor components vanish. The
stress ¬eld thus has shear components σθz in radial planes parallel to the line
direction z, and σzθ in planes normal to the z-axis perpendicular to the radius
vector from the dislocation line. The immediate neighbourhood of the dislocation
line has to be excluded; within the core radius r0 the atomic displacements and
156 From the Atomic Scale to the Continuum

forces have to be treated explicitly. The elastic strain energy associated with the
screw dislocation is thus approximately
Gb2 r1
Eel = ,
ln (9.1.2)
4π r0
where r1 is the external diameter of the crystal. Real crystals are not isotropic and
in consequence the stress ¬eld is more complex than in the simple model, but the
energy estimate (9.1.2) is still useful.
A comparable model for an edge dislocation has a radial displacement b in a
cylindrical ring. Now the stress ¬eld has both dilatational and shear components.
The largest normal stress acts parallel to the slip vector b. In an isotropic medium
the estimate for the elastic energy for the edge dislocation model is larger than for
a screw dislocation by 1/(1 ’ …), where … is Poisson™s ratio; this result re¬‚ects the
additional stresses in the edge dislocation case. From Section 5.2.5 typical values of
… are around 0.3. Thus the energetics of deformation will not be heavily dependent
on the dislocation con¬guration.

9.1.2 General transport properties
We can describe the transfer of physical properties through a solid matrix through
the action of ˜carriers™ that are associated with speci¬c properties. Electrical
conductivity depends on the transfer of charge by electrons and holes and
the conduction of heat depends on lattice vibrations (quantised as phonons).
We therefore have to take account of carrier motion (diffusion), the in¬‚uence
of external ¬elds and scattering of the carriers from the atomic lattice or its
We introduce the local concentration of carriers fk (x) for a state k, with velocity
vk . The total rate of change of the carrier concentration, which determines the local
transport properties, is given by
‚fk ‚fk ‚fk ‚fk
= + + , (9.1.3)
‚t ‚t ‚t ‚t
diff ¬eld scatt

and will vanish in a steady-state situation.
The ¬rst component on the right hand side of (9.1.3), ‚fk /‚t|diff , describes the
migration of the carriers by diffusion. In the presence of a temperature gradient the
diffusive term
‚fk ‚fk
= ’vk · ∇fk = ’vk · ∇T . (9.1.4)
‚t ‚T

The second contribution ‚fk /‚t|¬eld arises mainly from the interaction of the
electromagnetic ¬eld with charged carriers, such as electrons, holes or ions, and is
of particular signi¬cance for electrical conductivity.
The ¬nal contribution from scattering, ‚fk /‚t|scatt , can be described through a
9.1 Transport properties and material defects 157

relaxation time that is inversely proportional to a weighted integral of scattering
probability over all directions. The deviations from the reference state f0 are
determined by
‚fk 1 1
fk ’ f0 ,
= gk = (9.1.5)
‚t „ „
with an exponentially decaying solution
gk (t) = gk (0)e’t/„. (9.1.6)
The mean free path for the carriers Λ associated with the scattering is given by
Λ = „vk . Where several different classes of scattering phenomena are present,
such as the various classes of defects, the effects can be described approximately
as due to the addition of contributions with different relaxation times.
In a perfect lattice the dominant scattering contribution comes from thermal
vibrations of the lattice. This effect is very important in metals at low temperatures
with a low concentration of defects. However, for the semi-conducting materials
of the silicate mantle at relatively high temperatures the in¬‚uence of defects is

9.1.3 Atomic diffusion
The driving force for the diffusion of an atomic species A comes from gradients in
the local concentration cA, i.e., the number of atoms per unit volume. The local
¬‚ux of the species A,
JA = ’DA∇cA, (9.1.7)
where DA is the diffusivity of species A, with units m2 s’1 . Mass conservation
= ’∇ · JA = ∇ · (DA∇cA) . (9.1.8)
For a dilute concentration of species A, the diffusivity is only weakly dependent on
cA and then we obtain a diffusion equation (Fick™s equation),
= DA∇2cA. (9.1.9)
The concentration cA can be related to the chemical potential μA through the
atomic fraction of the species NA = cA/n, where n is the total number of atoms
per unit volume. The chemical potential
μA = μ0 + RT ln[γANA] = μ0 + RT ln , (9.1.10)
where μ0 is the chemical potential in the standard state and γA is the activity
coef¬cient for species A.
158 From the Atomic Scale to the Continuum

Direct diffusion through the lattice carries a high energy penalty. It is much more
ef¬cient for diffusion to occur by the migration of vacancies. Even so an energy
barrier has to be overcome to shift an atom from one site to another that requires
thermal activation to accomplish. Similar considerations apply for other classes of
interaction, including the migration of charged species where electrical neutrality
has to be maintained. In general, the diffusivity has a temperature dependence that
takes the Arrhenius form
D = D0 exp[’H— /RT ], (9.1.11)
in terms of an activation enthalpy H— , for absolute temperature T and gas constant
R. The activation enthalpy H— = E— + pV — is a function of pressure, since the
effect of pressure is to reduce the number of vacancies and to increase the potential
barrier between lattice sites.
An alternative expression for the diffusivity is in terms of the homologous
temperature (T/TM ), the ratio of the temperature T to the melting temperature TM ,
D = D0 exp[’aTM /T ], (9.1.12)
where the pressure dependence of diffusivity is included through the variation of
the melting temperature with pressure.

9.2 Lattice vibrations
Atoms in a crystal vibrate about their equilibrium positions, and analysis of such
vibrations allows a link to be made between thermal and elastic properties. The
low-frequency parts of the vibrational spectrum are equivalent to elastic waves; the
higher-frequency parts correspond to thermal vibrations.
We de¬ne the lattice displacement usl as the displacement of the sth atom from
its equilibrium position in the unit cell with origin vector l. The kinetic energy
associated with such vibrations is
K= 2 Ms|‚tusl | , (9.2.1)

where Ms is the mass of the sth atom. The potential energy of the distorted
crystal can be expanded about the equilibrium position, where there are no atomic
V = V0 + us l . (9.2.2)
‚usl ‚us l
ss ,ll

The equation of motion of the atoms can therefore be written as
Ms‚ttusl = ’ us l = ’ Gss ;ll · usl , (9.2.3)
‚usl ‚us l
sl sl
9.2 Lattice vibrations 159

ω 3n-3 optical modes

3 acoustic modes


. 6
( 16)