49

GEOL0026: Session C; Equations of State

I am indebted to my colleague Prof. Lidunka Vočadlo, from whom I have “borrowed”

much of this material which used to form part of the module “C365: Physics of

Planetary Interiors”. I have, however, made some minor corrections and amendments to

the text, added some additional material in italics and also added 3 short sections at the

end of the notes.

Equations of State

An equation of state, EOS, describes how the volume, V, or density , of a material

system varies as a function of pressure, P, and temperature, T; it therefore allows us to

determine a mineral property ( or V) at depth (P and T).

The simplest example of an EOS, the EOS for an ideal gas, is given by:

Eq. 75

where R is a constant, the gas constant; R = 8.31451 J K-1 mol-1. In other words, for a

molar volume of ideal gas, V, experiencing an external pressure, P, there will be an

associated temperature of the gas given by PV / R.

More complex EOS are vital to planetary scientists, because in order to interpret seismic

data to give Earth structure, or to rationalise planetary densities from moments of inertia

and mass data, it is necessary to either:

(a) "compress and heat" mineral data to the PT state of a planetary interior, or

(b) "decompress and cool" planetary data to ambient conditions.

The simplest isothermal EOS for a solid is just the definition for incompressibility (also

called the bulk modulus), K:

Eq. 76

The simplest isobaric EOS for a solid is just the definition for the thermal expansion

coefficient, :

Eq. 77

For most materials, the effect of pressure within planetary interiors is far greater than

that of temperature; it is therefore easier to consider initially an isothermal EOS and

then add on a thermal expansion correction. To this end, we shall first concentrate only

on the effect of pressure, i.e., on isothermal equations of state.

50

Isothermal Equations of State

We shall consider two types of isothermal EOS: infinitesimal EOS and finite

strain EOS. The first treats volume expansion in terms of infinitesimally

increasing increments in V, i.e., integrating over a volume or pressure range; the

second treats volume expansion in term of finite differences in strain, i.e.,

considering the subsequent volume change as the original volume plus a little

bit. Although the mathematics gets a bit tricky, especially for finite strain theory,

the essential points of both methods are outlined below.

Infinitesimal Equations of State:

Isothermal EOS begin in their simplest form with the incompressibility, K. If we

assume K is a constant (which is only true for small P), then integration of Eq.

76 gives:

Eq. 78

therefore,

Eq. 79

so

Eq. 80

and similarly

Eq. 81

where K0 = K (0,T) is the incompressibility at P = 0, T = constant.

However, if K was a constant, remaining unchanged with increasing pressure,

then as P → ,  → , which we know is impossible; indeed, both seismology

and atomistic analysis show us that K increases with P, i.e., dK/dP > 0. In other

words, the more you squash something, the harder it is to squash it further. At an

atomic level, the more squashed a material becomes, the larger the repulsive

forces, holding the atoms apart, become to resist the pressure.

If we decide that it is not valid to assume that K in Eq. 76 is a constant, we can

repeat the whole integration process again, but at the next level of

approximation; writing dK/dP = K0′, we can put

Eq. 82

i.e., K increases linearly with pressure. After integration, this gives:

Eq. 83

or,

51

Eq. 84

This is the Murnaghan Integrated Linear Equation of State, MILEOS.

However, this equation is also still approximate since K' itself is also a function

of P, i.e., K'  constant, especially at very high pressures, or for “squashy”

materials with a small K. In principle, we could just repeat the above integration

procedure with:

Eq. 85

But this is by no means a simple calculation, and, for minerals, K" is very

difficult to measure experimentally, and may, or may not, itself vary with P.

Consequently, this led to the development of finite strain theory, which more

readily takes into account the variation of K with pressure.

Finite Strain Equations of State:

Finite strain theory is quite complex in detail and is based on the analysis of a

deformed body undergoing strain.

We can illustrate the essence of the difference between “infinitesimal strain

methods” and “finite strain methods” by considering a simple 2-dimensional

example.

Suppose we take a square of material of unit side, and then deform it so that it

changes into a rectangle (1 + ε1) by (1 + ε2), as shown below.

Following the method used in Session A, we can write the strain tensor as

      2 1 0 0  

The area of the square before it was

deformed was 1 unit squared.

After deformation the area is

(1 + ε1) x (1 + ε2) = 1 + ε1 + ε2 + ε1ε2

In infinitesimal stain theory, we ignore all terms in ε1ε2, i.e. we ignore the

shaded part of the diagram. The relative change in area, ΔA/A, is then equal to

ε1 + ε2 = the sum of the diagonal terms (this is called “the trace”) of the strain

matrix. Clearly as the strains become larger this approximation becomes worse.

In 3-dimensions a similar argument applies to the ways in which we model the

