w = 0. (7.6.19)

The three components of the momentum equation are

‚2u

1 ‚p

’2 ©v = ’ +ν 2,

ρ ‚x ‚z

‚2v

1 ‚p

2 ©u = ’ + ν 2,

ρ ‚y ‚z

‚p

0= , (7.6.20)

‚z

and hence

‚p ‚p0 ‚p ‚p0

= = 0, = = 0. (7.6.21)

‚x ‚x ‚y ‚y

The equations to be solved for the ¬‚ow are

‚2u

’2©v = ν 2 ,

‚z

‚2v

’2©(u0 ’ u) = ν 2 , (7.6.22)

‚z

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)

where

” = (ν/©)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.

8

Continuum Equations and Boundary

Conditions

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.

131

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)

‚t

this relation is equivalent to the Lagrangian form (2.5.7),

Dρ

+ ρ∇ · v = 0; (8.1.2)

Dt

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

‚vj

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

‚uk

σ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

D

ρ U = h ’ ∇x q + σ.∇x v, (8.1.9)

Dt

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

D

(U + 1 v2) = h ’ ∇x q + ∇x · [σ · v] + g · v,

ρ (8.1.11)

2

Dt

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

D

g · v = ’v · ∇x • = ’ •, (8.1.13)

Dt

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

satis¬ed.

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

Dt

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)

‚t

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)

2

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

energy.

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

sides.

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

L

dpb

= , (8.2.7)

Tb ”v

dTb

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)

‚t

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

dependence.

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)

0

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

relations

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

0

= ij + aijklekl, (8.3.9)

ij

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

j,

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

discontinuity

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 ^

j,

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

V V S

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

V V S

8.3 Continuum electrodynamics 139

which we can interpret as the energy conservation equation

‚W

+Q+ S · n dS = 0,

^ (8.3.21)

‚t S

where

1

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

Vj

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

‚H

∇— = + 2. (8.3.23)

‚t ‚t ‚t

We can eliminate ∇ — ‚tH between (8.3.22) and (8.3.23) to give

‚2

1 ‚

’∇ — ∇— E = (σE) + 2 (

˜ E). (8.3.24)

0

μ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

‚2

‚E 1

2

’∇(∇ · 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

0

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

‚2E

‚E

2

∇ E ’ μ0μσ ’ 0μ0 μ =

˜

‚t2

‚t

’∇(log μ) — (∇ — E) ’ ∇(∇(log ) · E). (8.3.28)

For a homogeneous dielectric the material gradients vanish and (8.3.28) takes the

form

‚2E

‚E

2

∇ 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

˜

equation

μ ‚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

2

∇ 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

μ2

∇2E + ω ’ iμ0μσω E,

˜ (8.3.32)

c2

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)

0ω

8.3 Continuum electrodynamics 141

We set n = n + in and then

σμ

n2 = n 2 ’ n 2

+ 2in n = μ + i . (8.3.36)

0ω

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

c

d= , (8.3.38)

ωn

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

1/2

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

then

1/2 1/2

iσμ 1+i σμ

˜ ˜

=√

n + in ≈ , (8.3.41)

0ω 0ω

2

with penetration depth

1/2

2

d= . (8.3.42)

σωμμ0

˜

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

∇ · B = 0, ∇— E = ’ ; (8.3.43)

‚t

the current density in a moving conductor is

j = σ(E + v — B).

˜ (8.3.44)

so that

1

∇ — B = σE + σv — B.

˜ ˜ (8.3.45)

μ0

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

1

E= ∇ — B ’ v — B, (8.3.46)

σμ0

˜

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)

to

‚B 1

∇2B,

= ∇ — (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)

2

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 · ∇)v = ∇σ + ρg + (∇ — B) — B

‚t

= ∇[σ ’ 1 B2] + μ0(B · ∇)B ’ ρ∇•, (8.3.50)

2

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

‚t

where the momentum ¬‚ux density tensor Πm now takes the form

ij

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

σμ2

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

0

2

the Poynting vector E — B/μ0. Thus we ¬nd that

2

‚ 1B

12

= h ’ ∇ · Qm.

2 ρv + ρU + ρ• + (8.3.54)

22

‚t μ0

where the energy ¬‚ux density vector Qm is

Qm = ( 1 ρv2 + ρU)v ’ v · σ ’ k∇T

