q

π/a

Figure 9.2. Dispersion curve of an in¬nite lattice with n atoms per unit cell, mapped into

half the ¬rst Brillouin zone. There are 3 acoustic modes (a longitudinal P mode and two

transverse S modes polarized at right angles), and 3n ’ 3 optical modes with a gap of

excluded frequencies at the edge of the Brillouin zone.

where we have introduced the tensor G summarising the atomic forces. Such forces

between atoms depend only on the relative position of l and l so that

Gss ;ll = Gss (h) where h = l ’ l. (9.2.4)

With this force representation the equations of motion can be cast as

Ms‚ttuT = ’ Gss (h) · us,l+h , (9.2.5)

sl

sh

and we can seek a solution in terms of functions that preserve the translational

invariance between different unit cells,

usl = eiq·l’iωtUs ,q . (9.2.6)

The possible values of the wave vector q are determined by the boundary conditions

on the crystal; the periodic nature of the crystal means that the wave vectors q can

be restricted to a Brillouin zone in wavenumber space (see, e.g., Poirier, 1991). The

frequency ω for the wavevector q must satisfy

¯

ω2MsUs,q = Gss (h)eiq·h Us ,q = Gss (q)Us ,q , (9.2.7)

s s

h

¯

where Gss is the Fourier transform of the force tensor G.

For N unit cells with n atoms there are Nn vectorial equations (9.2.3), but the

periodic form (9.2.6) for usl reduces the system to just n vector equations. The 3n

frequencies ω and the eigenvectors Us,q are to be found from the equation system

¯

Gss (q) ’ ωMsδss Us ,q = 0. (9.2.8)

s

160 From the Atomic Scale to the Continuum

The lowest frequency modes, with longest wavelength, correspond to movements

of the entire lattice and can be regarded as equivalent to elastic waves. The

three acoustic branches correspond to two transverse S waves with orthogonal

polarisation, and longitudinal motion in P waves. The remaining 3n ’ 3 optical

branches involve more localised disturbances with motion of neighbouring planes

out of phase for small |q|. The optical mode frequencies generally fall in the range

of infrared or visible light, and can lead to optical absorption. Near the edge of the

Brillouin zone, there is a gap in frequency between the optical and acoustical modes

(Figure 9.2). At high temperatures, the optical modes carry a signi¬cant part of

the thermal energy and their propagation speed depends on frequency; interactions

between the modes are also important.

For small wavenumbers |q|, the frequency of the acoustic modes is proportional

to the wavenumber and the modes are non-dispersive: the group velocity ‚ω/‚q is

equal to the phase velocity ω/q, which is the macroscopic phase velocity. Because

crystals are generally anisotropic, the phase velocity will be a function of the

direction of the wave vector q. As the boundary of the Brillouin zone is approached

there will be dispersion for the acoustic modes for any given direction as indicated

in Figure 9.2.

Lattice speci¬c heat

The thermal excitation of lattice waves makes a distinct contribution to the speci¬c

heat of the solid. We have shown in (9.2.7) that we are able to reduce the lattice

vibration equations to a set of 3n independent harmonic oscillators for each wave

vector q. Using a Bose“Einstein quantum representation for each oscillator, the

average total energy of the system for temperature T is given by

hωq ehωq /kB T

NVL 3

E= ,

d q hω /k T (9.2.9)

(2π)3 e q B ’1

P

where kB is the Boltzmann constant and VL is the volume of a unit cell. The

summation over P includes all modes. We introduce the density of states g(ω) as

the number of modes per unit frequency range, and can then write

hω ehω/kB T

E = NVL dω g(ω) hω/k T , (9.2.10)

B ’1

e

The speci¬c heat is found by differentiating the energy with respect to temperature,

(hω/kBT )2 ehω/kB T

1 ‚E

CV = = dω g(ω)kB , (9.2.11)

2

NVL ‚T ehω/kB T ’1

A useful, but restricted, approximation can be made with a highly simpli¬ed

representation of the density of states g(ω) under the assumptions:

(i) only the acoustic modes need to be considered, and all are assigned the same

constant sound velocity vm;

9.2 Lattice vibrations 161

(ii) the Brillouin zone, which limits the range of allowed values of q, is replaced by

a sphere of the same volume in wavenumber space.

There is therefore a maximum wavenumber for the lattice waves, the Debye

wavenumber qD such that

1/3

(2π)3 6π2

3

4

3 πqD = , qD = .

i.e., (9.2.12)

VL VL

The associated cut-off frequency

’1/3

1 2

1/3

ωD = qDvm, with vm = 3 +3 , (9.2.13)

v3 vS

P

where vP and vS are the P and S wavespeeds. The Debye lattice spectrum is

therefore assumed to be quadratic in frequency

3nNVL 2

gD(ω)dω = ω dω, 0 ¤ ω ¤ ωD. (9.2.14)

2π2v3

m

With the substitution z = hω/kBT and the approximation (9.2.14) for the spectrum

g(ω) we can express the speci¬c heat as

3 ˜D /T

z4 ez

T

Cv = 3nNkB 3 ,

dz (9.2.15)

2

˜D ez ’ 1

0

where the Debye temperature ˜D = hωD/kB. When the temperature is large, i.e.,

T ˜D, the upper limit on the integral is small and the integrand can be expanded

in powers of z to give

3

˜D /T

z4 ˜D

1

dz 2 = , (9.2.16)

3

z T

0

so that we obtain the classical value for the speci¬c heat

Cv = 3nNkB, T ˜D, (9.2.17)

the Dulong and Petit Law. For very low temperatures, the upper limit of the

integral can be regarded as tending to in¬nity, and so the integral tends to a constant

(4π4/15). The speci¬c heat is then

3

12π4 T

Cv ∼ nNkB . (9.2.18)

15 ˜D

Despite the drastic simpli¬cations the Debye formula (9.2.15) works well for most

solids. The relation g(ω) ∝ ω2 should hold near ω = 0, when the material

behaves like an elastic continuum, but the sharp cut-off at ωD is not justi¬ed.

Exact calculations show that the lattice spectrum has a large spread with several

peaks corresponding to the very different propagation velocities of modes with

different polarisations. A slightly better physical approximation is to use a sum

162 From the Atomic Scale to the Continuum

of Debye-type contributions for each of the acoustic branches, but the simplicity of

the single Debye temperature is then lost.

The vibrational motion of the atoms leads to a thermal contribution to the

pressure, that can be estimated from the Debye treatment as

3 ˜D /T

z3

γ γ T

Pth (V, T ) = E D= 9nNkBT ,

dz z (9.2.19)

V V ˜D e ’1

0

where γ is the Gr¨ neisen constant, which under the Debye assumptions is

u

equivalent to γth introduced in (6.1.13).

The optical modes can be signi¬cant at high temperatures. The frequency

variation across the optical modes is limited. With the simple approximation of

a constant frequency ωE we can get a rough approximation to their contribution to

the speci¬c heat,

(˜E/T )2 e˜E /T

Copt = 3nNkB , (9.2.20)

v 2

e˜E /T ’ 1

where ˜E = hωE/kB.

9.3 Creep and rheology

The ultimate control on the behaviour of Earth materials under deformation comes

from processes at the atomic level, and the cumulative effects of the atomic level

interactions manifest themselves through the nature of the constitutive relation.

Many different mathematical forms for constitutive relations have been devised

to provide a macroscopic description of behaviour as different classes of atomic

processes become dominant (see, e.g., Evans & Kohlstedt, 1995). Such constitutive

relations have been deduced from experimental tests that, by necessity, are of