volume change as pressure is applied to the material.

52

Consider a solid body undergoing such strain as a result of some applied

pressure; if the initial co-ordinates (or state) defining its volume are xi, and after

deformation the final co-ordinates are Xi, then there is some displacement vector

associated with this strain, ui. The deformation itself may be defined in two ways

as follows:

1) The Lagrangian scheme, where the final co-ordinates are written in

terms of the initial co-ordinates and the displacement vector thus:

Eq. 86

and therefore

Eq. 87

where  = l / l, i.e., the fractional change in length, l, in each of three

co-ordinate components.

2) the Eulerian scheme, where the initial co-ordinates are written in

terms of the final co-ordinates and the displacement vector:

Eq. 88

and therefore

Eq. 89

where, for compression,  = - l / l (or  = −V/V, with  = 3)

When considering infinitesimal strain, the above two schemes are equivalent,

i.e., for infinitesimally small volume changes, the associated strain may be

defined from either the starting or finishing co-ordinates; but for finite strain, the

above two schemes give differing results, and we therefore have to state which

scheme is being used in any subsequent analysis. In this instance we shall be

considering the Eulerian scheme which is more meaningful for finite strain

conditions, since every quantity is expressed in terms of the coordinates of the

points in the strained solid, which are the experimentally accessible ones.

Using the Eulerian scheme, and a lot of mathematical manipulation involving an

expansion of the Helmholtz free energy, F, about a factor f, called the

compression of the solid, to second order, i.e., F  af2, where f = -(/3 – ½ 2/9)

 -/3  -3/3 (where “≈” applies to infinitesimal strains), and after appropriate

substitutions for K and P we get the following:

Eq. 90

This is the 2nd order Birch-Murnaghan equation of state, BMEOS.

We can now substitute this into the definition of K (Eq. 91)

53

Eq. 91

which gives

Eq. 92

and

Eq. 93

so for f → 0, i.e., very small strains, K'0 = 4, as observed in many experiments.

Intuitively, we would expect the result from finite strain theory for finite strains

approaching zero to compare well with infinitesimal strain theory; indeed when

we put K' = 4 in the MILEOS and set f = 0 in the 2nd order BMEOS, we see that

the MILEOS overestimates the predicted pressure by about 3% for a  / 0 ~

1.25 (the ratio inferred at the core-mantle boundary).

This approach can be extended to include an expansion of the Helmholtz free

energy about the strain factor to third order (i.e., F  af2 + bf3) which, after

tedious calculation, leads to:

Eq. 94

This is the 3rd order Birch-Murnaghan equation of State.

If K0' = 4, then this equation reduces to the 2 nd order BMEOS, but the 3rd order

BMEOS must be used when dealing with very high pressures where dK/dP

varies significantly with pressure.

Finally, there is also a 4th order BMEOS (F  af2 + bf3 + cf4), which allows for

K" to vary as a function of pressure, and is therefore appropriate when dealing

with extremely condensed materials (possibly relevant in the deep interiors of

the gas giant planets), although such a variation in K" is rarely measured

experimentally. This very complex EOS takes the form:

Eq. 95

where