2

1 1

+ 2 B — (v — B) + (∇ — B) — B. (8.3.55)

σμ2

μ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

‚B

= ∇ — (v — B) = (B · ∇)v ’ (v · ∇)B ’ (∇ · v)B. (8.3.56)

‚t

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

EXAMPLE: MAGNETOHYDRODYNAMIC WAVES

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

neglected.

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

0

pertubations satisfy the magnetic equations

‚B

∇ · B = 0, = ∇ — (v — B0 ), (8.3.e2)

‚t

the continuity condition

‚ρ

+ ρ0 ∇ · v = 0, (8.3.e3)

‚t

and the Navier“Stokes equation

‚v

= ’c2 ∇ · p + B0 — (∇ — B ).

ρ0 (8.3.e4)

0

‚t

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)

0

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

yield

ρ0 c2 (p · v)p ’ v = ’B0 — (p — B ) = ’p(B0 · B ) ’ B (p · B0 ). (8.3.e8)

0

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)

‚k

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

0

so that the magnetohydrodynamic waves reduce to ordinary sound waves. In the second

class,

√

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

with

√

v = ’B / ρ0 , (8.3.e15)

such Alfv´ n waves are exact solutions of the equations and do not depend on the

e

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

D

ρ U = h ’ ∇x · q + σ · ∇x v, (8.4.1)

Dt

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)

x

‚t

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

z

dζ

T (z) = a + b . (8.4.7)

k(ζ)

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

z

dζ/k(ζ)

zu

T (z) = Tu + (Tl ’ Tu) . (8.4.8)

zl

dζ/k(ζ)

zu

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

ζ

dζ

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

d

T (z) = hrLe’z/L + qm,

q(z) = ’k(z) (8.4.13)

dz

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.

EXAMPLE: STEADY-STATE PROBLEMS FOR A CYLINDER OR SPHERE

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

2

T ( ) = A + B ln ’ hc .

2k

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

r2

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

d2T

dT

γ = 2, (8.4.15)

dz dz

where γ = ρCpvzg/k = vzg/κH . The temperature solution is

1

eγz,

T (z) = c ’ d (8.4.16)

γk

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

‚2TG

‚TG

’ κ 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 ω

∞ ∞

ei(qz’ωt)

1

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)

e

2π ’∞

1/2

1 2 /4κt

e’z

= H(t), (8.4.22)

4πκt

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

z

= 1 T0 1 + erf √ , (8.4.25)

2

4κt

where the error function erf(x) is de¬ned as

x

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

choice

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

z

√

= T0 erf . (8.4.28)

4κt

EXAMPLE: DIURNAL AND ANNUAL TEMPERATURE FLUCTUATIONS

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

dz2

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

∞

n2π2κt

z 2 nπz

T (z, t) = Ta + exp ’ ,

sin (8.4.29)

L2

L nπ L

n=1

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

∞

n2π2κt

kTa

Q0(t) = 1+2 exp ’ . (8.4.30)

L2

L

n=1

Part II

EARTH DEFORMATION

9

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

strain-rate).

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

153

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

deformation.

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

(a)

n

l

b

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

Gb

σzθ = σθz = , (9.1.1)

2πr

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

imperfections.

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

diff

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

k

determined by

‚fk 1 1

fk ’ f0 ,

= gk = (9.1.5)

k

‚t „ „

scatt

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

dominant.

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

requires

‚cA

= ’∇ · JA = ∇ · (DA∇cA) . (9.1.8)

‚t

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

‚cA

= DA∇2cA. (9.1.9)

‚t

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

γAcA

μA = μ0 + RT ln[γANA] = μ0 + RT ln , (9.1.10)

A A

n

where μ0 is the chemical potential in the standard state and γA is the activity

A

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

2

1

K= 2 Ms|‚tusl | , (9.2.1)

sl

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

displacements,

‚2V

uT

V = V0 + us l . (9.2.2)

sl

‚usl ‚us l

ss ,ll

The equation of motion of the atoms can therefore be written as

‚2V

Ms‚ttusl = ’ us l = ’ Gss ;ll · usl , (9.2.3)

‚usl ‚us l

sl sl

9.2 Lattice vibrations 159

ω 3n-3 optical modes

}gap

3 acoustic modes

L

T1

T2