limited duration. In consequence the strain rates in the laboratory are rarely smaller

than 10’5 s’1 , which is still many orders of magnitude larger than the fastest natural

strain rates in rocks (around 10’9 s’1 ).

9.3.1 Crystal elasticity

The elastic properties of single crystals can be related to the statistical properties

of multi-atom assemblages through the same concept of strain energy as we have

used in the continuum description in Chapter 4, and in the treatment of the

potential energy of lattice vibrations in Section 9.2. The natural coordinates are

the positions of the atoms rather than their reference positions. Thus an Eulerian

description is appropriate and the Cauchy strain (2.2.27) needs to be used for

large deformation. This semi-classical treatment provides a good description of

behaviour in conjunction with the periodic nature of the atomic lattice. The

¬‚uctuations in the quantum phase vector in the periodic lattice depend on the

inverse elastic moduli for the material; these ¬‚uctuations are manifest through

9.3 Creep and rheology 163

reduction of the intensity in the X-ray diffraction peaks associated with different

lattice parameters (see, e.g., Chaikan & Lubensky, 1995). Careful X-ray intensity

measurements can therefore provide some constraints on elastic properties of the

lattice.

The in¬‚uence of crystal defects can be incorporated into the description of the

elastic state through the way that the strain energy is modi¬ed. The Eulerian stress

tensor in a material with defects can be expressed in terms of the Helmholtz free

energy F as

‚F ‚F

σij = ρ F ’ n δij + , (9.3.1)

‚n ‚A

A n

where n is the number of particles per unit volume and A is the displacement

gradient (5.1.1). The contribution hij = ‚F/‚A)n can be recognised as the normal

stress term, thermodynamically conjugate to strain, as in (6.2.11). The ¬rst term on

the right hand side of (9.3.1) allows for the presence of vacancies and interstitials:

‚F

= μ, (9.3.2)

‚n A

where μ is the chemical potential. The stress tensor can therefore change in an

isothermal state as a consequence of changes in either hij or μ. In the periodic

atomic lattice hij couples to the separation between layers, but μ is linked to

changes in density (Chaikan & Lubensky, 1995) and merely changes the isotropic

(pressure) component of stress.

9.3.2 Deformation behaviour

The deformation of materials becomes more pronounced as the melting

temperature TM is approached, and this dependence of the creep of a solid can

be represented in the general form (see, e.g., Evans & Kohlstedt, 1995)

m

σ d0 TM

n

™ =B exp ’g , (9.3.3)

G d T

with a dependence of the strain rate ™ on the grain size d, through the granularity

index m. g is a dimensionless quantity, and the stress σ is normalised by the shear

modulus G. If n = 1 we have a Newtonian relation, but generally n > 1 and we

have power-law creep; a typical value for n would be 3“5.

The apparent viscosity at constant stress for the general relation (9.3.3) is

therefore

σ Gn d m

σ TM

·σ = = exp g . (9.3.4)

™ (σ) Bσ d0 T

Alternatively the viscosity at constant strain rate is

1/n m/n

σ( ™ ) GB d TM

·™ = = exp g , (9.3.5)

d0 nT

™ ™™

164 From the Atomic Scale to the Continuum

which is proportional to σ(1’n)/n. The sensitivity to temperature will be lower at

constant strain rate for n > 1.

Creep deformation occurs through the transport of matter by the movement of

defects in the lattice. Diffusion controlled creep involves the motion of atoms

and vacancies in opposite directions and is Newtonian (n = 1). In contrast the

movement of dislocations through the lattice leads to non-Newtonian behaviour.

Diffusion creep

Under weak stress diffusion creep is a favoured mechanism. Different styles of

diffusion can be recognised depending on whether diffusion is faster through the

matrix or along the grain boundaries. Diffusion of the slowest ion along the fastest

path determines the nature of the diffusion and thus the granularity index m. For

matrix diffusion (Nabarro“Herring creep) m = 2, and for grain-boundary creep

(Coble creep) m = 3. Under compression the formation of vacancies in a lattice

will be suppressed, but under tension such vacancy formation will be enhanced.

There will therefore be a gradient in the concentration of vacancies in a solid under

stress leading to a ¬‚ow of vacancies from (9.1.7), with a consequent reverse ¬‚ow

of atoms. The strain rate is approximately inversely proportional to temperature, so

that the viscosity from (9.3.4) will be close to a linear function of T . The activation

energy for this class of diffusion is just that for self-diffusion in the lattice (Poirier,

1991).

When the diffusion through the lattice is very slow, the dominant mode of

diffusion will be the migration of atomic species along the grain boundaries. In

such Coble creep the viscosity is again approximately linear in T , and in order for

grains to remain in contact there needs to be some grain-boundary sliding.

Dislocation creep

For higher stresses, the deformation of a crystal is controlled by the movement

and interactions of dislocations and the thermal effects become more pronounced.

(a) Edge (b) Screw

kink

jog

E

n n

glide

A

plane

b b

l l

Figure 9.3. The migration of dislocations is controlled by kinks and jogs. (a) Migration of

an edge dislocation in the glide plane by the migration of kinks does not involve diffusion

processes. Jog motion out of the glide plane requires a diffusive component, with the

absorption (A) or emission (E) of vacancies. (b) Screw dislocations can slip on any plane

containing the line direction.

9.3 Creep and rheology 165

Dislocations move by the energy-ef¬cient process of the progressive breaking and

reforming of atomic bonds, and are the carriers of plastic deformation through a

solid.

The effect of deformation on a crystal is to force migration of existing

dislocations, that then interact with possible annihilation of the defects. For larger

strains, dislocations will be created and migrate leading to complex dislocation

tangles. Each of the processes of creation, migration and annihilation requires

work from the applied stresses.

The net local stress acting on a dislocation is the superposition of the externally

applied stresses and internal stresses due to the in¬‚uence of other dislocations,

defects and interfaces. The internal stress will normally resist the migration of

the dislocation and may increase with time, the process of hardening, or decrease

over time (recovery) depending on the way that the microstructure is modi¬ed.

There are two basic types of dislocation movement, glide in which the dislocation

moves in the surface de¬ned by its line and the Burgers vector, and climb in which

the dislocation moves out of the glide surface. For a density of mobile dislocations

ρd of a similar type with average dislocation velocity v the creep rate is

¯

™ = ρdb¯,

v (9.3.6)

where b is the length of the Burgers vector. The dislocation velocity is related

linearly to the deviatoric stress ”σ

v = k ”σ exp[’Q/RT ],

¯ (9.3.7)

with strong temperature dependence because of thermal activation barriers. The

hardening of materials is reduced at higher temperatures because dislocation

pile-ups can be avoided.

In silicate minerals, only small portions of the dislocation, called kinks and jogs,

are mobile at any moment. Both represent abrupt changes in the line direction;

kinks lie within the glide plane, and jogs take the dislocation locally out of the glide

plane. Thus dislocation migration in the glide plane occurs by kink migration, and

does not require atomic diffusion. Whereas climb out of the glide plane through the

movement of jogs needs atomic diffusion into the extra half plane, either along the

dislocation core or through the lattice (Figure 9.3a). The screw component can slip

along any glide plane containing the line direction (Figure 9.3b). The migration

of a kink has to overcome the resistance of the crystal structure, since there is an

energy barrier to be overcome to transfer the dislocation position from one location

to another. However, this energy barrier is much less than that required to force

atoms past each other so that the shear resistance is orders of magnitude less than

it would be in a perfect crystal. The mobility and the number of jogs and kinks will

be in¬‚uenced by the concentration of point defects in the crystal.