54

)24/143()7()8/3()8/3( 0000 +−+= KKKK

The 4th order BMEOS reduces to the 3rd order BMEOS if  is set to zero;

however, for this to be the case, K"  0 since rearranging Eq. 95 gives:

Eq. 96

Therefore, the use of the 3rd order BMEOS implicitly assumes a non-zero value

for K", even if it is not measured.

Logarithmic Equations of State

It has been known for some time that the convergence of the polynomial

expansion for Eulerian strains is poor at very large strain values when K'  0. In

other words, the Birch-Murnaghan formalism breaks down at high pressures.

Poirier and Tarantola (1998) used the Hencky logarithmic strain to derive a new

equation of state. The Hencky strain is the same as the Eulerian strain for small

strains and exhibits better convergence behaviour at large strains. Thus an

equation of state derived in this way should offer better performance at high

pressure than the BMEOS.

The Hencky strain, H, is defined as,

Eq. 97

where dl / l is the fractional change in length under uniaxial compression.

Hence, for deformation under hydrostatic pressure,

Eq. 98

As with the derivation of the Birch-Murnaghan EOS, the pressure is written in

terms of the free energy derivative with respect to the strain,

Eq. 99

and the free energy is expanded as a series in powers of n,

Eq. 100

leading to a logarithmic equation of state:

Eq. 101

When this expansion is truncated at N = 2, we find;

55

Eq. 102

and so

Eq. 103

This is the 2nd order Logarithmic Equation of State

For N = 3,

Eq. 104

This is the 3rd order Logarithmic Equation of State

The 3rd order LNEOS is equivalent to the 3rd order BMEOS for small strains, but

diverges from it at large strains. BMEOS3 can give negative pressures for large

strains with K' < 4, whereas the logarithmic formulation avoids this.

As with the BMEOS, the expansion can be continued, yielding for N = 4, the 4th

order LNEOS:

Eq. 105

where

See Vocadlo et al. (2000) for a comparison of the 4th order LNEOS with BM

Equations of State.

However, it has recently been recognised (e.g., Stacey, 2000) that the finite

strain approach to equations of state is fundamentally flawed. Holzapfel

observed that K' should tend to K' = 5 / 3 (the value for a free-electron Fermi

Gas) under high compression. Few equations of state satisfy this condition,

56

although Shanker (1999) offers a 'universal' equation of state derived with this

condition in mind.

In Summary:

The simplest isothermal EOS is the definition for incompressibility (bulk

modulus). Using infinitesimal strain theory, integration of K yields an EOS

which implies that material is infinitely compressible, which we know is not the

case. Therefore, K must vary with P (as the more compressed matter is, the

harder it is to squash), and under this next level of approximation we arrive at

the MILEOS which therefore has a non-zero value for K', but K" = 0. An

advantage of the MILEOS is that it is mathematically more tractable than the

BMEOS, so it is easy to obtain either P from a given V, or V from a given P.

However, K' itself is known to vary with P, especially at high pressures, and

therefore we turn to finite strain theory to account for this pressure dependence

of K'. Under the Eulerian scheme, where the initial volume is given in terms of a

fraction of the strained volume, and expanding the Helmholtz free energy to

second order in the strain factor, we arrive at the 2nd order BMEOS. This gives

P, K and K0' as functions of  / 0. If the strain factor is vanishingly small, then

K' --> 4, i.e., approaching the no strain P = 0 case. For  / 0 = 1.25, the pre- dicted pressure differs from that derived from the MILEOS by only 3%. K" = 0.

When dealing with very high pressures, expanding the Helmholtz free energy to

third order in the strain factor results in the 3rd order BMEOS. This is equivalent

to the 2nd order BMEOS when K'= 4, and to the 2nd order BMEOS plus terms in

K" when K  4; the 3rd order BMEOS implicitly assumes a non-zero value of K".

At extremely high pressures, when K" is known for the system under

investigation, expanding the Helmholtz free energy to fourth order in the strain

results in the 4th order BMEOS. When  = 0, 4th order BMEOS = 3rd order

BMEOS; but setting  = 0 implies a non-zero value for K".

Problems with the formulation of the BM equations of state can be dealt with by

employing the Hencky strain to yield high-order logarithmic equations of state,

which are better suited to the compressive regime at the base of the mantle and

in the Earth's core. Be aware, though, when choosing an equation of state, that

most finite-strain equations of state are thermodynamically flawed at very high

compressions (of order teraPascals). When considering the behaviour of material

inside gas giants (Jupiter's core pressure is ~5 TPa), brown dwarfs and stars,

then relativistic nuclear equations of state (e.g., the Thomas-Fermi EOS) are

more appropriate.

The summary in summary: If K' = 4, then use the MILEOS or the 2nd order

BMEOS; if K'  4, use MILEOS or 3rd order BMEOS; if K"  0, use 4th order

BMEOS. However, the higher order the EOS that is used, the more complicated

the calculation. The 4th order LNEOS is favoured for studies of the deep mantle

and core.

57

Thermal Equations of State:

Thermal EOS describe how a solid reacts to changes in temperature. When a

solid is heated, the associated change in thermal pressure of the system results in

thermal expansion. The thermal pressure is the increase in internal pressure

caused by heating at constant volume, e.g., the pressure you would feel if

someone was heating up a balloon, and you were trying to hold the balloon at

constant volume; this given by:

Eq. 106

where γ is the thermodynamic Grüneisen parameter. Integrating at constant V

and  gives us:

Eq. 107

Therefore:

Eq. 108

where U is the change in internal energy (1st law thermodynamics).

This is the Mie-Grüneisen equation of state.

In real systems such changes in thermal pressure would cause thermal

expansion.

The change in internal energy is difficult to measure, but by substituting Eq. 106

into Eq. 107 we get:

Eq. 109

At high T, KT ~ constant, but  =  (P), i.e., it varies with pressure; this

variation of  with P is directly related to the variation of K with T, and is given

by:

Eq. 110

These effects are more readily expressed via the Anderson-Grüneisen

parameter, which takes both an isothermal and adiabatic form.

a) The isothermal Anderson-Grüneisen parameter, T, is defined by:

Eq. 111

and since:

58

Eq. 112

therefore:

Eq. 113

so:

Eq. 114

i.e., this shows us how thermal expansion varies with relative volume change

and, therefore, with pressure.

At low pressures T is approximately constant, although structure dependent, but

at very high pressures, T varies with P also, via:

Eq. 115

where  is an empirically derived constant,  = V / V0, and the subscript T0

represents the value of the isothermal  measured at P = 0.

Measuring the Anderson-Grüneisen parameter is quite difficult; if it is to be done by

diffraction methods it requires a lot of accurate data to be collected as a function of both

P and T; a recent example of such work at UCL, in which we attempted to measure this

quantity for (Mg,Fe)O, is given in the Electronic Supplement (file MgFeO.pdf).

b) When dealing with adiabatic systems, e.g., seismic structure and EOS, the

