. 7
( 16)



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)
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
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)
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
E= ,
d q hω /k T (9.2.9)
(2π)3 e q B ’1

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
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)
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
(2π)3 6π2
3 πqD = , qD = .
i.e., (9.2.12)
The associated cut-off frequency
1 2
ωD = qDvm, with vm = 3 +3 , (9.2.13)
v3 vS
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)

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
Cv = 3nNkB 3 ,
dz (9.2.15)
˜D ez ’ 1

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
˜D /T
z4 ˜D
dz 2 = , (9.2.16)
z T

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
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
γ γ T
Pth (V, T ) = E D= 9nNkBT ,
dz z (9.2.19)
V V ˜D e ’1
where γ is the Gr¨ neisen constant, which under the Debye assumptions is
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
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:
= μ, (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)
σ d0 TM
™ =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
σ 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,
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
n n
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
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

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

Temperature [oC]

Pressure [GPa]

Lower Mantle Outer Core Inner Core

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)
ha = hb + 1 (pa ’ pb)(1/ρa + 1/ρb). (9.4.4)
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
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
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]
100 200 300 400 500

Temperature [˚C]

L St

High-P Clino-enstatite
Continental Shield Fe
1000 Geotherm [unquenchable]

Spinel +
Low-P Clino-enstatite
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


Elastic constant [GPa]





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

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
H = Ven + ’ (9.5.2)
|ri ’ rj|
i<j i=1

The external potential Ven arises from the interaction of electrons with the M
atomic nuclei,
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
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
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
resemble those for a single electron. The Kohn“Sham equations for a set of orbitals
{φi} are
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
ρ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
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
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
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,
φ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
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
D(q) = exp(iq · rij). (9.5.11)

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
1 ”V
F= + kBT ln 2 sinh , (9.5.12)
2K V 2kBT

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
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
1 ”V ¯
γhωq 1 coth
= = γE(T ). (9.5.14)
K V 2kBT

Here E(T ) is the energy in the lattice modes at temperature T and the Gr¨ neisen
constant γ relates the dilatation to the mean thermal energy density
”V ¯
= KγE(T ). (9.5.15)
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
(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
mi 2 = fi(r1, r2, . . . , rn), i = 1, . . . , n , (9.5.16)
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

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

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

γ =0
J ˜ 1.0

γ ˜ 1.4

γ ˜ 0.5
J = 1.8

γ ˜ 3.4

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

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



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


. 7
( 16)