The initial application of signi¬cant stress can produce signi¬cant transient

effects through the migration of dislocations to grain boundaries. Steady-state

rheology is governed by the rate-limiting step in the sequence of generation,

166 From the Atomic Scale to the Continuum

movement and elimination of the dislocations that have a strong dependence

on the local crystal environment. Generation of dislocations may arise at

boundaries, intersections of other dislocations or localised defects. Dislocations

may be annihilated by reaction with dislocations of opposite Burgers vector or by

interaction with a boundary in the process of recrystallisation. In the steady state

the dislocation density ρd is proportional to the square of the deviatoric stress

ρd = C(”σ)2, (9.3.8)

where C is nearly independent of temperature. Thus for a single style of

dislocations we can combine (9.3.6), (9.3.7) and (9.3.9) to obtain a simple

power-law creep model

™ = Ckb(”σ)3 exp[’Q/RT ], (9.3.9)

with an exponent n = 3. Evans & Kohlstedt (1995) present a representative sub-set

of constitutive equations for dislocation-controlled creep under a variety of other

scenarios.

9.4 Material properties at high temperatures and pressures

The conditions prevailing in the deep interior of the Earth lie far from standard

laboratory conditions. Although very signi¬cant progress has been made in

achieving very high pressures, it is particularly dif¬cult to achieve high pressure

and calibrated uniform temperatures across a specimen.

The pressure and temperature regime that needs to be explored to study the

Earth™s interior is illustrated in ¬gure 9.4. The pressure variation, here for model

AK135 (Kennett et al., 1995), depends on the radial density distribution as shown

in (6.1.2), but is not very sensitive to the choice of Earth model. However, the

temperature distribution is much less certain and estimates of temperature at great

depth have considerable uncertainty (1000 K or more), particularly in the core.

The challenge is therefore to achieve pressures of 150 GPa and temperatures

approaching 4000 K so that the full range of conditions in Earth™s mantle can

be investigated. Pressure represents the force applied per unit area, and thus

ultra-high pressures require either or both of very large force and a very small area

of application.

9.4.1 Shock-wave techniques

The ¬rst technique to achieve pressures comparable to those at the centre of the

Earth (370 GPa) was the use of shock waves induced by explosions or the impact of

high velocity projectiles. The impact or explosion produces a compressional shock

front that passes through the material faster than the local sound wavespeed, until it

is re¬‚ected by the far side of the target to return as a rarefaction. Shock waves heat

as well as compress the target material. Temperatures can reach higher than 2000

9.4 Material properties at high temperatures and pressures 167

5000 400

4000

300

Temperature [oC]

Pressure [GPa]

3000

T

200

p

2000

Lower Mantle Outer Core Inner Core

100

1000

MAP DIAMOND CELL

0 0

0 1000 2000 3000 4000 5000 6000

Depth [km]

Figure 9.4. The pressure variation from the AK135 model (Kennett et al., 1995) and the

likely temperature band for the interior of the Earth as a function of depth. The regimes

accessible to multi-anvil (MAP) and diamond cells with laser heating are also indicated.

K for silicate minerals, and melting can sometimes occur under shock conditions.

The experiments measure the shock-wave velocity U and the particle velocity u,

imparted to the particles of the material, and the temperature. A common means

of achieving high incident projectile velocities is the use of a two-stage light-gas

gun in which a large primary projectile imparts its momentum to a small ¬‚yer plate

that impinges on the target. The experimental conditions can be varied by adjusting

the projectile velocity, but only one shock state can be investigated at a time. As

a result a very substantial investment is required to map out the relation between

shock-induced pressure and density known as Hugoniot curves.

The passage of a shock-front changes the state of the material and the properties

in front of the shock (a) and after its passing (b) are related by the conservation of

mass, momentum, internal energy (or enthalpy) across the pressure discontinuity

of the shock-front. The density ρ, pressure p and speci¬c internal energy e and

enthalpy h are thus related by the Rankine“Hugoniot conservation equations:

ρb = ρa(U ’ ua)/(U ’ ub), (9.4.1)

pb = pa + ρa(U ’ ua)(ub ’ ua), (9.4.2)

ea = eb + 1 (pa + pb)/(1/ρa ’ 1/ρb), (9.4.3)

2

ha = hb + 1 (pa ’ pb)(1/ρa + 1/ρb). (9.4.4)

2

168 From the Atomic Scale to the Continuum

If an arti¬cially porous sample with ρa considerably lower than the normal crystal

density is employed as the target, large increases in internal energy can be achieved

at a given pressure.

The initial effect of the shock is to produce uniaxial compression, but at high

shock-induced pressures there will be bifurcation of the shock wave with two

successive shocks. The ¬rst achieves the uniaxial compression to the Hugoniot

elastic limit and the second slower shock wave brings the material to its ¬nal

state. For material with phase transitions there will be further shock bifurcations to

produce multiple shock-fronts (Ahrens, 1987).

The Hugoniot trajectory in shock pressure-density space needs to be mapped to

isothermal (constant T ) or isentropic (adiabatic) conditions for comparisons with

other classes of information (see Poirier, 1991, §4.7.3). The intermediary is the

Mie“Gr¨ neisen equation

u

pH ’ pK = γρ(eH ’ eK), (9.4.5)

which relates the pressures and energies along the Hugoniot (pH, eH) to another

thermodynamic state for the same density. The Gr¨ neisen ratio γ is given by

u

1 ‚p ±th KS

γ= = , (9.4.6)

ρ ‚E ρCp

ρ

where, as in (6.1.15), ±th is the coef¬cient of thermal expansion, KS is the adiabatic

bulk modulus and Cp is the speci¬c heat at constant pressure. The ratio γ is

frequently assumed to depend only on density.

Ahrens & Johnson (1995a,b) provide a summary of shock-wave data for many

classes of rocks and minerals. The properties of silicate minerals at high pressure

can be estimated quite well from those of mixtures of the relevant oxides.

Shock-wave experiments provide direct calibration of pressure and such results are

therefore important for pressure measurements in more indirect methods.

9.4.2 Pressure concentration by reduction of area

The simplest device to subject samples to pressure is a piston cylinder apparatus in

which uniaxial load is applied directly to relatively large specimens. In this way

pressures of the order of 7 GPa can be achieved, comparable to conditions at 100

km depth.

A more uniform loading and much higher pressures can be achieved in a

multi-anvil press. An externally applied uniaxial load is transmitted to a small

octahedral sample assembly through a secondary set of anvils. Six anvil wedges

press on the faces of eight cubes, typically made of tungsten carbide for hardness.

These cubes squeeze the octahedral sample assembly nestling in the inner space

created by truncating the cube corners. A ceramic octahedron acts as the

pressure-transmitting medium. With large presses (nominal ratings up to 5000

tons) pressures approaching 30 GPa can be reached for sample volumes of several

9.4 Material properties at high temperatures and pressures 169

cubic millimetres. This means that the pressure conditions throughout the mantle

transition zone to about 800 km are accessible; mantle temperatures can also be

matched. There is a trade-off between the size of the sample and the pressure that

can be achieved. The largest presses tend to be used for synthesising signi¬cant

quantities of high-pressure materials or in-situ measurements such as thermal and

electrical conductivity.

Studies with the multi-anvil press have played a major role in elucidating the

phase changes in the mineral assemblage through the upper mantle and the mantle

transition zone (Figure 9.5). Recent developments have included establishing

multi-anvil presses at synchrotrons so that the high-intensity X-ray beams can be

used to examine the properties of the material during deformation. Thus, e.g.,