Anderson-Grüneisen parameter is defined thus:

Eq. 116

and is related to the isothermal parameter via:

Eq. 117

where  is the Grüneisen parameter.

Equations of State at P and T

An all-encompassing EOS at pressure and temperature, giving V(P,T), may be

obtained by first using an isothermal EOS to calculate V(P,0); the Anderson- Grüneisen parameter may then be used to calculate (P) from (0). Thus, finally

we may calculate V(P,T) from V(P,0) and (P).

59

Shock Equations of State

An alternative way to the methods described above of obtaining EOS is to use

shock experiment data. A shock experiment can be performed by firing a gun at

a target in a controlled manner and analysing the resulting particle and shock

wave velocities. In this way, instantaneous short-lived high pressures may be

simulated (P > 350 GPa), and the effects on the material observed.

A projectile is fired at a target and sends shock waves through it with velocity us.

The large macroscopic particles within the target move with velocity u (u < us).

The compressive wave travels through the target, hits the rear and is reflected as

a dilation, which then blows the sample to bits.

Using suitable methods (electric counting, optical path changes, etc.) u and us

can be measured. Using the laws of conservation of mass, momentum and

energy, gives:

Eq. 118

Eq. 119

and

Eq. 120

therefore the change in internal energy is given by:

Eq. 121

The shock event is neither adiabatic nor isothermal; it is not adiabatic because it

is non-reversible, since the sample gains kinetic energy as well as potential

energy, and it is not isothermal because the sample heats up.

The locus of the P vs.  points produced by a series of shock experiments is

called the Hugoniot.

The 'mixed phase' region is interpreted as being due to 'phase changes'. The

pressure and resulting microstructures bear no relation to equilibrium phase

changes, as the kinetics are too sluggish, but are often due to the production of a

high density glass phase.

Once the Hugoniot data has been obtained, it must be reduced to either an

isothermal or adiabatic EOS for it to be useful. The increase in internal energy

caused by shocking the sample gives rise to a temperature increase and therefore

a thermal pressure. To convert Hugoniot data to an adiabatic EOS, the thermal

pressure must be removed. This is done via the Mie-Grüneisen equation (seen

earlier):

Eq. 108

60

This is often large, so an accurate value of  is essential for this EOS to be

successful.

To obtain an isothermal EOS, the adiabatic EOS is converted via

thermodynamic identities.

References

Jackson. I., and S.M. Rigden (1996) Analysis of P-V-T data: Constraints on the

thermoelastic properties of high-pressure minerals. Phys. Earth. Planet. Int. 96,

85-112.

Poirier. J.-P., and A. Tarantola (1998) A logarithmic equation of state. Phys.

Earth. Planet. Int. 109, 1-8.

Shanker. et al. (1999) On the universality of phenomenological isothermal

equations of state for solids. Physica B 271, 158-164.

Stacey. F.D. (1995) Theory of thermal and elastic properties of the lower mantle

and core. Phys. Earth. Planet. Int. 89, 219-245.

Stacey. F.D. (2000) The K-primed approach to high-pressure equations of state.

Geophys. J. Int. 143, 621-628.

Stacey. F.D., and D.G. Isaak (2001) Compositional constraints on the equation

of state and thermal properties of the lower mantle. Geophys. J. Int. 146, 143- 154.

Vocadlo. L., J.-P. Poirier, and G.D. Price (2000) Grüneisen parameters and

isothermal equations of state. Am. Min. 85.

Williams. Q., and E. Knittle (1997) Constraints on core chemistry from the

pressure dependence of the bulk modulus. Phys. Earth. Planet. Int. 100, 49-59.

61

Some further notes (by IGW)

1/. Equations of State at P & T

The derivatives of the isothermal incompressibility, KT, and the isobaric volume

coefficient of thermal expansion, αP, are linked via the relationship

P T TT P T K KP         =        2 1

(e.g. Bina and Helffrich, Annual Rev. Earth Planet. Sci., 20, 527-552, 1992) and thus

only one of these two derivatives need be measured. As mentioned above (pages 57 &

58), the temperature dependence of KT is commonly expressed in terms of the

isothermal Anderson-Grüneisen parameter, δT ; this dimensionless quantity may be

written in a number of ways, possibly the most useful being

δT = -(∂lnKT/∂lnV) = -(1/αPKT)(∂KT /∂T)

or

δT = (∂lnαP/∂lnV) = -(KT/αP)(∂α/∂P)

(see e.g. Bina and Helffrich, 1992; Poirier, Introduction to the Physics of the Earth’s

Interior). Knowledge of the first derivatives of K and α through δ is, therefore, an

essential minimum requirement if we wish to extrapolate accurately data obtained at

modest P/T to the conditions existing in the Earth’s lower mantle.

We often have experimental data that are scattered in P and T, as shown in the diagram.

The procedure that we can use to fit such PVT values is outlined below.

We choose a suitable reference

temperature, To.

This then allows us to find Vo(T) by

thermally expanding up the vertical

axis at zero pressure, using, e.g. from

Vo(T) = Vo(To)[1 + αo(T-To)]

We can then calculate the equations of the different P-V isobars, for example using the

BM-3. To do this, however, we need to know how Ko and Ko' vary with temperature, i.e.

we need values of (dKo/dT) (and, ideally, of dKo'/dT). If we have enough P-V data

points, we can use non-linear, least-squares fitting to adjust all six parameters

Vo(To), αo, Ko(To), Ko'(To), (dKo/dT) and (dKo'/dT)

to get the best fit to the data.

This provides us with a method by which to determine δT, the isothermal Anderson- Grüneisen parameter. Note that unless the data are very precise, covering a wide PT