the density can be monitored across a phase transition from the change in unit

cell size. The elastic properties of unquenchable high-pressure phases that cannot

be brought back to normal conditions have been measured with a combination

of synchrotron X-ray methods and ultrasonic interferometry. With a suitably cut

transducer both P and S waves can be injected into the specimen through the sample

holder. The change of length of the specimen can be monitored during the course

of the experiments and hence the timing of the multiply re¬‚ected waves can be

unravelled to provide wavespeed information.

Depth [km]

MgSiO3

100 200 300 400 500

1500

X

Ortho-enstatite

β+

Temperature [˚C]

L St

High-P Clino-enstatite

Continental Shield Fe

1000 Geotherm [unquenchable]

Spinel +

Stishovite

Low-P Clino-enstatite

500

5 10 15

Pressure [GPa]

Figure 9.5. The pressure and temperature (P,T ) phase diagram for the MgSiO3 composi-

tion. In-situ measurements have revealed the elastic properties of the unquenchable phase

high-pressure clino-enstatite and the small discontinuity in elastic properties from ortho-

enstatite is compatible with the seismic X discontinuity beneath continental shields, but

not the shallower Lehmann discontinuity (L). The movement of the phase boundary with

increasing Fe content is indicated. [Courtesy of J. Kung.]

170 From the Atomic Scale to the Continuum

Higher pressures can be achieved in diamond anvil cells in which two anvils

made from superhard materials squeeze a very small sample (10“500 μm in

diameter). In this way, with relatively modest forces, it is possible to attain

pressures up to several hundred GPa. Heating can be provided externally (up to

1000 K) or internally by laser heating. The diamond anvils are transparent to a

wide range of the electromagnetic spectrum (X-rays, visible and infrared radiation)

and so samples can be monitored during an experiment. The laser heating can

be induced by directing beams from high-intensity infrared lasers through the

anvils with a focus at the specimen. Temperatures can be determined from the

local Planck emission spectrum with a 1 μm resolution. It is dif¬cult to avoid

temperature gradients across the specimen away from the directly heated region,

but now uniformity can be achieved for samples as large as 100 μm. Fortunately

the diamonds themselves stay cool because of their high thermal conductivity.

The pressure in the diamond cell is determined indirectly from the shift in the

¬‚uorescence line emitted by tiny ruby chips embedded in the sample. The pressure

scale is calibrated by making simultaneous measurements of the density of metals

in the cell, using in-situ X-ray diffraction measurements, and the shift of the ruby

¬‚uorescence line. The pressure is derived from the density through the equation

of state for the metal derived from shock-wave experiments, with corrections for

thermal pressure (typically a few GPa).

Although some materials survive the transition from the diamond anvil cell

to laboratory conditions without loss of structure or change in chemistry, many

high-pressure states are metastable and need to be characterised in their high

pressure and temperature state. The relatively small size of the diamond

cell assembly means that rapid X-ray diffraction studies can be made at high

temperatures using a synchrotron beam. The geometry and intensity of the recorded

diffraction patterns provide information on the parameters of the periodic cell

structure and also on the actual charge density, which is valuable information for

computational modelling. The pressure“density relationships can then be used to

infer parameters such as the bulk modulus and its pressure derivative.

Despite the small size of the samples in the diamond anvil cell, it is possible

to make direct measurements of elastic wavespeeds for single crystals. Brillouin

spectroscopy detects the frequency shift of electromagnetic waves induced by

thermally excited phonons (vibration modes of the crystal). These shifts can be

directly related to elastic wavespeeds. It is also possible to carry out ultrasonic

interferometry on samples in the diamond cell using waves delivered to the

specimen through one of the anvils (Jacobsen et al., 2002). From the interference

pattern of re¬‚ected waves the travel times of P and S waves can be determined, and

then converted to wavespeed with knowledge of the sample size under pressure

(Figure 9.6). In either case the full set of anisotropic elastic moduli can be

determined from multiple experiments for different crystal orientations. Pressure

derivatives can also be found using different experimental conditions.

9.5 Computational methods 171

c11

300

Elastic constant [GPa]

200

c12

100

c44

0 2 4 6 8 10

Pressure [GPa]

Figure 9.6. The variation of the three independent elastic moduli of cubic ferropericlase

(Mg,Fe)O with pressure to 10 GPa determined by ultrasonic interferometry. Ferropericlase

coexists with silicate perovskite (Mg,Fe)SiO3 in the lower mantle.

Pressures beyond the range of the conventional multi-anvil press (80 GPa)

have been achieved with a three stage system based on the use of a superhard

nanodiamond aggregate (T. Irifune - personal communication). The same material

offers the possibility of enlarged diamond anvils. In each case this means that larger

samples can be investigated at high pressures.

9.5 Computational methods

Since the conditions deep in the Earth are dif¬cult to reach with experimental

methods, a useful alternative source of information comes from computational

methods applied at the atomic level. There are two broad classes of method

depending on the treatment of electronic structure. In electronic structure

calculations the Schr¨ dinger equations for the quantum behaviour of a set of atoms

o

are solved at some degree of approximations to provide direct access to material

properties. On the other hand the explicit solution of the Schr¨ dinger equations

o

can be avoided by using interatomic potentials that represent the energy of the

system through nuclear coordinates.

9.5.1 Electronic structure calculations

Such ¬rst-principles methods make no initial assumptions regarding the physical

nature of a mineral system, instead they rely on basic quantum mechanics and the

172 From the Atomic Scale to the Continuum

atomic number of an atom. The basic building blocks are nuclei and electrons. In

principle, if we place magnesium, silicon and oxygen atoms in a suitable box in the

proportions two magnesium atoms to one silicon atom and four oxygen atoms, with

the temperature and pressure set to room conditions, it should be possible to obtain

the magnesio-silicate mineral forsterite as a ¬nal state. However, this is currently

too complex a problem to be solved directly so some level of information needs to

be supplied about the likely structure.

The energy of a quantum mechanical system of electrons and atoms can be found

from the solution of the Schr¨ dinger equations for the electron states

o

Hψ(r1, r2, · · · , tN) = Eψ(r1, r2, · · · , tN), (9.5.1)

where r is the coordinate of the ith electron, H is the Hamiltonian and E is the

energy. The Hamiltonian operator consists of a sum of three terms, the in¬‚uence of

the atomic nuclei (Ven), interactions between electrons (Vee) and the kinetic energy

of the system, so that

N N

1

∇2.

1

H = Ven + ’ (9.5.2)

i

2

|ri ’ rj|

i<j i=1

The external potential Ven arises from the interaction of electrons with the M

atomic nuclei,

M N

Zp

Ven = ’ , (9.5.3)

|ri ’ Rp|

p i

where the charge on the nucleus at Rp is Zp. In addition the in¬‚uence of electronic

spin needs to be considered explicitly in some circumstances, such as the presence

of a magnetic ¬eld.

The coupled Schr¨ dinger equations for the N electrons (9.5.1) cannot be solved

o

directly because of the complexity and large dimension of the equation. The

representation of the 3N-dimensional wavefunction ψ scales as gN, where g is

the number of degrees of freedom for a single electron. Nevertheless, there are a

number of approximations based on the combination of self-consistent solutions

for a single-electron Schr¨ dinger equation with summation over all electrons in the

o

system. The best approximations using correlated methods scale as N7 and can be

used to determine chemical properties of small molecules (Harrison, 2003) but are

not feasible for larger systems.

Fortunately it is not necessary to solve for the 3N-dimensional wavefunction

ψ in order to extract information on the energetics of the system. The electron