range, it may well be necessary, for example, to assume (dKo'/dT) = 0.

62

2/. Simple Isobaric Equations of State

At e.g. P = 0, the simplest equation of state is obtained by integrating the expression for

the thermal expansion coefficient. From

dT dV V 1 =

we obtain

 = dTV dV 

i.e.

CdTV += )ln(

If we now choose a reference temperature, To, at which V = Vo, we obtain

CdTV T To += )ln(

and hence C = ln(Vo)

Thus

=      T To o dT V V ln

(2.1)

We cannot go any further unless we know how α varies with T. If α is temperature

independent, i.e. α = αo, then

)(ln oo o TT V V −=      

or

)( oo TTo eVV −= 

If the temperature difference is not too large, then )(1 )( oo TT TTe oo −+=−  , and so

 )(1 ooo TTVV −+= 

an expression with which you are doubtless familiar.

If α cannot be assumed to be independent of T, then an expression such as

α = ao + a1(T-To) + ……

must be used instead on the right-hand side of equation (2.1).

Often, however, for practical purposes, one can dispense with this treatment and simply

fit the V(T) data directly to a polynomial such as

V = Vo + b1(T) + b2(T

2) + b3(T

3) + ……

However, if you do this you, should be very cautious if you extrapolate the data

beyond the range of measurement – polynomials can go “haywire” very easily!

63

3/. More Physically Meaningful Isobaric Equations of State

More physically meaningful interpretations of V(T) data may be obtained using Grüneisen

approximations for the zero pressure equation of state (see Wallace, Thermodynamics of

Crystals, pub Dover, 1998; a truly awesome book!) in which the effects of thermal

expansion are considered to be equivalent to elastic strain induced by thermal pressure.

It has been found empirically that the thermodynamic Grüneisen parameter, γ, given by

γ = αKV/CV

(where α is the thermal expansion coefficient; K, the incompressibility (bulk modulus);

V, the volume and CV, the specific heat at constant volume) is reasonably constant.

Rearranging this equation and integrating with respect to T, taking γ as constant, we obtain

 = VdTKdTCV 

but

dT dV V 1 =

and so

  =      = KdVVdT dT dV V KdTCV 1 

But, the left-hand side of this equation  = )T(EdTCV 

and so, if we assume that K is independent of temperature (which is a bit iffy!), we then

find that

0 0 V K )T(E )T(V += 

(3.1)

Wallace, gives a better expression for V(T) “to second order”:

0 0 )( )( )( V TbEQ TEV TV + − =

(3.2)

where Q=V0K0/ and b=(K0'-1)/2;  is a Grüneisen parameter (assumed constant), K0 and

K0' are the incompressibility and its first derivative with respect to pressure respectively

at T = 0 and V0 is the volume at T = 0. The internal energy, E(T), may be calculated using

the Debye approximation (Session B notes, or Cochran’s book) from:

 −       = T x D B D e dxxT TnkTE   0 3 3 1 9)(

(3.3)

where n is the number of atoms in the unit cell, kB is Boltzmann's constant, and D is the

Debye temperature.

Thus, if we fit equations (3.1) or (3.2) to V(T) data obtained from diffraction experiments,

we are able not only to model the expansion curve over a very wide temperature range,

but also to determine useful material parameters such as D and K0'.

64

An example of this work is shown below, taken from a study of ε-FeSi in which we first

used this approach (Vočadlo, Knight, Price and Wood; Phys. Chem. Minerals, 29, 132- 139, 2002). Note that if you read this paper, you will find that some of the symbols used

in the equations are different from those given above (I have changed them here so that

they are consistent with those used in Session B).

A further example of this type of interpretation of thermal expansion data, in which the

analysis is extended to the axial expansions, as well as the volumetric expansion, can be

found in the paper on post-perovskite structured CaIrO3 by Lindsay-Scott et al. (file

CaIrO3.pdf in the Electronic. Supplement).

65

Exercises:

1/. Show that

  d dP dV dP V =−

2/. Assuming that the incompressibility is given by K = Ko + Ko 'P, derive the

Murnaghan Integrated Linear Equation of State, expressing the equation in terms of P

and V rather than P and ρ.

3/. If Ko = 200 GPa and Ko ' = 4, what is the difference in pressure predicted by the

Murnaghan Integrated Linear EoS and the Birch-Murnaghan 3rd order EoS for

compactions V/Vo = 0.9, 0.8, and 0.7?

4/. Calculate the thermal expansion coefficient of a material that has been compressed to

80% of its original volume, assuming that the ambient pressure value of α is 3 x 10-5 K-1

and that the Anderson-Grüneisen parameter, δT =1.47.

5/. (i) The thermal expansion coefficient, α, is defined as

dT dV V 1 =

The table on the right shows experimental data

for the change in V with T for ε-FeSi, between

275 K and 607 K.

Using Excel (or equivalent) plot a graph of

V vs T, fit a straight line to the data, and hence

find the value of α at 300 K (assuming that α

is independent of temperature).

(ii) The incompressibility (bulk modulus), K, of a material is defined by

dV dP VK −=

The table on the right shows experimental data

for the change in V with P for ε-FeSi, between

0 GPa and 3.006 GPa.

Plot a graph of V vs P, fit a straight line to the data,

and hence find the value of K at 0 GPa (assuming

that K is independent of Pressure).

66

(iii) The thermodynamic Grüneisen parameter, γ, is given by

γ = αKV/CV

(see Lecture notes for definitions of the symbols used).

The heat capacity of ε-FeSi at 300 K and 0 GPa is 48.11 J mol-1 K-1, and there are 4

formula units in each unit cell.

Calculate the molar volume of ε-FeSi (Avogadro’s number = 6.022 x 1023 per mole) and

hence find the value of the thermodynamic Grüneisen parameter at 300 K and 0 GPa.

(iv) Explain briefly the additional factors that would have to be taken into consideration

when applying the analysis used in parts (i) and (ii) to materials in planetary interiors.

How might these affect the value of the thermodynamic Grüneisen parameter?

6/. Mineralogical Models

In this final exercise, taken from Dr Vočadlo’s course (C365, Physics of the Earth &

Planetary Interiors, we shall investigate whether the PREM seismological model

Preliminary Reference Earth Model) is consistent with an upper mantle mineralogical

model of predominantly olivine and garnet, assuming a given geotherm.

Method:

We need to find to show how a) density and b) seismic velocities of olivine and pyrope- rich garnet vary as a function of pressure and temperature, and therefore, for a given

geotherm, how they will vary as a function of depth. We can then compare our

calculations with PREM to test the validity of this mineralogical model.

Mineralogical Data (Remember, if necessary, to change the units!)

Olivine Garnet

ρ0,300 /g cm-3 3.353 3.705

K0,300 /GPa 128 169.4

K' 4 4

α0,300 /x10-5 K-1 2.66 2.36

δt 6.59 6.27

dK/dT /GPa K-1 0.0223 0.02513

γ 1.25 1.5

dKs/dT /GPa K-1 0.018 0.0195

μ0,300 /GPa 78.1 92.6

dμ/dP 1.71 1.56

dμ/dT /GPa K-1 0.0136 0.0102

67

Geotherm Data (Stacey, 1977)

Depth /km P /GPa T /K

0 0 300

72 2.2 1035

122 3.5 1266

172 5.4 1480

222 7.1 1670

272 8.5 1823

322 10.5 1945

372 11 2025

422 14 2085

472 16 2130

522 18 2165

572 20 2197

622 21.6 2225

670 23.8 2250

771 28.3 2290

871 32.8 2335

971 37.3 2372

1071 41.9 2407

1171 46.5 2445

1271 51.2 2483