density ρe(r) uniquely determines the positions and charges of the nuclei and

thus the Hamiltonian operator H (9.5.2). The energy functional contains three

elements, the kinetic energy T , the electrostatic potential W (from both nuclei

9.5 Computational methods 173

and electronic charge density) and a residual contribution Exc associated with

correlations between electrons and exchange effects associated with electron spin:

E[ρe] = T [ρe] + W[ρe] + Exc[ρe]. (9.5.4)

The density functional approach of Kohn and Sham exploits this property to break

the N-body Schr¨ dinger equation into a set of N coupled equations that closely

o

resemble those for a single electron. The Kohn“Sham equations for a set of orbitals

{φi} are

12

2∇ + Ven + Vee[ρe(r)] + Vxc[ρe(r)] φi = iφi, (9.5.5)

where the electron density ρe(r) and energy E are given by

N N

2

ρe(r) = |φi(r)| , E[ρe] = i. (9.5.6)

i i

The exchange correlation potential Vxc is the functional derivative of Exc with

respect to ρe, i.e., Vxc = ‚Exc/‚ρe.

The set of non-linear equations (9.5.6) describes the behaviour of non-interacting

˜electrons™ in an effective local potential that actually has an unknown component

Vxc. If the exact functional were known, then there would be a complete

correspondence between the many-body situations and the uncoupled equations.

For any imperfect representation of the density functional the energy estimate E[ρe]

will be an upper bound on the true ground state energy E0. The development of

excellent approximations to the density functional allows unbiased and predictive

studies of a wide range of materials, and they are often referred to as ab initio

methods.

The exact form of the exchange-correlation potential is not known, but two

styles of approximations have proved to be quite effective (see, e.g., Harrison,

2003). The local density approximation (LDA) de¬nes the exchange correlation

as a function of the electron density at a point, and works well for most

insulating systems. However, the computation of energy differences between rather

different structures can have signi¬cant errors; thus the binding energy tends to

be overestimated and energy barriers in diffusion are likely to be too small. A

more sophisticated approximation, the generalised gradient approximation (GGA)

includes information from the gradient of the electron density as well, which is very

effective for metals but not for strongly correlated oxides such as FeO or Fe3 O4 .

In our sketch of the theory we have not included the in¬‚uence of electron

spin, but both spin up and spin down states need to be considered separately

for understanding magnetic effects that are important for minerals containing

transition metals. The density functional theory as described above is a ground state

approach, based on an athermal state at absolute zero Kelvin without any zero-point

energy. The absence of thermal excitations means that perfect periodicity can

be assumed for the lattice of a crystal, and hence an in¬nite medium can be

174 From the Atomic Scale to the Continuum

simulated with a single primitive unit cell subject to periodic boundary conditions.

For the orthorhombic symmetry of forsterite we need to employ four formula

units and thus 28 atoms. Experimental work cannot attain direct correspondence

with this athermal state, but in many cases room-temperature thermal excitation

of the electrons is small and electronic band structure and material properties

can be compared between experiment and the calculations from ¬rst principles.

Temperature effects tend to be more pronounced for magnetic and anisotropic

systems.

The ef¬cient description of the φi is an important part of the implementation

of density functional theory, and these are usually assembled from a superposition

of basis functions. Many approaches use plane waves that would correspond to

free electron states; since solutions can be improved through the addition of further

plane waves.

Only a fraction of the full set of electrons is involved in chemical bonding

and hence it is possible to develop a scheme of pseudopotentials in which the

interactions of valence electrons and nuclei are represented by a weak effective

potential. The pseudopotential is much smoother than the bare core potential

around a nucleus and so not only are fewer electrons considered, but also fewer

basis functions are needed. Eliminating the tightly bound electrons in the system

can make signi¬cant reductions in the number of effective φi, e.g., a drop from

280 to 48 electrons for forsterite. Careful use of pseudopotentials cuts down on the

computational effort, so that large mineral systems such as clays can be studied, or

alternatively supercells with hundreds of atoms can be used to describe scenarios

without full periodicity, as in thermal phenomena or crystals containing defects.

9.5.2 Atomistic simulations

The potential energy of an assembly of N atoms U(r1, r2, . . . , tN) is expressed as a

function of the position vectors {rj} from the various nuclei. In this approximation

the electronic structure of the system is subsumed into the potential function. U is

commonly expanded into a series of different levels of interatomic interaction

N N N N N

1 1

U= φij(ri, rj) + φijk(ri, rj, rk) + · · · , (9.5.7)

2 3

i=1 j=1 i=1 j=1 k=1

where the symbol indicates that multiple counting of terms is avoided. The two-

body functions φij depend only on the positions of pairs of atoms (i, j), whereas

the three-body terms depend on the atomic triplet (i, j, k). The two-body term φij

can be represented in terms of an electrostatic component and a residual potential

Vij that has both attractive and repulsive components,

qiqj

φij(ri, rj) = + Vij(|ri ’ rj|), (9.5.8)

|ri ’ rj|

9.5 Computational methods 175

where qi is the charge on the ith atom (or ion). Analytic functions are used to

represent Vij such as the Buckingham potential, suitable for ionic solids,

V(r) = Ae’r/ρ ’ Cr’6; (9.5.9)

the repulsive component A exp(’r/ρ) represents the overlap of atomic energy

shells, and the longer range attractive term Cr’6 includes induced dipole effects.

The parameters A and C are empirical and need to be adjusted to observed

crystal properties. Many-body effects can be included by incorporating, e.g.,

angle-dependent forces or triple dipole terms (Catlow, 2003).

Methods which rely on parametric ¬ts are valid for the range of physical

parameters for which the ¬t is made and extrapolation outside this range has to

be undertaken with caution. Thus, for example, the presence of a high-pressure

phase transition could only be predicted with con¬dence if the physical parameters

for both the phases were included in the ¬tting procedure.

9.5.3 Simulation of crystal structures

The basic tool kit of ab initio methods allow the construction of lattice energy

and electronic charge density for a given set of atomic positions and species. If

spin-polarisation effects are included then magnetic moments can also be found.

Such calculations have demonstrated that the high-temperature polymorph of iron,

with a face-centred cubic structure, has a complex magnetic ground state that gives

rise to anomalous thermal expansion.

For understanding Earth materials further properties are desirable, including

the equilibrium volume, which require the forces acting on the atoms and the

15 8.0

P S

100 GPa

14 7.5

Wavespeed [km/s]

Wavespeed [km/s]

100 GPa

13 7.0

12 6.5

11 6.0

0 GPa

0 GPa

10 5.5

100 010 001 100 111 001 010 001 100 010 001 100 111 001 010 001

Propagation direction Propagation direction

Figure 9.7. Variation of (a) longitudinal and (b) shear wave speeds of (Mg,Fe)(Si,Al)O3

perovskite as a function of propagation direction and pressure calculated using ab initio

methods [after Li et al. (2005)].

176 From the Atomic Scale to the Continuum

stress on the cell. These properties can be extracted from the ¬rst derivatives

of the energy. The ¬rst derivative with respect to atomic positions provides the

forces, and the tensor derivative with respect to small strain yields the stress.

The equilibrium volume is best determined by ¬nding the minimum of the

energy“volume relationship for variable cell volume.

The computation of the stress for a strained con¬guration enables the direct

evaluation of the corresponding elastic constants (Kiefer et al., 2001). Alternatively

the elastic constants can be derived from the variation of strain-energy density with

lattice strains (Steinle-Neumann et al., 1999).

An example of the use of density functional theory with a plane wave basis set,

using the generalised gradient approximation, is provided by the work of Li et

al. (2005) who have examined the elasticity of ferro-magnesium perovskite with

allowance for aluminium substitution for silicon. This is the most abundant mineral

on Earth through the lower mantle, but there is little experimental constraint on

properties. The elastic properties are illustrated in Figure 9.7 as a function of

propagation direction and pressure for one con¬guration of aluminium substitution.

An extensive set of calculations suggests that the spin state of the iron atoms has a

smaller in¬‚uence on the elastic properties than the inclusion of aluminium instead

of iron.

9.5.4 Finite temperature

An important aspect of material properties is the in¬‚uence of temperature and in

particular the temperature derivatives of material properties such as elastic moduli

and wavespeeds. The methods for ¬nite temperature are based on thermodynamic

properties in terms of the Helmholtz free energy derived from statistical mechanics.

(a) Lattice dynamics

This is a semi-classical approach that uses a representation of the unit cell in terms

of independent quantized harmonic oscillators whose frequencies change with cell

volume (Born & Huang, 1954) and so thermal expansion can be described. In

this quasiharmonic approximation the motions of the atoms are treated collectively

as lattice vibrations (phonons) whose frequencies ωq for wave vector q are

determined by solving

mω2(q)sk(q) = Dkl(q)sl(q), (9.5.10)

for lattice displacements s(q), where the dynamical matrix D(q) is de¬ned as

‚2E

D(q) = exp(iq · rij). (9.5.11)

‚si‚sj

ij

Here rij is the separation of atoms i and j, and ui, uj are their displacements

from the equilibrium position. For a unit cell with n atoms, there are 3n sets

of eigenvalues ω2(q) for a given wave vector q associated with 3n sets of

eigenvectors s(q) that describe the patterns of atomic displacement for the mode.

9.5 Computational methods 177

The vibrational frequencies can be calculated from ¬rst principles by looking at the

in¬‚uence of small displacements (Kresse et al., 1995).

Once the eigenfrequencies have been determined, the thermodynamical

properties can be found by using standard methods from statistical mechanics for

Bose“Einstein oscillators. The Helmholtz free energy, as a function of volume,

takes the form

2

hωq

1 ”V

F= + kBT ln 2 sinh , (9.5.12)

2K V 2kBT

q

where kB is the Boltzmann constant. The ¬rst term on the right-hand side of

(9.5.12) is the potential energy associated with the compressibility K of the solid as

an elastic continuum. The second term is the sum of the free energies in the lattice

modes.

If we make the simplifying assumption that a change of volume ”V produces

the same relative change of frequency for each mode

”ω ”V

=γ , (9.5.13)

ω V

then the condition for minimum free energy from (9.5.12) is

hωq

1 ”V ¯

γhωq 1 coth

= = γE(T ). (9.5.14)

2

K V 2kBT

q

¯

Here E(T ) is the energy in the lattice modes at temperature T and the Gr¨ neisen

u

constant γ relates the dilatation to the mean thermal energy density

”V ¯

= KγE(T ). (9.5.15)

V

The quasiharmonic approximation assumes that the modes can be treated

as independent. However, as the vibrational amplitudes increase at higher

temperatures phonon“phonon interactions become important and the independence

of the different lattice modes is destroyed. Up to the Debye temperature ˜D at

which all vibrational modes are excited, the assumption that the frequency of the

lattice modes is just a function of volume is reasonable. At ambient pressure ˜D

is in the range of a few hundred K for metals and 600“1000 K for typical silicate

minerals.

(b) Molecular dynamics

An alternative approach to the inclusion of the effects of temperature is provided

by the use of molecular dynamics techniques. This approach explores the time

evolution of one representation of a crystal system. The ergodic hypothesis of

statistical thermodynamics is then used to equate the time average of a physical

¯

property X to the ensemble average over con¬gurations X .

178 From the Atomic Scale to the Continuum

For the n atoms in the con¬guration, the time variation is found by integrating

the set of n coupled equations

‚2ri

mi 2 = fi(r1, r2, . . . , rn), i = 1, . . . , n , (9.5.16)

‚t

for the positions of the atoms {ri} controlled by the atomic forces f and how

these are speci¬ed. The numerical solution is normally achieved with periodic

boundary conditions on a supercell; 2 — 2 — 2 unit cells is usually suf¬cient for

the computation of thermodynamic variables. The fastest thermal vibration of the

lattice needs to be covered in a few time steps, so typical steps are of the order

of femto-seconds (10’15 s) with a total simulation duration of a few pico-seconds

(10’12 s). The initial conditions are the positions of the atoms and their velocities

drawn randomly from a Gaussian distribution corresponding to the temperature of

interest. To encourage the equation system towards convergence during the interval

of numerical integration, the temperature can be temporarily increased to overcome

energy barriers.

The equations of motion (9.5.16) conserve energy, and so the natural

thermodynamic ensemble is the microcanonical where the internal energy E is

¬xed. However, the normal experimental state has temperature controlled. To

undertake molecular dynamics simulations in the canonical ensemble, with n, V, T

¬xed, a thermostat has to be applied to the system, e.g., via friction terms applied

to the equations of motion (Allen & Tildesley, 1989). Free energies cannot be

obtained directly from molecular dynamics simulations, but the difference between

the free energy and that of a reference system (”F) can be found by thermodynamic

integration. ”F is the work done in a reversible and isothermal change from

the total energy of the reference system to that with the actual total energy; the

calculations need to be performed with care (de Wijs et al., 1998).

9.5.5 In¬‚uence of defects

As we have seen in Sections 9.1, 9.3 the presence of defects has a profound

in¬‚uence on the properties of materials and so there is need for computational

methods to include the presence of such defects. Two styles of approach are

currently employed: (a) the construction of a periodic array of defects, (b) the

defect is embedded in an in¬nite representation of the surrounding crystal.

(a) Supercell methods

The simplest procedure for calculating the energies of defects is to build a periodic

construct with a defect at the centre of a large ˜supercell™. The in¬‚uence of the

defect can then be determined by comparing the results of calculations for the

perfect lattice with those when the defect is present. The difference in energy

between the two states is the energy of formation for the defect, and calculations

may be performed for either constant volume or constant pressure.

Rather large supercells are needed if complex defects are being considered. Even

9.5 Computational methods 179

Region IIa

Region I

Region IIb

Defect extends to infinity

Figure 9.8. Schematic representation of the Mott-Littleton technique for treating the in¬‚u-

ence of defects, with an inner zone (Region 1) surrounding the defect, a transition zone

(Region IIa) and a perfect crystal extending to in¬nity (Region IIb). Only one octant is

shown.

with such large supercells there will be a contribution from defect interaction

as well as formation. However, this contribution can be turned to advantage

by calculating defect energy for different sizes of supercells to determine defect

interactions as a function of defect spacing.

A further complication that can arise when investigating different possibilities

for the presence of a class of defect in a crystal structure is the number of different

con¬gurations that have to be tested to ¬nd the structure with the smallest energy

of formation and hence highest likelihood.

(b) Embedded defects

The alternative approach is based on ideas introduced by Mott and Littleton, in

which an isolated defect or defect cluster is placed at the centre of an idealised

crystal structure extending to in¬nity (Figure 9.8).

In region I, immediately surrounding the defect, with 100“300 atoms a

minimisation procedure is used to bring the atomic forces to zero. In the interface

region IIa, the displacements are evaluated as the sum of the response to individual

defects including the polarisation effects of defect charge, together with short-range

interactions between the ions in region I and those in this interface zone. In the