a) Calculating the density profile as a function of depth:

Using the MILEOS, a density-depth profile may be obtained in two ways: (i) by

calculating density as a function of pressure at constant temperature, then as a function

of temperature at this pressure; (ii) by calculating density as a function of temperature at

zero pressure, and then as a function of pressure at this temperature.

(i) Going up in pressure first, and then temperature.

a) Obtain the density as a function of pressure at 300K via the MILEOS:

        P K K +1

=

0 K 1 0,300P,300 

68

b) Calculate the Thermal expansion coefficient as a function of pressure via the

Anderson-Grüneisen parameter:

            P,300 0,300 0,300P,300 T

=

c) Finally, calculate the density as a function of pressure and temperature via the

thermal expansion coefficient:

T)-(1

=

P,300P,300TP, 

(ii) Going up in temperature first, and then in pressure.

a) Obtain the density as a function of temperature via the zero pressure thermal

expansion coefficient:

T)-(1

=

0,3000,300T0, 

b) Then find the incompressibility as a function of temperature from:

T dT dK -K

=

K 0,300T0, 

(we assume that K' does not vary)

c) Finally, calculate density as a function of temperature and pressure via the MILEOS.

          P K K +1

=

T0 K 1 T0,TP, , 

Calculate the density profile with depth using both methods and compare the results;

by how much do they differ?

In general, the first method is probably considered to be

the more reliable – why?

b) Calculating the seismic velocity profiles:

Seismic velocities may be expressed in terms of the adiabatic incompressibility, the

shear modulus, and the density. The adiabatic incompressibility, KS, may be obtained

from the isothermal incompressibility, KT, as a function of pressure and temperature;

this may then be put into the Vp and Vs equations, along with the density profile

calculated previously. The stages in this process are as follows:

a) Get KT at 300 K as a function of pressure, i.e. find KT (P,300):

There are a variety of ways in which you can do this.

The basic assumption of the

MILEOS is that:

PKK

=

K TT 0,300300P +,

69

As an alternative, you should also evaluate KT (P,300) from the 2 nd-order Birch- Murnaghan EOS, by putting your calculated density profile into the equation below.

])/(5)/(7)[2/( 3/50 3/7 0,  −K

=

K TT 0,300300P

b) KS may now be calculated as a function of pressure, using the Grüneisen parameter,

from

T)+(1K

=

K P0,300Ts P,300P,300 

We have not covered this conversion from isothermal (constant temperature) to

adiabatic (no heat flows in or out of the body) parameters in the lectures – it is one of

the many uses of the Grüneisen parameter; you will find it discussed in, for example,

Poirier: Introduction to the Physics of the Earth’s Interior.

c) KS may now be calculated also as a function of temperature, via dKS/dT.

T dT dK -K = K s ss P,TP  300,

d) The shear modulus as a function of pressure may be found from dμ/dP:

P dP d +

=

0,300P,300  

e) and the shear modulus as a function of temperature may be found from dμ/dT:

T dT d -

=

P,300TP    ,

f) Finally (!) we can put KS, μ and ρ, all as a function of pressure and temperature, into

the expressions for VP and VS:

  TP TPs p 3 4 +K

=

V TP , ,,

  TP TP s

=

V , ,

c) Comparison with PREM

Having obtained the density and seismic velocity profile within the given geotherm for

both olivine and garnet, these may now be plotted as a function of depth together with

the PREM model for comparison.

The PREM data are given in the table below. They may also be obtained electronically

as an Excel spreadsheet from the GEOL0026 Moodle site or the Electronic Supplement

70

In the light of your calculations, comment on the possible mineral phases in the

upper mantle. Discuss briefly how the iron content of the minerals might affect

your results, and also the shortcomings in this type of calculation using the

MILEOS.

References:

Poirier, J. 1991 Introduction to the Physics of the Earth's interior.

Weidner, D.J. 1986 Mantle Model Based on Measured Physical Properties of Minerals.

In Chemistry and Physics of Terrestrial Planets, ed. S.K. Saxena.

Duffy, T.S. and Anderson, D.L. 1989 Seismic Velocities in Mantle Minerals and the

Mineralogy of the Upper Mantle. J. Geophys. Res. 94 1895-912.

The PREM Model

PREM data (from Poirier), units are as follows:

Depth in km

Density in Mg per cubic

metre

Incompressibility, K, in kbar (10 kbar = 1GPa)

Shear modulus, μ, in GPa

Seismic velocities, Vp and Vs in km/s

Depth

Density K μ Vp Vs

24.4 3.38 1315 68.2 8.112255 4.491939

40 3.38 1311 68 8.100089 4.485348

60 3.38 1307 67.7 8.085466 4.475443