outermost region IIb only the polarisation due to the net charge of the defect

con¬guration is considered.

10

Geological Deformation

In Part I we have concentrated on materials which in the aggregate, at least,

have relatively simple properties and focussed on the impact of the latest phase

of deformation. These concepts are of considerable utility in understanding the

nature of evolution of structure in rocks, but we have to take into consideration the

nature of the complex mineralogical composites whose properties depend strongly

on physical conditions.

We have abundant evidence of past deformation in rocks on many length scales.

The structures that are seen both in the ¬ne scale microfabric and on much larger

scales in ¬eld exposures frequently re¬‚ect the in¬‚uence of multiple phases of

deformation and spatially varying strain.

Different classes of rocks have a wide range of chemical and mineralogical

compositions, which have experienced physical and chemical environments that

themselves may have shaped the actual fabric. Rock masses contain materials with

a large variation in strength, determined by local mineralogy, which respond in very

different ways to the same stress ¬eld, leading to differences in local strain. Further,

the speci¬c temperature and pressure conditions due, e.g., to depth of burial will

affect the way in which the rock units behave.

The action of the stress system will lead to a pattern of deformation that will vary

from point to point as parts of the body move relative to each other, in part because

of the differences in mechanical properties of the different parts of the rock mass.

One of the dif¬culties that has to be faced in developing descriptions of

deformation is that most of the ¬‚ow laws which have been developed for geological

materials are based on low strain experiments, yet normally with a much higher

strain rate than in nature. It is now possible to examine high-strain conditions in

torsion, but still at high strain rate. The link from the laboratory to the conditions

prevailing in a major orogeny involves substantial extrapolation, in particular with

regard to rates of deformation.

180

10.1 Microfabrics 181

Figure 10.1. Example of dislocation structures responsible for plastic deformation in a

quartz foil. The width of the transmission electron microscope image is approximately

5 μm. Clear screw dislocations are seen at the bottom of the image in contrast to the

dominantly edge dislocations above. [Courtesy of J. Fitz Gerald.]

10.1 Microfabrics

10.1.1 Crystal defects

The crystal defect systems that we introduced in the previous chapter play an

important role in the microstructure of rocks through the way in which their

arrangement is in¬‚uenced by the thermal history.

For silicate minerals a particularly important class of point defects comes from

the diffusion of (OH) ions through the structure from the presence of water. The

Si“O“Si bridge bonds are replaced by Si“OH“HO“Si, so that the silicon“oxygen

Si“O bonding is replaced by the weaker OH“HO bond. The result is a signi¬cant

decrease in the mechanical strength of silicates, a greater ease of recrystallisation

to form new grains, and an enhanced rate of chemical reactivity.

The dislocation line defects have the most direct in¬‚uence on the mechanical

behaviour of crystals. Figure 10.1 shows part of a transmission electron micrograph

through a foil prepared from a quartz crystal deformed at 800—¦ C, and a strain rate

near 10’5 s’1 in creep deformation under a differential stress of 435 MPa. The

dislocations show points of rapid growth that are dragging the edge dislocations

behind them. Slip on crystallographic planes allows one part of the crystal to

undergo shear with respect to a neighbouring part, so that the shape of a grain

can change. Quite high densities of dislocations can occur in silicate minerals with

182 Geological Deformation

deformed materials reaching densities as high as 1016 m’2 at which points tangles

of dislocations largely neutralise slip. The slip systems operating within a crystal

are determined by energetic considerations and the ease of movement on a system

will depend on temperature and strain rate.

Within a crystal grain there are often subgrains with subtle differences in lattice

orientation, less than 5—¦ , which can be formed when dislocations group to form

planar walls that can move as a unit. After deformation the crystal will try to

reduce the dislocation density and arrays of dislocations can get organised into

planar networks.

10.1.2 Development of microstructure

The deformation of rocks frequently occurs at temperatures such that

recrystallisation can occur at the same time as the internal changes within the

crystal grains. The microstructures found in rocks are therefore a function of strain

rate, the temperature during the deformation process and, for silicates, the presence

or absence of water.

The ¬nal microstructure will be a function of the thermal history of the rock

mass. For example, we expect that rocks deformed at high strain rates or lower

temperatures will show strong evidence of the deformation with grains ¬‚attened

by the deformation, and internal structure within the grains. However, if the rock

is subsequently maintained at moderate temperatures there will be a process of

recovery as the internal energy of the crystal grains is reduced by the reduction

of dislocation density, with diffusion allowing dislocation climb to take place. In

addition, recrystallisation will occur whereby new crystal grains are developed with

internal strains removed. The orientation of these new grains will be signi¬cantly

different from that of the strain, so that high-angle boundaries will tend to form and

migrate through the material, removing most of the dislocations remaining from the

recovery stage. There may then be a process of grain growth as internal energy is

again reduced by the elimination of grain boundaries, as large grains absorb smaller

ones.

When deformation occurs at higher temperatures or lower strain rates, solid-state

diffusion can become important at the grain scale and dislocations can then climb

under the in¬‚uence of the stress regime causing the deformation. The result is

the migration of grain boundaries and possible growth of recrystallised zones and

the ¬nal microstructure will depend on the relative importance of a variety of

competing processes.

The majority of experimental apparatus for studying rock-deformation relies on

some class of compression through, e.g., a piston cylinder system. Such equipment

is of considerable value for many scenarios, particularly the in¬‚uence of pressure

and temperature on a rock, or the in¬‚uence of ¬‚uids. However, many natural

features are associated with very large shear strains that cannot be produced in

such a con¬guration. Experimental techniques have therefore been developed in

10.1 Microfabrics 183

A starting material c(0001)

E

γ =0

J ˜ 1.0

γ ˜ 1.4

B F

γ ˜ 0.5

J = 1.8

γ ˜ 3.4

C G

γ˜3

J = 3.8

D γ ˜ 13.5 H

γ ˜ 10

J = 3.3

Figure 10.2. Development of shear fabric in calcite (Carrara marble) subjected to large

¬nite torque at 500—¦ C, 300 MPa con¬ning pressure and shear strain rates of 1—10’3 s’1

to various amounts of ¬nite shear strain: (A“D) Cross-polarised transmitted light images

of the calcite microstructures; the amount of ¬nal strain is indicated for each sample and

represented by a non-deformed circle and the corresponding simple shear ellipse. (E“H)

Contoured c-axis pole ¬gures of calcite orientations measured by electron backscatter

diffraction analyses; the strength of the crystallographic preferred orientations is indicated

by the texture index γ. Data in (E) are taken from Pieri et al. (2001), and the data in (B“D)

and (F“H) from Barnhoorn et al. (2004). [Courtesy of A. Barnhoorn.]

184 Geological Deformation

which very large shear strains can be produced on the outer surface of a cylinder

subjected to continuous torques. Since the centre of the cylinder is ¬xed there is a

strong gradient in shear strain between the centre and the outside, but the outermost

5 per cent will have a similar shear con¬guration.

In such a way shear strains of more than 10 can be produced under elevated

temperature conditions, with consequent signi¬cant changes in microstructure. To

achieve the large shear strains in a reasonable time frame the strain rate needs to

be reasonably high (10’3 s’1 ), which is many orders of magnitude faster than in

actual geological processes. Nevertheless it is possible to study the evolution of

microstructure under severe strain. An example based on the work of Barnhoorn

et al. (2004) for large torques applied to Carrara marble is shown in Figure 10.2.

There is a progressive elongation of the grains as the torque progresses and many

approach the simple form of the elongated strain ellipsoid. The initial broad range

of grain orientations relative to the c-axis of calcite is progressively replaced by

coordinated grain orientations determined by the shear process. The sense of shear

and the amount of strain in natural rocks can often be deduced from the grain

elongation.

At higher temperatures and lower strain rates, recrystallisation causes

progressive destruction of such a fabric leading to the formation of a ¬ne-grained

ultra-mylonite, for which the deformation history can generally not be recognised.

Other features such as dyke displacement at the transition from the shear zone to

the surrounding rocks may then be used to deduce the shear sense and strain.

More complex textures arise when multiphase materials are subjected to severe

shear and this can lead to layering within mylonites.

10.1.3 Formation of crystallographically preferred orientations

As rocks are deformed they commonly develop a preferred orientation of

the crystallographic directions of the constituent minerals (cf. Figure 10.2).

The development of the preferred directions can be associated with two main

mechanisms: the rotation of grains or the process of recrystallisation. A rock has

to deform coherently with the grain boundaries remaining in contact, or cracks

will develop with ultimately brittle failure. One way in which a rock can submit

to a general shear strain with continuity across the grain boundaries is to have

homogeneous strain across the constituent grains. In order for this to be achieved by

slip within the grains, ¬ve independent slip systems have to act within each grain.

The condition for preservation of volume means that the trace of the incremental

strain must vanish, i.e., i eii = 0. In consequence only ¬ve of the elements of the

strain tensor eij can be independent. To produce this condition ¬ve independent

simple shears are needed, and hence the requirement of ¬ve independent slip

systems.

When solid-state diffusion is suf¬cient to allow climb, fewer slip systems are

needed. At such high temperatures some grain boundary sliding could also

10.1 Microfabrics 185

(a) (b) (c)

(d) (e) (f)

(g) (h) (i)

Figure 10.3. This sequence of frames simulates the process of dynamic recrystallisation in

quartz (Jessell & Bons, 2002). (a) Starting grain boundary geometry. The shading re¬‚ects

the internal energy of grains and subgrains, darker shading indicates higher energy. (b)

Situation after the imposition of a dextral shear of 0.65, modelled using 13 steps. (c“f) Post

deformation recovery and grain growth in intervals of 26 steps. (g) Representation of the

microstructure in the original state (a) with shading indicating lattice orientation. (h),(i)

Comparable lattice orientations for states (b),(f) showing the retention of the orientations.

[Courtesy of M. Jessell.]

contribute to crystal orientation. Other processes such as dissolution“precipitation

creep with solution of, e.g., calcite or quartz under pressure, can also move material

within a polycrystalline framework in the presence of some hydrous ¬‚uids to

develop fabric. Material is moved from the crystal faces under most compression to

those that are less compressed so that grains change shape without changing mass.

Numerical simulations of microstructures have proved a useful tool to guide

understanding of the way in which natural materials acquire their characteristics

(see, e.g., Jessell & Bons, 2002). Because of the different physical character of

the various processes, the topology and geometry of the polycrystalline aggregate

has to be carefully described and transferred to modules that represent the effect of

speci¬c grain-scale deformation processes. The simulation shown in Figure 10.3

186 Geological Deformation

for evolution of microstructure in quartz combines a ¬nite-element description of

deformation linked to a front-tracking model for grain boundary migration driven

by surface and defect energy terms.

Crystallographically preferred orientation of crystals is generally indicative of

prior deformation history, since it can often survive post-deformation annealing

processes. The presence of such preferential orientations has a signi¬cant impact

on the physical properties of the rock mass, especially with regard to the anisotropy,

e.g., for seismic wave propagation. The orientation of olivine is commonly invoked

as a major cause for the development of seismic anisotropy in the lithosphere that is

manifested through a time separation between S wave pulses of different orientation

or azimuthal variations in the phase velocities of surface waves.

The processes of diffusion and grain-boundary sliding should not lead to the

development of preferred orientations, but there can be dissolution and precipitation

processes that are related to a preferred crystallographic orientation

10.2 Macroscopic structures

There are three ways in which rocks can undergo large deformations:

(1) by ¬‚ow with normally an inhomogeneous strain ¬eld,

(2) by buckling or bending, with de¬‚ection of the rock layers, so that

considerable overall shortening can be achieved with only moderate

internal deformation

(3) by slip of one part of the body past another along fault surfaces or in narrow

zones, with little deformation away from the slip surfaces.

The deformation observed in rock outcrops is generally complex, re¬‚ecting

multiple phases of inhomogeneous strain that have interacted with the relative

strengths of different components of the rock. The outcome can be complicated

folding such as is illustrated in Figure 10.4 with an outcrop from the Proterozoic

Figure 10.4. Example of complex folding structures in the Proterozoic fold belt of the King

Leopold Range, Western Australia

10.2 Macroscopic structures 187

rocks of the King Leopold fold belt in northwestern Australia. Such convoluted

structures with tight folds and rapid changes of orientation imply substantial ¬nite

compression.

10.2.1 Multiple phases of deformation

Although we have developed the theoretical formulation of inhomogeneous strain

in Chapter 2, most of our prior examples have been based on the simpler condition

of homogeneous strain through a body. Once the strain ¬eld becomes both

¬nite and inhomogeneous a much wider range of behaviour is possible. The

mathematical development can be described through the use of a spatially varying

deformation gradient F and, as in (2.1.5), the in¬‚uence of multiple phases of

deformation is described by the ordered product of the deformation gradients F1F2.

A useful way to follow the nature of such complex deformation is to look at

the development of the local Eulerian ellipsoid generated from a sphere in the

original reference state (Section 2.2.2). This process is illustrated for multiple

two-dimensional deformation in Figure 10.5. A set of small circles is drawn on the

surface of the undeformed material, and the action of progressive deformation can

A B

C D

Figure 10.5. Illustrations of multistage inhomogeneous deformation leading to

concentrated belts of strong deformation indicated by the distortion of the initially circular

markers. (A),(B) Superimposition of oblique tighter folding (at a 60—¦ ) angle on a simple

fold. (C),(D) Simple shear superimposed on an inhomogeneous shear pro¬le.

188 Geological Deformation

be seen through the spatial behaviour of the Eulerian ellipses in the inhomogeneous

strain ¬eld.

Two simple cases are shown that lead to complex patterns of behaviour with

only two phases of deformation. In the upper row of Figure 10.5 A-B we consider

the application of a sinusoidal variation of strain, followed by a further sinusoidal

variation with shorter wavelength at an angle of 60—¦ to the original. In A the strain

is concentrated on the ¬‚anks of the fold, although the centre is strongly displaced,

there is only modest distortion of the circles near the axis. However, with the

superimposition of the second ¬nite strain in B, the pattern is much more complex.

Narrow bands of severe deformation develop and also the deformation is suf¬cient

that the circles have been distorted to non-elliptical shapes.

In the lower row, Figure 10.5 C“D, we consider the imposition of shear strain. At

¬rst a concentrated band of shear deformation is applied with a hyperbolic tangent

(tanh) pro¬le in C, so that the strain fades away towards the edges. This is followed

by a simple shear at right angles in D. The combination signi¬cantly modi¬es the

direction of the Eulerian triad in D compared to C, and some care is needed to infer

the original deformation ¬eld after removal of the homogeneous simple shear. The

situation becomes rapidly more dif¬cult when the second phase of deformation is

itself inhomogeneous.

The presence of a strain ¬eld will lead to the distortion of the original shapes

of objects in the rock. Over regions with reasonably homogeneous strain ¬elds,

the reactions of natural markers can provide useful measures of strain. Initially