90 3.37 1303 67.4 8.082781 4.472136

115 3.37 1287 66.5 8.03122 4.442177

185 3.36 1278 66 8.014124 4.432026

220 3.36 1270 65.6 7.989328 4.418576

220 3.44 1529 74.1 8.553865 4.641196

265 3.42 1579 75.7 8.699555 4.704732

310 3.49 1630 77.3 8.731373 4.706272

355 3.52 1682 79 8.815233 4.737424

400 3.54 1735 80.6 8.908935 4.77162

400 3.72 1899 90.6 9.13901 4.935062

450 3.79 2037 97.7 9.387111 5.07724

500 3.85 2181 105.1 9.646119 5.224816

550 3.91 2332 112.8 9.904919 5.371136

600 3.98 2489 121 10.15252 5.513802

635 3.98 2523 122.4 10.21748 5.545608

670 3.99 2556 123.9 10.26955 5.572489

670 4.38 2999 154.8 10.75145 5.944953

721 4.41 3067 163.9 10.91332 6.096354

771 4.44 3133 173 11.06865 6.242112

871 4.5 3303 179.4 11.24969 6.314006

971 4.56 3471 185.6 11.41873 6.379793

1071 4.62 3638 191.8 11.58007 6.443225

1171 4.68 3803 197.9 11.73211 6.502794

1271 4.73 3966 203.9 11.88801 6.565655

1371 4.79 4128 209.8 12.0241 6.618125

1471 4.84 4288 215.7 12.1662 6.675786

1571 4.9 4448 221.5 12.29014 6.723398

71

1671 4.95 4607 227.3 12.42161 6.77637

1771 5 4766 233.1 12.5491 6.827884

1871 5.05 4925 238.8 12.67179 6.876564

1971 5.11 5085 244.5 12.77917 6.917178

2071 5.16 5246 250.2 12.89643 6.963359

2171 5.21 5409 255.9 13.01188 7.008358

2271 5.26 5575 261.7 13.12729 7.05357

2371 5.31 5744 267.5 13.24168 7.097651

2471 5.36 5917 273.4 13.35671 7.141951

2571 5.41 6095 279.4 13.47301 7.186453

2671 5.46 6279 285.5 13.59114 7.231139

2771 5.51 6440 290.7 13.68295 7.263513

2871 5.56 6537 293.3 13.70794 7.263043

2891 5.57 6556 293.8 13.71244 7.262703

2891 9.9 6441 0.00E+00 8.066016 0

2971 10.02 6743 0.00E+00 8.203378 0

3071 10.18 7116 0.00E+00 8.360728 0

3171 10.33 7484 0.00E+00 8.511708 0

3271 10.47 7846 0.00E+00 8.656669 0

3371 10.6 8202 0.00E+00 8.79644 0

3471 10.73 8550 0.00E+00 8.926541 0

3571 10.85 8889 0.00E+00 9.051313 0

3671 10.97 9220 0.00E+00 9.167737 0

3771 11.08 9542 0.00E+00 9.28004 0

3871 11.19 9855 0.00E+00 9.384546 0

3971 11.29 10158 0.00E+00 9.485432 0

4071 11.39 10451 0.00E+00 9.578931 0

4171 11.48 10735 0.00E+00 9.67008 0

4271 11.57 11009 0.00E+00 9.75455 0

4371 11.65 11273 0.00E+00 9.836867 0

4471 11.73 11529 0.00E+00 9.913952 0

4571 11.81 11775 0.00E+00 9.985171 0

4671 11.88 12013 0.00E+00 10.05582 0

4771 11.95 12242 0.00E+00 10.12144 0

4871 12.01 12464 0.00E+00 10.18726 0

4971 12.07 12679 0.00E+00 10.24917 0

5071 12.13 12888 0.00E+00 10.30771 0

5150 12.17 13047 0.00E+00 10.35404 0

5150 12.76 13434 156.7 11.02979 3.504364

5171 12.77 13462 157.4 11.03872 3.510807

5271 12.83 13586 160.3 11.07029 3.534707

5371 12.87 13701 163 11.10602 3.558808

5471 12.91 13805 165.4 11.1362 3.579354

5571 12.95 13898 167.6 11.16139 3.597511

5671 12.98 13981 169.6 11.18631 3.614728

5771 13.01 14053 171.3 11.20592 3.628608

5871 13.03 14114 172.7 11.22458 3.640608

5971 13.05 14164 173.9 11.2385 3.650434

6071 13.07 14203 174.9 11.24772 3.658113

6171 13.08 14231 175.5 11.25565 3.662981

6271 13.09 14248 175.9 11.25893 3.665752

6371 13.09 14253 176.1 11.26153 3.667836

51作业君版权所有