5

Reviews in Mineralogy & Geochemistry Vol. 71 pp. 99-128, 2010 Copyright © Mineralogical Society of America

First Principles Quasiharmonic Thermoelasticity of Mantle Minerals Renata M. Wentzcovitch Department of Chemical Engineering and Materials Sciences and Minnesota Supercomputer Institute University of Minnesota Minneapolis, Minnesota, 55455, U.S.A. [email protected]

Zhongqing Wu Department of Physics & Astronomy University of Southern California Los Angeles, California, 90089-0484, U.S.A.

Pierre Carrier Department of Computer Science University of Minnesota Minneapolis, Minnesota, 55455, U.S.A. INTRODUCTION Thermodynamic (Anderson 2005) and elastic properties (Musgrave 1970) of minerals provide the fundamental information needed to analyze seismic observations and to model Earth’s dynamic state. The connection between pressure, temperature, chemical composition, and mineralogy that produce seismic velocity gradients, heterogeneities, and discontinuities, can be established with knowledge of thermoelastic and thermodynamic equilibrium properties of single phases and their aggregates. There is a broad-based need in solid Earth geophysics for these thermoelastic properties to model the Earth. However, the materials and conditions of the Earth’s interior present several challenges. The chemical composition of the Earth’s mantle is complex with at least five major oxide components and tens of solid phases. Today, this type of challenge is more effectively addressed by a combination of experimental and computational methods. Experiments offer accurate information at lower pressures and temperatures, while computations offer more complete and detailed information at higher pressures and temperatures, where experimental uncertainties are large and conditions difficult to control in the laboratory. The overwhelming success of density functional theory (Hohenberg and Kohn 1964; Kohn and Sham 1965) combined with the quasiharmonic approximation (QHA) (Born and Huang 1956; Wallace 1972) for computations of thermodynamic properties of major mantle minerals has been reviewed in this book (Wentzcovitch et al. 2010). Although these properties were cited and compared with experiments only at ambient conditions, in general, quasiharmonic calculations perform even better at higher pressures (see references in Wentzcovitch et al. 2010). In this paper we review the extension of these calculations to high pressure and high temperature elasticity. There are important exceptions, such as CaSiO3-perovskite, whose high temperature structure is stabilized by anharmonic fluctuations (Stixrude et al. 1996). Such 1529-6466/10/0071-0005$05.00

DOI: 10.2138/rmg.2010.71.5

100

Wentzcovitch, Wu, Carrier

cases cannot be addressed using the QHA and one must resort of molecular dynamics (MD) in some form. Great progress has been made in the last decade in this area of first principles computations of elastic constants, whether by MD (Oganov et al. 2001; Li et al. 2005, 2006b,c, 2009; Stackhouse et al. 2005) or QHA calculations (Karki et al. 1999; Wentzcovitch et al. 2004, 2006). This progress resulted from the maturing of theoretical and computational methods (Car and Parrinello 1985; Giannozzi et al. 1991; Wentzcovitch and Martins 1991; Payne et al. 1992; Baroni et al. 2001) and from the availability of stable and well-tested public domain software (Payne et al. 1992; Kresse and Furthmüller 1996; Gonze et al. 2002; Giannozzi et al. 2009) for materials computations. The purpose of our effort has been to provide the most significant and accurate information needed to understand the fine- and large-scale structure of the mantle. A wide variety of geophysical issues of fundamental importance remain to be resolved, e.g., the geodynamic implications of upper mantle anisotropy, the causes and significance of an upper mantle low velocity zone, the significance of seismic discontinuities and high velocity gradients in the transition zone for possible mantle stratification, and whether heterogeneity in the lower mantle is primarily of chemical/mineralogical versus thermal in origin. Experimental elasticity and thermodynamic data for upper mantle and transition zone minerals are more numerous. The pressures and temperatures in the upper mantle (including the transition zone) can be more closely approached experimentally and with relatively high accuracy. While it is possible to produce extremely high pressures and temperatures in the laser heated diamond anvil cell (DAC), elasticity of relevant phases at these conditions are less well established experimentally for various reasons, including uncertainties in pressure scales, thermal gradients across laserheated samples, and other technical difficulties in characterizing the most extreme pressures and temperatures. Therefore, the emphasis of our effort has been on lower mantle phases at extreme conditions. This is the pressures and temperatures regime where theory can make the greatest impact in augmenting experimental data and producing essential information.

THEORETICAL BACKGROUND Several fundamental formal papers on elasticity theory have been written throughout the years (Barron and Klein 1965; Thomsen 1972; Davies 1974). There are also excellent textbooks on more fundamental aspects and one usually starts research on elasticity by studying them (Musgrave 1970; Nye 1985). Here we present a distilled version of the basic formalism for calculations of thermoelastic properties of materials. Because our calculations are performed numerically and access directly the free energy computed by first principles, the need for analytical treatment and approximations is minimal. The tricks are in the numerical implementation, not in the analytical development. The simplicity of the basic formalism contrasts with elaborate formal developments addressing analytically the consequences of the free energy expansion (Thomsen 1972; Davies 1974), by necessity truncated, in a power series of strains. Formal developments originated in the practical impossibility of measurements to access directly the free energy, but only its derivative with respect to volume, i.e., pressure, or the elastic constants via velocity measurements. Reviews of static (Karki et al. 2001a) and semi-empirical (Stixrude and Lithgow-Bertelloni 2005) approaches to elasticity of minerals discuss several useful details of elasticity calculations. Here we focus on details of high-pressure and high-temperature QHA calculations that have not been explicitly presented yet. We treat a crystal as a homogeneous, anisotropic medium and assume that stress and strain are anisotropic. The propagation of acoustic waves in an aeolotropic medium is described by Cristoffel’s equation (Musgrave 1970): [Cijkl n j nl − ρv 2 δik ] pk = 0

(1)

Thermoelasticity of Mantle Minerals

101

Cijkl are the components of the (4th rank) elastic constant tensor in Cartesian coordinates, with i, j, k, and l as one of the Cartesian indexes x, y, or z, nj’s are the direction cosines of the propagation direction, and pk’s are the components of the polarization vector. This equation is true at any temperature or pressure and the eigen-solutions are three waves, one quasi-longitudinal and two quasi-transversal. Purely longitudinal and transversal polarizations occur only in isotropic materials. We will return to this point at the end of this section. These calculations of Cijkl require the definition of a finite strain induced by a large compression and of an additional superposed infinitesimal strain. The large finite strain of choice is the frame invariant Eulerian strain Eij : Eij =

1 ( fij + f ji − fik f jk ) 2

(2)

with Xi − ai = fij X j (3)   where a is the position vector of a point in a material in the unstrained state and X is the position vector of the same point after the material is compressed. Eij is invariably assumed to be isotropic in finite strain expansions of the free energy used to derive equations of state. This is assumed here as well. Eij = − f δij f =

1  V0   2  V 

2/3

(4a )  − 1 

( 4 b)

where V0 is a reference volume and f > 0 (< 0) corresponds to compression (extension). The Helmholtz free energy of a compressed state is then expanded in a Taylor series in these strains. F ( f ) = a + bf + cf 2 + df 3 + ...

(5)

This expansion includes a linear term, meaning, the reference state with V = V0 is not the zero pressure state. The temperature dependence of this free energy is temporarily omitted. The second infinitesimal strain, uij, is a Lagrangian strain associated with the deformation of the pre-compressed state: xi − Xi = uij X j

(6)

The choice of Lagrangian versus Eulerian strain is indifferent. Lagrangian strains seem more natural in practical calculations. Its frame invariant form is: eij =

1 (uij + u ji − uiku jk ) 2

(7a )

1 (uij + u ji ) 2

(7b)

or simply eij =

Superposition of a deformation on pre-compressed states changes Equation (5) into: F ( f , eij ) = a(eij ) + b(eij ) f + c(eij ) f 2 + d (eij ) f 3 + ...

(8)

It is implicit in this formula that volume is expressed also in terms of two independent strains, f and eij. Therefore, V0 in Equation (4b) is V0(0,eij).

Wentzcovitch, Wu, Carrier

102

The isothermal stress-strain coefficient tensor at a certain pressure and temperature T Cijkl ( P, T ) is: T C ijkl ( P, T ) =

1  ∂ 2G  V  ∂eij ∂ekl

    P ,T

(9a )

with (in contracted notation) P(V , T ) T = −

∂F (V , T ) ∂f ∂f ∂V

− eij ,T

∂F (V , T ) ∂eij ∂eij ∂V

(9b) f ,T

For volume conserving combinations of eij, the second term on the right-hand side of Equation (9b) essentially vanishes. The first term is essentially independent of eij for any infinitesimal strain. Then, Equation (9a) becomes: T C ijkl ( P, T ) =

1  ∂2F  V  ∂eij ∂ekl

  + Pδijkl   f ,T

(10a )

with δijkl =

1  ∂ 2V  V  ∂eij ∂ekl

 1  = ( 2δij δkl − δik δ jl − δil δ jk )   f ,T 2

(10 b)

For these strains, one can also write changes in the free energy in the following convenient form (Wallace 1972): δF V

(11a )

f ,T

1 T = σij eij + C ijkl eij ekl 2

δG V

1 T eij ekl = C ijkl 2

(11b)

or, because of Equations (9a,b):

P ,T

Although the natural variables of the Gibbs free energy are temperature, T, and stresses, σij (see Stixrude and Lithgow-Bertelloni 2005), the use of strains, eij, is possible as long as one can calculate this free energy. This is possible because pressure, i.e., the first term in the righthand side of Equation (9b) is basically independent of eij. Our implementation of quasiharmonic elasticity calculations proceeds as follows: a) equilibrate the crystal structure at several pressures (several f ’s) (Wentzcovitch et al. 1993); b) apply a series of combinations of strains, eij, positive and negative, as many as necessary, for each equilibrium structure and relax the atomic coordinates with fixed cell shape; c) compute the vibrational density of states for equilibrium and strained configurations (Baroni et al. 2001); d) the quasiharmonic free energy is then computed for all configurations generated: F ( f , ei , j , T ) = E ( f , ei , j ) +

  ω ( f , e )   1 ∑ ωq ( f , ei, j ) + kBT ∑q ln  1 − exp  − qk T i, j   2 q B   

(12)

where E(f,eij) is the static total energy from first principles calculations; e) fit Equation 8 for each strain eij and temperature T. The temperature grid is as dense as desired. The pressure grid (read compressive Eulerian strain grid) is sparse. For example, to cover the entire pressure range of the mantle, up to ~135 GPa, the calculation in MgSiO3- perovskite required 10-15 pressures −20 GPa < P < 200 GPa. This enlarged pressure interval warrants smoothness of

Thermoelasticity of Mantle Minerals

103

free energy derivatives at the end of the desired pressure interval and provides information on expanded volumes necessary to describe the state of the system at low pressures and high temperatures. Once these free energies are available one can proceed with the finite difference calculation of the second derivatives in Equation (9a) or Equation (10). Thermodynamic properties are obtained simultaneously from F(V(f),T). This procedure for calculating high-P, T elasticity is analogous to that of static elasticity calculations (Wentzcovich et al. 1995; Karki et al. 1997, 2001a), although in static calculations one obtains the first principles stress tensor (Nielsen and Martin 1985) directly. High temperature isothermal stresses can also be calculated (Davies 1974): σTij ( P, T ) =

1 ∂F V ∂eij

(13) f ,T

The isothermal constants given by Equation (9a) are relevant for comparisons with measurements in static compression experiments only. The time scale of deformation in seismic events, or in laboratory measurements of acoustic wave velocities using Brillouin scattering or resonant ultra-sonic spectroscopy, is much shorter than that of thermal diffusion for relevant length scales. Therefore, one needs to compute instead the adiabatic elastic constants. Standard algebraic manipulations involving changes of variables give: S T C ijkl ( P, T ) = C ijkl ( P, T ) +

VT λ ij λ kl CV

(14a )

with  ∂S ( P, T )  λ ij ( P, T ) =    ∂eij  f ,T

(14b)

where S(P,T) is the entropy. Similarly, adiabatic stresses become: σ Sij ( f , T ) =

1 ∂U V ∂εij

= σTij ( f , T ) + S, f

ST λ ij CV

(15)

In Cartesian notation, there are 81 independent elastic constants, but this number is reduced to 21 by the requirement that Cijkl are symmetric with respect to interchanges (i,j), (k,l), and (ij,kl). This allows the replacement of a pair of Cartesian indices ij by a single index α, the Voigt index, according to: 11 → 1, 22 → 2, 33 → 3, 23 or 32 → 4, 13 or 31 → 5, 12 or 21 → 6. Simultaneously, one should also replace Cartesian strains by Voigt strains, εα, such that: eii = εα, eij = eji = εα/2, and σij = σji = σα, with the relationship between ij → α shown above. With this replacement, one can re-write Equations (11a,b) using also Voigt indices varying from 1 to 6: δF V

(16a )

f ,T

1 T = σα εα + Cαβ ε α εβ 2

δG V

1 T = Cαβ ε α εβ 2

(16 b)

and

P ,T

We are interested also in computing elastic properties of polycrystalline aggregates. This is an extensive subject thoroughly discussed in textbooks dedicated to elasticity (Musgrave 1970). Besides, we follow a treatment of this topic that is standard and extensively used in mineral physics. Therefore, we simply mention it here.

104

Wentzcovitch, Wu, Carrier

A purely isotropic material is viewed as a polycrystalline aggregate without preferred orientation of grains. The average elastic constants of this material are defined by the relationship between macroscopic stresses and strains. There are basically two ways to think about and calculate this relationship: i) assuming the stress is uniform across grains and strain is not or ii) assuming strain is uniform and stress is not. The uniform stress case corresponds to a state with internal relaxation at grain boundaries and gives lower bounds for the elastic constants (Musgrave 1970). This is the so-called Reuss average. The second case is the Voigt average and gives upper bounds. The average of the two averages is the Hill or VoigtReuss-Hill (VRH) average (Watt et al. 1976). The Voigt and the Reuss averages consist in computing the orientational averages of the elastic constant tensor, C, and compliance tensor, S = C−1, respectively. The average elastic constant tensor so computed, contains only two independent elastic constants with their respective bounds: KVRH ± δKVRH and GVRH ± δGVRH, with δMVRH = ½(MV – MR). Elastic constant tensors with preferred lattice orientations can also be calculated this way. For a detailed presentation of this topic, see Musgrave (1970). Isotropic elastic constants with tighter bounds can be obtained using the Hashin-Shtrickman variational principle (Hashin and Strikman 1962). Results for several crystalline systems have been summarized by Watt (1980). Finally, in a isotropic medium there are two distinct wave velocities obtained from Crystoffel’s equation: the longitudinal or primary velocity, VP , and the twofold degenerate shear velocity, VS: VP =

4 K+ G 3 ρ

and

VS =

G ρ

(17)

from which the bulk sound velocity can be defined: Vϕ =

K 4 = VP2 − VS2 ρ 3

(18)

These are the basic ingredients of our calculations. We now summarize results on the high temperature elastic properties of the major lower mantle phases, MgO, MgSiO3-perovskite, and MgSiO3-post-perovskite, and some geophysical consequences derived from these results. Although there is by now a considerable number of papers on high temperature elasticity of mantle phases computed using MD (Oganov et al. 2001; Li et al. 2005, 2006b,c, 2009; Stackhouse et al. 2005), here we restrict ourselves to review quasiharmonic calculations only. These calculations offer elastic properties in a continuum of pressure and temperature, which is important when computing gradients of these properties. This is quite common in geophysics, particularly when trying to distinguish thermal versus compositional effects, as will be pointed out below.

ELASTICITY OF LOWER MANTLE PHASES Earth’s internal structure, temperature, and mantle phases have been summarized in the previous paper (see Figs. 1, 2, and 3 in Wentzcovitch et al. 2010). We have been primarily concerned with the elastic properties of lower mantle phases: (Mg,Fe,Al)(Si,Al)O3-perovskite, (Mg,Fe,Al)(Si,Al)O3-post-perovskite, (Mg,Fe)O, and CaSiO3-perovskite. Here we review the high temperature quasiharmonic elasticity of MgSiO3-perovskite, MgSiO3-post-perovskite, and MgO. Some properties of the (Mg,Fe) phases have been reviewed by Hsu et al. (2010) in this volume. Static elastic properties of the perovskite and post-perovskite containing aluminum and iron have been reported by Caracas and Cohen (2007). High temperature properties of CaSiO3perovskite cannot be investigated using the QHA. It is a very anharmonic phase and its high temperature structure is stabilized by anharmonic fluctuations (Stixrude et al. 1996; Caracas et

Thermoelasticity of Mantle Minerals

105

al. 2005). Its elastic properties have been investigated by molecular dynamics (Li et al. 2005) and by using a free energy expansion in terms of its structural parameters (Stixrude et al. 2007). In first approximation Earth is spherical and homogeneous but there are important deviations of both. Seismic tomography (Woodhouse and Dziewonski 1984; Grand 1994; van der Hilst et al. 1997; Masters et al. 2000) reveals the 3D velocity structure of the mantle, local and global. The starting point for understanding these observations is knowledge of individual phases’ elastic properties at in situ conditions. These calculations have used the local density approximation (LDA) (Ceperley and Alder 1980; Perdew and Zunger 1981) combined with quasiharmonic theory (Born and Huang 1956; Wallace 1972) and density functional perturbation theory (DFPT) (Baroni et al. 2001) for phonon calculations. Variable cell shape MD is used for structural equilibration at desired pressures (Wentzcovitch 1991). The same pseudopotentials of Troullier and Martins (1991), Vanderbilt (1990) and Von Barth and Car (unpublished), cited in the appendix of Wentzcovitch et al. (2010), were used here. Computations were performed using codes of the Quantum ESPRESSO distribution (Giannozzi et al. 2009).

MgO Periclase (MgO) was the first system to have its high temperature elastic properties investigated by this fully first principles approach (Karki et al. 1999). The original results, including some small discrepancies from experiments, have been reproduced since then. MgO is the most anisotropic phase of the lower mantle (Wentzcovitch et al. 2006) and its anisotropy is strongly pressure dependent. Its thermodynamic properties (Karki et al. 2000b), including anharmonic effects (Wu et al. 2008; Wu and Wentzcovitch 2009) were summarized in Wentzcovitch et al. (2010). Although temperature effects on its elasticity at ambient conditions are known to be substantial and to counteract the effect of pressure (Isaak et al. 1989; Chen et al. 1998), elasticity at deep lower mantle conditions have not been measured so far. A posteriori inspection of the quasiharmonic thermal expansivity (see Wentzcovitch et al. 2010) indicates that this is a good approximation for MgO at 0 GPa up to ~1,100 K and even better at geophysically relevant conditions (Wu and Wentzcovitch 2009). The elastic constants of this materials with rocksalt structure (see Wentzcovitch et al. 2010), C11, C12, and C44, were obtained up to 150 GPa and 3,000 K (Karki et al. 1999) by calculating the free energies for strained lattices. For cubic systems, volume conserving strains (tetragonal and trigonal) allow the determination of C44 and of Cs = (C11 − C12)/2, for which isothermal and adiabatic values are identical. Isothermal or adiabatic values of the bulk modulus, K = (C11 + 2 C12)/3, are then used together to obtain the corresponding values of C11 and C12. The predicted ambient values and their initial pressure and temperature gradients agree quite well with measurements (Isaak et al. 1989; Sinogeikin and Bass 1999) (see Table 1, Figs. 1 and 2). As seen in Table 1, except for C11, the predicted Cij at ambient conditions are smaller than the experimental values. This is typical of DFT calculations that overestimate the equilibrium volume (see Table 1), in this case ~0.5%. So far, LDA has proved to be the best functional for calculations of structural and elastic properties of minerals at finite temperatures (Wentzcovitch et al. 2010) but the small volume overestimation causes the elastic constants to be underestimated. C44 in particular is the most underestimated. The reason is unclear but this result is reproduced still today. However, the temperature dependence of C44 is best reproduced (see Fig. 1). The temperature gradients of Cij at 0 GPa are well predicted for all Cij, but above ~1,250 K the QHA gradients, especially of C11 and C12 start deviating from experimental values. As discussed in Wentzcovitch et al. (2010) this is caused by unharmonic effects. The predicted cross P-T derivatives of Cij are smaller (with opposite sign) than those obtained by experiments on MgO to 8 GPa and 1600 K (Chen et al. 1998). However, QHA results are consistent with earlier data to 0.8 GPa and 800 K (Spetzler 1970) (Table 1) and potential-induced breathing (PIB) model calculations (Isaak et al. 1990).

Wentzcovitch, Wu, Carrier

106

Table 1. Adiabatic elastic moduli (M) of MgO and their pressure and temperature derivatives at ambient conditions. The modulus c110 corresponds to elastic modulus for the longitudinal wave along [110] direction. Numbers in parentheses are experimental uncertainties. From Karki et al. (1999). c12

c11

c44

c110

KS

G

LDA+QHA (V0 = 18.81Å3/cell) M (GPa)

298

94

147

344

162

128

∂M/∂P

9.56

1.45

1.03

6.39

4.15

2.44

−0.0598

0.0089

−0.0088

−0.0343

−0.0140

−0.0216

0.56

−0.06

0.20

0.45

0.14

0.44

297.9(15)

95.8(10)

154.4(20)

351.3(22)

163.2(10)

130.2(10)

9.05(20)

1.34(15)

0.84(20)

6.04(37)

4.0(1)

2.4(1)

−0.0585

0.0075

−0.0126

−0.0381

−0.0145

−0.024

−1.3(4)

5.1(24)

−0.2(3)

1.7(7)

3.0(15)

−1.8(10)

0.1(4)

0.1(3)

0.1(1)

0.0(1)

0.1(3)

0.l(2)

∂M/∂T (GPa/K) ∂M2/∂P∂T (10−3/K)

Experiments (Vexp = 18.69Å3/cell) M (GPa) * ∂M/∂P * ∂M/∂T (GPa/K) + ∂M /∂P∂T (10 /K) 2

*

#

-3

$

Sinogeikin et al. (1998), + Isaak et al. (1989), # Chen et al. (1998), $ Spetzler (1970)

The pressure-induced change of sign in the single-crystal anisotropy, as expressed by the anisotropy factor: A=

(2C44 − C11 + C12 ) C 11

(19)

is characteristic of MgO (Sinogeikin and Bass 1999; Duffy et al. 1995). This quantity expresses the breakdown of the Cauchy relation that holds for isotropic materials. Our results show that the strong pressure dependence of the anisotropy in MgO is preserved at high temperatures (Fig. 3). Temperature effects counteract pressure effects and are monotonically suppressed with increasing pressure. MgO is the most anisotropic phase of the lower mantle at core-mantle boundary (CMB) conditions, as will be seen below. Therefore, lattice-preferred orientation (LPO) in MgO is a possible cause of anisotropy observed in this region. The isotropic longitudinal (VP) and shear (VS) wave velocities of MgO along several isotherms are shown in Figure 4. The pressure dependence of VS agrees very well with 300 K experiments up to ~100 GPa (Murakami et al. 2009). In the lower mantle, the properties of MgO are modified by the presence of iron (see Hsu et al. 2010). The wave velocities of ferropericlase, (Mg,Fe)O, are smaller and the longitudinal velocity is affected by the high spin to low spin crossover of iron at lower mantle pressures (Wentzcovitch et al. 2009; Hsu et al. 2010). In fact all thermodynamic properties are affected by this crossover (Wu et al. 2009). The seismic parameters: RS / P =

∂(ln Vφ ) ∂(ln VS ) ∂(ln ρ) , Rϕ / S = , and Rρ / S = ∂(ln VP ) P ∂(ln VS ) P ∂(ln VS ) P

(20)

express the relative magnitude of lateral (isobaric) variations in VS, VP , Vφ (Eqn. 17 and 18), and density ρ in the mantle. Seismic tomography (Dziewonski and Woodhouse 1987; Grand 1994; Masters et al. 2000), as summarized by Karato and Karki (2001), indicates that RS/P increases

Figure 1. Temperature dependence of the elastic moduli of MgO at 0 GPa. Experimental data (Isaak et al. 1989) are shown as open circles.

Figure 2. Pressure dependence of the elastic moduli of MgO along various isotherms. Experimental data are shown as solid squares (Sinogeikin and Bass 1999).

Thermoelasticity of Mantle Minerals 107

Figure 3. Pressure dependence of calculated anisotropy factor of MgO along 300 (lowest curve), 1,000, 2,000, and 3,000 (uppermost curve) K isotherms. Dashed lines correspond to lowpressure (18.6 GPa) data at ambient temperature (Sinogeikin and Bass 1999). The zero pressure data at 300, 1000, and 1800 K are indicated by diamonds (Isaak et al. 1989). From Karki et al. 1999.

Figure 4. Pressure dependence of velocities of MgO along various isotherms. Experimental data are denoted by solid triangles (Li et al. 2006a), open circle (Sinogeikin and Bass 1999), solid circles (Murakami et al. 2009), and dashed lines (Zha et al. 2000), respectively.

108 Wentzcovitch, Wu, Carrier

Thermoelasticity of Mantle Minerals

109

from more than ~2.0 to more than 3.0 from the top to the bottom of the lower mantle. Possible causes of lateral velocity variations include variations of temperature, composition, phase, etc. High temperature experiments at ambient pressure have yielded RS/P ~1.5 for thermally induced v (Isaak et al. 1989). These QHA velocities constrain “thermal” RS/P in MgO to vary from ~1.4 at the top to ~1.9 at the bottom of the lower mantle. Earlier PIB calculations estimated RS/P ~ 2.5 in this region (Isaak et al. 1990). As will be seen below, lateral temperature change in the bottom of the mantle can induce the post-perovskite transition (Tsuchiya et al. 2004; Oganov and Ono 2004; Wentzcovitch et al. 2010). Lateral temperature and phase change can increase considerably this parameter (Wentzcovitch et al. 2006). For a more extensive discussion of lateral variations in the mantle see Stixrude and Lithgow-Bertelloni (2007).

MgSiO3-perovskite Owing to its abundance in the Earth lower mantle (LM) (see Fig. 1 and 2 in Wentzcovitch et al. 2010) the elastic properties of MgSiO3-perovskite (perovskite henceforth), more precisely (Mg,Fe)SiO3, determine to a great extent the properties of this region. In situ conditions, with temperatures varying from ~1,900 K to 4,000 K, are still challenging for elasticity measurements, even though elasticity measurements in the entire pressure range of the mantle is now possible (Murakami et al. 2009). Therefore, these first principles calculations have been essential to advance understanding of the state of the lower mantle. Single crystal and aggregate elastic properties of perovskite and of other mantle phases have been used to predict acoustic velocities of hypothetical lower mantle aggregates. Comparison with the spherically averaged seismic velocity profiles of this region (Dziewonski and Anderson 1981) have then made to test the plausibility of homogeneous mineralogical models (Stixrude et al. 1992; da Silva et al. 2000; Karki et al. 2001b; Wentzcovitch et al. 2004; Xu et al. 2008). These comparisons aim to resolve, for instance, the problem of compositional stratification in the mantle, with implications for the style of thermochemical convection that has been operating in the Earth. Perovskite has orthorhombic (Pbnm) symmetry and 9 independent elastic constants (Cij henceforth), 3 diagonal (C11, C22, C33), 3 off-diagonal (C12, C13, C 23), and 3 shear (C44, C55, C66) Cij. Its adiabatic Cij at ambient conditions are summarized and compared to experimental data in Table 2. As can be seen, the LDA/QHA (Wentzcovitch et al. 2004) equilibrium volume at ambient conditions is ~1.2% larger than the experimental volume (Karki et al. 2000b). This “expanded” volume causes substantial underestimation of Cij. The thermodynamics properties of perovskite were summarized in Wentzcovitch et al. (2010) and the source of this discrepancy does not appear to be the QHA. The generalized gradient approximation (GGA) (Wang and Perdew 1991; Perdew et al. 1996) overestimates this volume even further (Carrier et al. 2007) (on these two issues, see discussion in Wentzcovitch et al. 2010). However, this discrepancy should decrease with pressure and these predictions should be more reliable, especially predictions of pressure and temperature gradients. Cij at LM conditions are shown in Figure 5 compared with results of MD simulations using the GGA (Oganov et al. 2001). Both sets of results “shifted” pressure to bring ambient condition values into agreement with experimental values. The overall good agreement between these results is more than coincidence. Pressure shift to correct problems with DFT results is not a good practice (Wu et al. 2008) since the problem caused by DFT decreases with pressure (increasing density). A better perspective of the discrepancy between calculations and measurements of aggregated elasticity and velocities is shown in Figure 6, where no shift is applied to QHA/LDA results. In the overall scale of changes in Cij caused by pressure and temperature, as well as experimental uncertainties, the discrepancy is relatively small, and the accuracy of these results is still useful. Without losing perspective of these errors and uncertainties, it is still possible to extract sensible information about the state of the lower mantle. Figure 7 displays the adiabatic bulk and shear moduli, KS and G, for isotropic and homogeneous mineralogical models along a standard adiabatic geotherm (Brown and Shankland 1981) that nearly coincides with the isentropes

0.0009

0.0015

0.0003

0.0002

0.0004

−0.016

0.0004

−0.002

0.0001

0.002

0.00

0.004

−0.001

−0.006

0.0006

0.0012

0.001

−0.002

(4) ∂M2/∂P∂T

−0.014

−0.018

3.11

5.07

−0.021

−0.020

−0.030

2.3

3.7

264, 253*,#

(5) ∂M/∂T

−0.031

−0.009

2.3

3.4

146

−0.015

−0.025

3.2

4.5

147

(4′) ∂M/∂T

−0.073

1.4

2.7

144

−0.027

−0.072

0.7

2.0

147

−0.058

1.3

2.5

186

(4) ∂M/∂T

(4′) ∂M/∂P

5.0

9.4

204

4.4

5.0

3.1

485

(5) ∂M/∂P

8.7

7.6

(4) ∂M/∂P

537

482

(3) M (162.3Å3)

0.0005

−0.028

−0.01

−0.0231

2.0

1.0

2.1

177, 173*, 175+,#

0.00002

−0.0001

−0.0003

0.020

0.07

11.04, 10.88# 10.86*

0.00001

−0.0001

−0.0002

0.003

0.028

6.57, 6.56+, 6.49*,6.53#

6.51

(1) M (164.1Å )

(2) M

10.89

496

560

444

491

…

175

428

474

c22

… 257

185 203

c33

… 144

165 186

c44

… 135

138 153

c55

… 128

121 134

c66

… 145

127 144

c12

…

173

142 156

c13

…

198

238 263

c23

…

456

162 179

KS

…

524

6.32 6.53

10.58 10.94

G

…

477

VS

VP

…

3

…

c11

Table 2. Adiabatic elastic constants, their pressure and temperature derivatives, and wave velocities of MgSiO3-perovskite predicted by LDA+QHA: (1) results at ambient conditions (Wentzcovitch et al. 2004); equilibrium volume is listed in parentheses in first column; (2) range of values obtained in static density functional theory calculations (Wentzcovitch et al. 1998; Kiefer et al. 2002); (3) experimental values obtained by Brillouin (Yeganeh-Haeri 1994), ultrasonic measurements (+) (Sinelnikov et al. 1998), (*) (Li and Zhang 2005), and (#) (Sinogeikin et al. 2004) the experimental equilibrium volume is listed at parentheses of first column; (4) at ambient conditions and (4′) at 300 K and 100 GPa; (5) experimental values at ambient conditions obtained by ultrasonic measurements (Li and Zhang 2005). From Wentzcovitch et al. (2004).

110 Wentzcovitch, Wu, Carrier

Figure 6. Pressure dependence of (a) bulk modulus and shear modulus, (b) the isotropic longitudinal and shear wave velocities of MgSiO3. The experimental data are from Sinelnikov et al. (1998) (empty quinquangular), Sinogeikin et al. (2004) (plus), Yeganeh-Haeri et al. (1994) (empty up triangles), Li and Zhang (2005) (open circles) and Jackson et al. (2005) (solid squares). Symbols connected by dashed (1,500 K) and solid (3,500 K) line are molecular dynamic simulation results of Oganov et al. (2001) (empty symbols) and Marton and Cohen (2002) (solid symbols).

Figure 5. Pressure dependence of the adiabatic elastic constants of MgSiO3. The line’s thickness represents uncertainties caused mainly by the use of the local density approximation (LDA). Full (dashed) lines correspond to results within (outside) the PT regime of validity of the QHA. Measurements are represented by full circles at 0 GPa (Yeganeh-Haeri et al. 1994). Full symbols at P > 0 GPa are the results of Oganov et al. (2001) at 38 and 88 GPa (squares: 1500 K, diamonds: 2500 K, up triangles: 3500 K). From Wentzcovitch et al. (2004).

Thermoelasticity of Mantle Minerals 111

112

Wentzcovitch, Wu, Carrier

Figure 7. Bulk (KS) and shear (G) moduli for (Mg(1−x)Fex)SiO3-perovskite and pyrolite with 20% < Vfp < 30% (ferropericlase = Mg(1−y) FeyO) along the Brown and Shankland (1981) geotherm. 0.0 < x < 0.12 and 1 < y/x < 4. The effect of iron on the elastic moduli of perovskite and ferropericlase were extracted from Kiefer et al. (2002) and Jackson (1998), respectively. From Wentzcovitch et al. (2004).

of MgSiO3 and MgO (Karki et al. 2000a,b). The mineralogical models consist of pure perovskite, (Mg(1-x)Fex)SiO3, with x~0.12 and pyrolite (Jackson 1998), which contains ~20 vol% of ferropericlase, (Mg(1-y)Fey)O, with y~0.18. Iron partitioning (y/x) between these phases was allowed to vary in the range 1
Thermoelasticity of Mantle Minerals

113

of the possibilities. Indeed, we know today that there is a new phase present in the deep mantle, i.e., MgSiO3-post-perovskite (to be discussed next). The shear modulus of this phase is larger than that of perovskite by ~20 GPa at D′′ conditions, which extends up to ~300 km above the CMB. PREM does not show the discontinuity most likely associated with this phase change. It has wide topography (~300 km wide) and the positive velocity jump in G associated with this phase change is smoothed by PREM. This topic is reviewed next.

MgSiO3-post-perovskite Post-perovskite is almost unanimously accepted to be the predominant phase of MgSiO3 just above the core mantle boundary (CMB), the D′′ region. Its discovery (Murakami et al. 2004; Oganov and Ono 2004; Tsuchiya et al. 2004) has had multidisciplinary impact in geophysics. High temperature elasticity calculations of this phase produced revealing results that, when compared with seismic data, suggest its presence in the D′′ layer. Therefore, in this section will make reference to geophysical issues more than in any other. The importance of these calculations to geophysics will be abundantly demonstrated here. Post-perovskite has orthorhombic structure (Cmcm) and 9 independent Cij, like perovskite (Tsuchiya et al. 2004; Wentzcovitch et al. 2006). They are quite different from perovskite’s (see Table 3). Post-perovskite has a layered structure (see Fig. 6b in Wentzcovitch et al. 2010), expands anisotropically, and has complex pressure- and temperature-dependent elastic behavior. However, its aggregate moduli do not differ so much from those of perovskite. Thermodynamic properties of both phases, KS, are quite similar beyond ~80 GPa (Tsuchiya et al. 2005). The shear modulus of post-perovskite, G (Fig. 8a), is larger and has larger pressure and temperature gradients. This result and the fact that post-perovskite is ~1.5% denser than perovskite at deep lower mantle conditions determine the velocity contrasts between these phases. The shear velocity, VS, of post-perovskite and its gradients (Fig. 8b) are larger than those of perovskite (Fig. 6b). Because of G, VS of post-perovskite is also larger and has larger gradients. In contrast, the bulk velocity, Vφ (Eqn. 18), is smaller than that of perovskite because of their similar KS and post-perovskite’s larger density. Figure 8c shows the velocity jumps across the calculated phase boundary (see Fig. 21 in Wentzcovitch et al. 2010) (see Fig. 8c caption). These results are consistent with the increase in seismic velocities observed ~200–300 km above the core–mantle boundary in certain places (Lay and Helmberger 1983) but most easily detected beneath regions of past subduction, presumably colder places, such as beneath Central America (Lay and Helmberger 1983; Wyssession et al. 1998; Garnero 2004). There, VS > VP and ΔVS ~ 2–3% is observed, but this observation is clearly a regional property of a notably heterogeneous layer. Another baffling property of D′′ revealed by global tomographic models (Woodhouse and Dziewonski 1984; Grand 1994; Masters et al. 2000) is the anti-correlation between lateral (isobaric) heterogeneities in Vφ and VS. The likely causes are usually addressed by comparing the seismic parameters, RS/P, Rφ/S, and Rρ/S (see Eqn. 20) to theoretical or experimental predictions of these ratios at relevant conditions (Karato and Karki 2001). It is known that in the shallow lower mantle RS/P ~ 2.3 and Rφ/S ~ 0.0, whereas in D′′ RS/P ~ 3.4 and Rφ/S ~ −0.2 (Masters et al. 2000). Velocity anomaly ratios produced by isobaric temperature changes in pure post-perovskite and perovskite aggregates are displayed in Figure 9a,b along with the seismic parameters extracted from Karato and Karki (2001). In pure perovskite aggregates, RS/P increases with pressure and temperature but reaches at most ~2.3 at 135 GPa and 4,000 K, whereas Rφ/S is approximately pressure independent and slightly decreases with temperature to reach ~0.16 at similar conditions (see Fig. 2a). In post-perovskite RS/P decreases with pressure, but because of its larger (∂G/∂T)P, it increases more rapidly with temperature to reach ~ 2.8 at the same conditions. Rφ/S is smaller than that of perovskite, ~0.1, at these conditions (see Fig. 9a,b). It has been argued that anelasticity, anisotropy, and lateral variations in calcium content in the deepest mantle might be necessary to produce these large RS/P and negative values for Rφ/S (Karato and Karki 2001). Figure 9c compares the seismic parameters with the computed ratios of velocity anomalies

255

238

231

2.2

−0.028

0.00030

(2) M

(3) M

(4) ∂M/∂P

(5) ∂M/∂T

(6) ∂2M/∂P∂T

c66

c55

(1) M

0.00017

0.00061

(6) ∂2M/∂P∂T

4.5

0.00014

−0.045

2.7

326

352

296

−0.037

6.4

900

−0.083

(3) M

888

1095

(4) ∂M/∂P

1119

(2) M

c22

(5) ∂M/∂T

874

1146

(1) M

c11

0.000063

−0.017

3.6

685

656

655

BS

0.00026

−0.069

6.3

1131

1139

1077

c33

0.00032

−0.030

2.2

284

294

276

G

-0.000063

0.0011

2.9

498

454

539

c12

0.0000035

−0.00029

0.032

14.0

14.0

13.9

VP

0.000089

0.022

2.3

486

418

436

c13

0.0000053

−0.0003

0.022

7.3

7.4

7.2

VS

−0.00013

−0.0054

2.5

536

507

469

c23

0.0000012

−0.000073

0.022

11.3

11.1

11.1

VΦ

0.00076

−0.011

3.1

343

311

311

c44

Table 3. Elastic properties in GPa, km/s, and K of MgSiO3 (1) perovskite and post-perovskite ((2), (4), (5), (6)) at 125 GPa and 2500 K and (3) 140 GPa and 4000 K. Diagonal cijs are typically underestimated by ~2%, off diagonal ones by ~1.5%, and shear ones by less than 1%. BS and G are Voigt-Reuss-Hill averaged adiabatic bulk and shear moduli. Typically they are underestimated by 2% (random deviations in individual cij are averaged out). VP, VS and VΦ are compressional, shear, and bulk velocities. [From Wentzcovitch et al. 2006.]

114 Wentzcovitch, Wu, Carrier

Thermoelasticity of Mantle Minerals

115

Figure 8. Pressure dependence of (A) aggregate adiabatic bulk (BS) and shear (G) moduli and density (ρ) of post-perovskite along the 300 K (uppermost), 1000 K (second from top), 2000 K (thrid from top), 3000 K (fourth from top) and 4000 K (fifth from top) isotherms; (B) isotropic longitudinal (VP = [(BS + 4G/3)/ρ]1/2), shear (VS = (G/ρ)1/2) and bulk (VS = (B/ρ)1/2) wave velocities. Filled circles in (A,B) are PREM values (Dziewonski and Anderson 1981) for comparison; (C) velocity jumps across the previously obtained phase boundary (Tsuchiya et al. 2004). Thick black lines represent the jumps at the center of our phase boundary with a DFT related uncertainty of 10 GPa. Shaded areas are possible values throughout the boundary uncertainty domain (Tsuchiya et al. 2004)). Open triangles and squares are respectively results from MD simulations at 135 GPa and 4,000K and 136 GPa and 3,000 K (Stackhouse et al. 2005). The vertical dashed line indicates approximately the topmost location of the D′′ discontinuity. From Wentzcovitch et al. (2006).

116

Wentzcovitch, Wu, Carrier

Figure 9. Seismic parameters RS/P = (∂lnVS/∂lnVP)P, Rφ/S= (∂lnVφ/∂lnVP)P, and Rρ/S = (∂lnρ/∂lnVS)P at several temperatures for post-perovskite (A) and perovskite (B) (sequence same as in Fig. 8 but now starting from 1000 K). These values are presented for comparison. Filled circles with error bar (± 20%) are from MD simulations (Stackhouse et al. 2005). (C) Contribution of phase transition to these parameters with uncertainties (shaded blue areas) related to phase boundary computation (same as in Fig. 8C). Seismic values extracted from (Master et al. 2000) and uncertainties (gray shaded areas) were summarized in (Karato and Karki 2001). From Wentzcovitch et al. 2006.

caused by the post-perovskite transition along our phase boundary (Fig. 21 in Wentzcovitch et al. 2010). Very large values for RS/P (~6) and negative values for Rφ/S (−0.5) result from this phase change. Therefore, lateral variation in phase abundances enhances the seismic parameters in the correct directions. The presence of secondary phases, such as ferropericlase and CaSiO3, should decrease these ratios. Besides, the real multicomponent system should produce a PT domain for co-existence between these phases, decreasing these anomalies. Nevertheless, lateral variation in the abundances of these phases produces robust results in the right direction. As pointed out from the outset, the topography of D′′ is consistent with a solid-solid phase transition with positive Clapeyron slope (4-10 MPa/K) induced by lateral temperature variations

Thermoelasticity of Mantle Minerals

117

(Sidorin et al. 1999). Interestingly, the correlation between geographic location and nature of Vφ and VS anomalies is also consistent with this assumption if post-perovskite is present in D′′. Beneath the Central Pacific Vφ is faster and VS is slower than their spherical mantle averages, whereas beneath the Circum-Pacific the opposite holds (see Masters et al. 2000 and references therein). These regions are generally considered to be respectively hotter and colder than average. The positive Clapeyron slope of this transition implies that hotter regions should contain less post-perovskite and therefore have faster Vφ and slower VS, whereas colder regions should be enriched in the latter and have slower Vφ and faster VS, as observed. Lateral variations in phase abundances alone do not account for the complex structure and properties of D′′. However, this perspective should be useful as a reference model from which interpretations of deviations might be sought. These results indicate that Rρ/S, the third parameter in Equation (20), caused by the postperovskite transition, should be positive and perhaps increase with depth in the lowermost mantle (see bottom row in Fig. 9), unless anelasticity or chemical heterogeneities occur simultaneously. Some 3D density models do not support this prediction (Ishii and Tromp 1999), whereas others do (Romanowicz 2001). Estimates of this parameter obtained by joint inversions of seismic and geodynamic data also tend to offer a positive ratio (Forte et al. 1994). It appears that until a consensus on the 3D structure of ρ is reached, this issue will remain open. It has been pointed out that an average excess density of 0.4% in the lowermost mantle is possible (Masters and Gubbins 2003), in agreement with our expectations. Another pressing issue is the anisotropy of D′′, a case of boundary layer anisotropy. Flow in these regions have a significant horizontal component and can align grains with preferred orientation if some of the phases present are in the dislocation creep regime. Despite multiple anisotropy styles observed (Wysession et al. 1998) and the numerous possible sources of anisotropy (Kendall and Silver 1998), evidences suggest that anisotropy in certain places of D′′ could result from lattice preferred orientation in largely strained aggregates produced be mantle circulation. Beneath regions of past and present subduction, such as the Caribbean, where colder than average and largely deformed paleo-plates are expected to reside and horizontal flow is expected, transverse anisotropy with VSH > VSV is generally observed (VSH(SV) is the velocity of horizontally (vertically) polarized shear waves propagating horizontally). It is plausible that in post-perovskite the primary slip system involves (010) (silica “layers” reside in this plane). Lateral material displacement could then align mainly the silica layers parallel to the horizontal plane. However, this simple picture did not provide a satisfactory explanation for the anisotropy in D′′ (Tsuchiya et al. 2004). Figure 10 shows “shear wave splittings” at conditions of the thermodynamic postperovskite phase boundary (see Figure 21 in Wentzcovitch et al. 2010) for some preferred alignments of the crystalline axes of MgO, perovskite, and post-perovskite in the vertical direction of transversely anisotropic aggregates. Such aggregates have a particular crystalline axis oriented vertically and random orientation of axes in the horizontal plane (see Fig. 10 caption). It is seen that at relevant conditions (i) the vertical alignment of [001] in postperovskite produces the largest positive (VSH − VSV) splitting; (ii) the shear wave splittings in perovskite and post-perovskite have similar magnitudes regardless of orientation; (iii) vertical alignment of perovskite’s [100] produces positive (VSH − VSV); and (iv) horizontal alignment of MgO’s {100}, its primary slip plane at high Ps (Yamazaki and Karato 2002), produces considerably larger (VSH − VSV). Despite lack of direct information on the slip systems of these phases at D′′ conditions, it appears that unlikely that post-perovskite can be a more significant source of anisotropy in D′′ than perovskite. Despite being less abundant (20 vol%) periclase is the most anisotropic phase. It is also weaker than the other phases. It is therefore likely to undergo more extensive deformation (Yamazaki and Karato 2002) and be a more important source of anisotropy in D′′.

118

Wentzcovitch, Wu, Carrier

Figure 10. Shear wave splittings, (VSH − VSV), in MgSiO3 post-perovskite and perovskite at the post-perovskite phase boundary conditions with respective uncertainty. VSH(SV) is the velocity of a horizontally (vertically) polarized shear wave propagating horizontally in aggregates with transverse anisotropy. Such aggregates have one crystalline axis aligned vertically and randomly oriented of axes in the horizontal plane. For [001] aligned vertically, VSH = (N/ρ)1/2 and VSV = (L/ρ) where, N = (1/8) (C11 + C22) – (1/4)C12 + (1/2) C66 and L = (1/2)(C44 + C55) (Wentzcovitch et al. 1998). Velocities for aggregates with [100] and [010] aligned vertically can be obtained by cyclic permutation of indices. From Wentzcovitch et al. (2006).

SELF-CONSISTENT QHA The QHA, as discussed here and in the previous article (Wentzcovitch et al. 2010), is a simple and powerful method for evaluating free energies using phonon frequencies obtained by first principles (Baroni et al. 2010 in this volume). There is still another useful consequence of this approximation not often appreciated: it can predict also complex crystal structures at high pressures and temperatures using only results obtained for equilibrium (eij = 0) configurations (Carrier et al. 2007). For clarity sake we repeat here the expression of the QHA free energy: F (V , T ) = E (V ) +

  ω ( V)   1 ωq (V ) + kBT ∑ ln  1 − exp  − q  ∑  2 q kBT   q  

(21)

U(V) is the DFT static energy versus volume obtained after full structural relaxation under isotropic pressure and ωq(V) is the phonon spectrum for these fully relaxed structures. Crystal structure parameters and phonon frequencies are uniquely related to volume since their values under pressure have been determined by static calculations and are given at Pstat(V) only. At T=T′, the pressure P′(V,T′), contains the following contributions: P′ = Pstat + PZP + Pth

(22)

where, Pstat, PZP, and Pth are the negative volume derivatives of the three terms in the righthand side of Equation (21), respectively. No further structural relaxation is performed at T′. The function V(P′,T′) is obtained by inverting P′(V,T′). Therefore, if V(P′,T′) = V(Pstat),

Thermoelasticity of Mantle Minerals

119

structural parameters at P′,T′ are equal to those at Pstat. This is the statically constrained QHA (sc-QHA). The validity of this approximation can be tested by comparing its predictions with some available high P,T experimental data. If the result is favorable, it asserts the validity of the sc-QHA in that P,T range. In general, this is expected if the temperature is in the domain of validity of the QHA. This can be assessed by a posteriori inspection of the predicted QHA thermal expansivity (see Fig. 8 in Wentzcovitch et al. 2010). MgSiO3-perovskite is one of the materials most studied at high pressures and temperature. There is a wealth of crystallographic data (Fiquet et al. 1998; Funamori et al. 1996; Ross and Hazen 1989, 1990; Utsumi et al. 1995) available up to 60 GPa and ~2,670 K. Its orthorhombic structure has 20 atoms per primitive cell and 10 degrees of freedom: 3 lattice parameters and 7 internal parameters. Figure 11 (see caption) displays experimental lattice parameters (data points) obtained at various conditions plotted against volume, irrespective of pressure and temperature, and compares them with parameters determined by static LDA calculations (full line) (Carrier et al. 2007). Clearly all experimental lattice parameters obtained at various

Figure 11. Comparison between lattice parameters predicted by the QHA and experimental data as a function of volume. The LDA and LDA + zero-point-motion-energy equilibrium volumes are also compared to the experimental equilibrium volume at 0 GPa. Experimental data are from Figuet et al. (1998), Funamori et al. (1996), Ross and Hazen (1989, 1990), Utsumi et al. (1995), Wang et al. (1994). Temperatures vary from 295 to 1024 K for Wang et al. (1994) from 293 to 2668 K for Fiquet et al. (1998), from 298 to 1173 K for Utsumi et al. (1995), from 293 to 2000 K for Funamori et al. (1996) and from 77 to 400 K for Ross and Hazen (1989, 1990). Solid lines are from static calculations. From Carrier et al. (2007).

120

Wentzcovitch, Wu, Carrier

P,T’s lie on or very close to the theoretical lines. Although VLDA < Vexp at 0 GPa, measured and computed structural parameters at the same volume agree very well. The effect of zero point motion energy on the calculated zero pressure volume is indicated as well. It increases volume significantly and should be taken into account in the performance analysis of exchangecorrelation energy functionals. Closer examination of Figure 11 at small volumes indicates a small but systematic deviation of experimental data from the theoretical lines. Some of these data points include very high temperature data. For instance, both data of Fiquet et al. (1998) and of Funamori et al. (1996) for the lattice parameter a are slightly larger (~0.3% or 0.01 Å) than the QHA prediction. The opposite is observed for b (~ −0.1% or −0.005 Å). No clear discrepancies are noticeable in the lattice parameter c. Since these points are the highest T data sets, one might suspect the validity of the QHA is questionable, but this is not the case. The origin of these discrepancies can be traced to deviatoric thermal stresses. Pstat in Equation (22) is isotropic by construction, since structures were optimized under hydrostatic conditions (Wentzcovitch et al. 1993). However, PZP and Pth are not necessarily isotropic. Deviatoric thermal stresses (in Voigt notation) are: δσTi (T , P ) = P +

1 ∂F V ∂ei

(23) f ,T

These stresses along [100], [010], and [001] are shown in Figure 12. Negative (positive) values indicate the structure should expand (contract) along the corresponding direction. One can see

Figure 12. Deviatoric stresses (Eqn. 23) along a, b, and c axes. From Carrier et al. (2007).

Thermoelasticity of Mantle Minerals

121

that there is a natural relation between the high T deviations seen in Figure 11 and the deviatoric stresses in Figure 12. These deviatoric stresses can be relaxed to first order if one knows the compliance tensor: κij (T , P ) = Cij−1 (T , P )

(24)

Since Cij are available for perovskite (Fig. 5), deviatoric stresses in Equation (23) can be relaxed by application of the following strains: εi ( P, T ) = κij ( P, T )δσ j ( P, T )

(25)

Figure 13 displays the resulting corrections on the lattice parameters as a function of volume at various temperatures, combined with the experimental data of Fiquet et al. (1998) and Funamori et al. (1996). The agreement between predicted and measured experimental data improves considerably despite fluctuations in the experimental data. The relaxed volume, V(1)(P,T), under “hydrostatic” conditions differs from Vstat(P). Therefore, phonon frequencies ωq(V(1)), might be recalculated and also become temperature dependent. This dependence is

Figure 13. Corrected lattice parameters [a−ε1a],[b−ε2b], and [c−ε3c] where εi’s are given by Equation (25), compared with the data of Funamori et al. (1996) data between 400 and 1500 K (dots) and Fiquet et al. (1998) data above 1500 K (triangles). Solid lines are from theory. As indicated by the Fiquet et al. (1989) data, experimental errors increase with temperature. From Carrier et al. (2007).

122

Wentzcovitch, Wu, Carrier

not caused by the usual phonon-phonon interaction but by the non-isotropic nature of the thermal pressure. This procedure should be repeated until deviatoric thermal stresses vanish, at which point all thermally induced forces should also vanish. This is the self-consistent QHA (sc-QHA). Corrections of elastic constants to first order are also possible and have been carried out for perovskite and post-perovskite (Carrier et al. 2008) at D′′ conditions (P ~ 120 GPa and T ~ 3,000 K). For these orthorhombic materials, the diagonal and off diagonal components of cij expanded in a Taylor series of strains (in Voigt’s notation) are: 6

∂Cij

k =1

∂ε k

Cij ( P, T , ε) = Cij ( P, T , 0) + ∑

ε k + ...

(26)

P ,T

or ∂Cij k =1 m =1 ∂P′ 6

6

Cij ( P, T , ε) = Cij ( P, T , 0) + ∑ ∑

P ,T

∂P′ ∂σ m

P ,T

∂σ m ∂ε k

ε k ...

(27)

P ,T

where P′ =

1 3 ∑σm 3 m =1

(28)

Therefore the first-order corrected elastic constants at the strains given by Equation (25) are: Cij ( P, T , ε) = Cij ( P, T , 0) +

∂Cij ∂P′

P ,T

δσΣ + ... 3

(29)

where δσΣ is the sum of deviatoric stresses: 3

δσΣ = ∑ δσ m

(30)

m =1

each given by: 6

δσ m = ∑ Cmk ε k

(31)

k =1

This correction requires knowledge of dCij/dP and of the deviatoric stresses, both of which are known from the sc-QHA calculation. Changes in the isothermal bulk and shear moduli of perovskite and post-perovskite after one cycle of thermal stress relaxation are displayed in Figure 14. In perovskite, these changes are within 0.5 GPa up to 150 GPa and 4,000 K. This is within the margin of uncertainty of the calculations. In post-perovskite the changes are a little larger and negative. The changes are largest at the highest P,T but still relatively very small compared with values shown in Figs. 6a and 8a.

SUMMARY We have presented the formalism used to compute thermoelastic properties of solids using the statically constrained QHA (sc-QHA). Combined with first principles LDA calculations it is a simple and accurate although computationally intensive method to compute thermoelastic properties of crystals from which elasticity of aggregates can be obtained within bounds. Results for MgO, MgSiO3-perovskite, and MgSiO3-post-perovskite and comparisons with high

Thermoelasticity of Mantle Minerals

123

Figure 14. Corrections to bulk and shear moduli of (a) pervskite (Wentzcovitch et al. 2004) and (b) postperovskite (Wentzcovitch et al. 2006). From Carrier et al. (2008).

pressure and high temperature experimental data have been reviewed to illustrate the accuracy possible to achieve. An interesting byproduct of QHA calculations is the prediction of high temperature crystal structures. This is possible because, by construction, the QHA relates crystal structure parameters and (non-interacting) phonon frequencies uniquely with volume and such structure-volume relationship is well established by static calculations. In non-cubic solids, sc-QHA calculations develop deviatoric thermal stresses at high temperatures. Thermal pressure is not isotropic. Relaxation of these stresses leads to a series of corrections to the structure that may be taken to any desired order, up to self-consistency. This is the self-consistent QHA (sc-QHA). We have shown how to correct elastic constants for deviatoric stresses generated by the sc-QHA. We have illus-

124

Wentzcovitch, Wu, Carrier

trated the procedure by correcting to first order the elastic constants of perovskite and post-perovskite, the major silicate phases of the Earth’s lower mantle. This correction is very satisfactory for obtaining the aggregate elastic constants and velocities of these minerals at in situ conditions of the lower mantle. This procedure can also be used to predict elastic constants in the presence of deviatoric stresses, or to correct elasticity measurements performed under non-hydrostatic conditions, as often happen in diamond-anvil cells, if deviatoric stresses are known. Several insights of geophysical significance have been obtained from high temperature elasticity calculations: 1) MgO is the most anisotropic phase in the lower mantle. Since it is also the weakest, it is a potentially important source of anisotropy in aggregates with lattice preferred orientation produced by mantle flow; 2) aggregates with typical pyrolite Mg/Si ratio, i.e., ~20 vol% of MgO and ~80 vol% of perovskite, appear to reproduce the elastic properties of the lower mantle better than the pure perovskite aggregate, at least down to ~1,600 km depth. This suggests that shallow to mid lower mantle has the same chemistry of the upper mantle (as far as Mg/Si ratio is concerned); 3) in the deep lower mantle, post-perovskite free aggregates appear to slowly develop elastic properties that deviate from those of the Preliminary Reference Earth Model (PREM). The elastic properties of post-perovskite suggest that this deviation is caused by its presence in the D′′ layer, the deepest ~300 km of the mantle; 4) the post-perovskite transition causes velocity jumps similar to those detected in some places in D′′ beneath subduction zones. This is consistent with expectations based on the post-perovskite phase boundary with positive Clapeyron slope; 5) velocity changes across the post-perovskite transition suggest that the anticorrelation between lateral bulk and shear velocity changes in D′′ could be caused by lateral changes in the perovskite/post-perovskite abundances. For almost two decades the D′′ region has remained an enigma. It is a complex layer at the interface of two chemically distinct regions, mantle and core. Post-perovskite should co-exist with other solid phases, and probably also with melts, in this region. There are still several mysteries, in the details, to be resolved; but, the insights offered by these calculations have advanced considerably knowledge of the puzzling D′′ layer. These successes show that first principles calculations of thermoelastic properties of mantle minerals are poised to contribute much more to our understanding of the planet in the years ahead.

ACKNOWLEDGMENTS This research was supported by NSF grants EAR 0810272, EAT 0635990, ATM 0428774 (VLab), and EAR 0757903.

REFERENCES Anderson GM (2005) Thermodynamics of Natural Systems. 2nd edition. Cambridge Univ. Press, Cambridge, UK Baroni S, Giannozzi P, Isaev E (2010) Density-functional perturbation theory for quasi-harmonic calculations. Rev Mineral Geochem 71:39-57 Baroni S, Dal Corso A, Giannozzi P, de Gironcoli S (2001) Phonons and related crystal properties from densityfunctional perturbation theory. Rev Mod Phys 73:515-562 Barron THK, Klein ML (1965) Second-order elastic constants of a solid under stress. Proc Phys Soc 85:523532 Born M, Huang K (1956) Dynamical Theory of Crystal Lattices. International Series of Monographs on Physics. Oxford at the Clarendon Press, Hong Kong Brown JM, Shankland TJ (1981) Thermodynamic parameters in the Earth as determined from seismic profiles. Geophys J R Astron Soc 66:579-596 Car R, Parrinello M (1985) Unified approach for molecular dynamics and density functional theory. Phys Rev Lett 55:2471-2474

Thermoelasticity of Mantle Minerals

125

Caracas R, Cohen R (2007) The effect of chemistry on the physical properties of perovskite and post-perovskite. In: Post-perovskite. The Last Mantle Phase Transition. Hirose K, Brodholt J, Lay T, Yuen D (eds) AGU Monograph Series 174:115-128 Caracas R, Wentzcovitch RM, Price GD, Brodholt J (2005) CaSiO3 perovskite at lower mantle pressures. Geohys Res Lett 144:L06306 Carrier P, Justo JF, Wentzcovitch RM (2008) Quasiharmonic elastic constants corrected for deviatoric thermal stresses. Phys Rev B 78:144302 Carrier P, Wentzcovitch RM, Tsuchiya J (2007) First principles prediction of crystal structures at high temperatures using the quasiharmonic approximation. Phys Rev B 76:064116 Ceperley DM, Alder BJ (1980) Ground state of the electron gas by a stochastic method. Phys Rev Lett 45:566569 Chen G, Liebermann RC, Weidner DJ (1998) Elasticity of single crystal MgO to 8 gigapascals and 1600 Kelvin. Science 280:1913-1916 da Silva CRS, Wentzcovitch RM, Patel A, Price GD, Karato S-I (2000) The composition and geotherm of the lower mantle: constraints from the calculated elasticity of silicate perovskite. Phys Earth Planet Inter 118:103-109 Davies GF (1974) Effective elastic-moduli under hydrostatic stress. 1. Quasiharmonic theory. J Phys Chem Solids 35:1513-1520 Duffy TS, Hemley RJ, Mao HK (1995) Equation of state and shear strength at multi-Mbar pressures: magnesiumoxide to 227 GPa. Phys Rev Lett 70:1371-1374 Dziewonski AM, Anderson DL (1981) Preliminary reference Earth model. Phys Earth Planet Inter 25: 297-356 Dziewonski AM, Woodhouse JH (1987) Global images of the Earth interior. Science 236:37-48 Fiquet G, Andrault D, Dewaele A, Charpin T, Kunz M, Haüsermann D (1998) P-V-T equation of state of MgSiO3 perovskite. Phys Earth Planet Inter 105:21-31 Forte AM, Woodward RL, Dziewonski AM (1994) Joint inversions of seismic and geodynamic data for models of 3 dimensional mantle heterogeneity. J Geophys Res 99:21857-21878 Funamori N, Yagi T, Utsumi W, Kondo T, Uchida T, Funamori M (1996) Thermoelastic properties of MgSiO3 perovskite determined by in situ x-ray observations up to 30 GPa and 2000 K. J Geophys Res 101:82578269 Garnero E (2004) A new paradigm for Earth’s core-mantle boundary. Science 304:834-836 Giannozzi P, Baroni S, Bonini N, Calandra M, Car R, Cavazzoni C, Ceresoli D, Chiarotti GL, Cococcioni M, Dabo I, Dal Corso A, de Gironcoli S, Fabris S, Fratesi G, Gebauer R, Gerstmann U, Gougoussis C, Kokalj A, Lazzeri M, Martin-Samos L, Marzari N, Mauri F, Mazzarello R, Paolini S, Pasquarello A, Paulatto L, Sbraccia C, Scandolo S, Sclauzero G, Seitsonen AP, Smogunov A, Umari P, Wentzcovitch RM (2009) QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials. J Phys Condens Matter 21:395502 Giannozzi P, de Gironcoli S, Pavone P, Baroni S (1991) Ab initio phonon dispersions in elemental semiconductors Phys Rev Lett 43:7231-7234 Gonze X, Beuken J-M, Caracas R, Detraux F, Fuchs M, Rignanese G-M, Sindic L, Verstraete M, Zerah G, Jollet F, Torrent M, Roy A, Mikami M, Ghosez Ph, Raty J-Y, Allan DC (2002) First-principles computation of material properties : the ABINIT software project. Comput Mater Sci 25:478-492 Grand SP (1994) Mantle shear structure beneath the America and surrounding oceans. J Geophys Res 99:1159111621 Hashin Z, Strikman S (1962) A variational approach to the theory of the elastic behavior of polycrystals. J Mech Phys Solids 10:343-352 Hohenberg P, Kohn W (1964) Inhomogeneous electron gas. Phys Rev 136:B864-B871 Hsu H, Umemoto K, Wu A, Wentzcovitch RM (2010) Spin-state crossover of iron in lower-mantle minerals: results of DFT+U investigations. Rev Mineral Geochem 71:169-199 Isaak DG, Andersen OL, Goto T (1989) Measured elastic moduli of single-crystal MgO up to 1800 K. Phys Chem Miner 16:704-713 Isaak DG, Cohen RE, Mehl MJ (1990) Calculated elastic and thermal properties of MgO at high pressures and temperatures. J Geophys Res 95:7055-7067 Ishii M, Tromp J (1999) Normal-mode and free-air gravity constraints on lateral variations in velocity and density of Earth’s mantle. Science 285:1231-1236 Jackson I (1998) Elasticity, composition and temperature of the Earth’s lower mantle: a reappraisal. Geophys J Int 134:291-311 Jackson JM, Zhang J, Shu J, Sinogeikin SV, Bass JD (2005) High-pressure sound velocities and elasticity of aluminous MgSiO3 perovskite to 45 GPa: Implications for lateral heterogeneity in Earth’s lower mantle. Geophys Rev Lett 32:L21305 Karato S-I, Karki BB (2001) Origin of lateral variation of seismic wave velocities and density in the deep mantle. J Geophys Res 106:B21771-B21783

126

Wentzcovitch, Wu, Carrier

Karki BB, Stixrude L, Clark SJ, Warren MC, Ackland GJ, Crain J (1997) Structure and elasticity of MgO at high pressure. Am Mineral 82:51-60 Karki BB, Stixrude L, Wentzcovitch RM (2001a) High-pressure elastic properties of major materials of earth’s mantle from first principles. Rev Geophys 39: 507-534 Karki BB, Wentzcovitch RM, de Gironcoli S, Baroni S (1999) First-principles determination of elastic anisotropy and wave velocities of MgO at lower mantle conditions. Science 286:1705-1707 Karki BB, Wentzcovitch RM, de Gironcoli S, Baroni S (2000a) Ab initio lattice dynamics of MgSiO3‑perovskite. Phys Rev B 62:14750-14756 Karki BB, Wentzcovitch RM, de Gironcoli S, Baroni S (2000b) High-pressure lattice dynamics and thermoelasticity of MgO. Phys Rev B 61:8793-8800 Karki BB, Wentzcovitch RM, de Gironcoli S, Baroni S (2001b) First principles thermoelasticity of MgSiO3‑perovskite: consequences for the inferred properties of the lower mantle. Geophys Res Lett 28:2699-2702 Kendall J-M, Silver PG (1998) Investigating causes of D′′ anisotropy. In: Core-Mantle Boundary Region. Geodynamics Series. Volume 28. Gurnis M, Wysession ME, Knittle E, Buffet B (eds) Am Geophys Union, Washington, DC, p 97-118 Kiefer B, Stixrude L, Wentzcovitch RM (2002) Elasticity of (Mg,Fe)SiO3-Perovskite at high pressures. Geophys Res Lett 29:14683 Kohn W, Sham LJ (1965) Self-consistent equations including exchange and correlation effects. Phys Rev 140:1133-1138 Kresse G, Furthmuller J (1996) Efficient iterative schemes for ab initio total-energy calculations using a planewave basis set. Phys Rev B 54:11169-11186 Lay T, Helmberger DV (1983) A lower mantle S-wave triplication and the shear velocity structure of D′′. Geophys J R Astron Soc 75:799-838 Li B, Zhang J (2005) Pressure and temperature dependence of elastic wave velocity of MgSiO3 perovskite and the composition of lower mantle. Phys Earth Planet Inter 151:143-154 Li BS, Woody K, Kung J (2006a) Elasticity of MgO to 11 GPa with an independent absolute pressure scale: Implication for pressure calibration. J Geophys Res 111:B11206 Li L, Brodholt J, Stackhouse S, Weidner DJ, Alfredsson M, Price GD (2005) Elasticity of (Mg, Fe)(Si, Al)O3perovskite at high pressure. Earth Planet Sci Lett 240:529-536 Li L, Weidner DJ, Brodholt J, Alfè D, Price GD (2006b) Elasticity of Mg2SiO4 ringwoodite at mantle conditions. Phys Earth Planet Inter 157:181-187 Li L, Weidner DJ, Brodholt J, Alfè D, Price GD (2009) Ab initio molecular dynamics study of elasticity of akimotoite MgSiO3 at mantle conditions. Phys Earth Planet Inter 173:115-120 Li L, Weidner DJ, Brodholt J, Alfè D, Price GD, Caracas R, Wentzcovitch R (2006c) Elasticity of CaSiO3 perovskite at high pressure and high temperature. Phys Earth Planet Inter 155:249-259 Marton FC, Cohen RE (2002) Constraints on lower mantle composition from molecular dynamics simulations of MgSiO3 perovskite. Phys Earth Planet Inter 134:239-252 Masters G, Gubbins D (2003) On the resolution of density within the Earth. Phys Earth Planet Inter 140:159167 Masters G, Laske G, Bolton H, Dziewonski A (2000) The relative behavior of shear velocity, bulk sound speed, and compressional velocity in the mantle: implications for chemical and thermal structure. In: Earth’s Deep Interior: From Mineral Physics and Tomography from Atomic to the Global Scale. Karato S-I, Forte AM, Liebermann RC, Masters G, Stixrude L (eds) American Geophysical Union, Washington DC. Geophysical Monograph Series 117:63-87 Murakami M, Hirose K, Kawamura K, Sata N, Ohishi Y (2004) Post-perovskite phase transition in MgSiO3. Science 304:855-858 Murakami M, Ohishi Y, Hirao N, Hirose K (2009) Elasticity of MgO to 130 GPa: Implications for lower mantle mineralogy. Earth Planet Sci Lett 277:123-129 Musgrave MJP (1970) Crystal Acoustics. Holden-Day, Boca Raton, FL Nielsen OH, Martin RM (1985) Quantum mechanical theory of stress and force. Phys Rev B 32:5780-3891 Nye JF (1985) Physical Properties of Crystals: Their Representation by Tensors and Matrices, 2nd ed. Oxford University Press, Oxford Oganov AR, Brodholt JP, Price GD (2001) The elastic constant of MgSiO3 perovskite at pressure and temperatures of the Earth’s mantle. Nature 411:934-937 Oganov AR, Ono S (2004) Theoretical and experimental evidence for a post-perovskite phase of MgSiO3 in Earth’s D′′ layer. Nature 430: 445-448 Payne MC, Teter MP, Allan DC, Arias TA, Joannopoulos JD (1992) Iterative minimization techniques for ab initio total-energy calculations: molecular dynamics and conjugate gradients. Rev Mod Phys 64:10451097

Thermoelasticity of Mantle Minerals

127

Perdew J, Zunger A (1981) Self-interaction correction to density-functional approximations for many electron systems. Phys Rev B 23: 5048–5079 Perdew JP, Burke K, Ernzerhof M (1996) Generalized gradient approximation made simple. Phys Rev Lett 77:3865-3868 Romanowicz B (2001) Can we resolve 3D density heterogeneity in the lower mantle? Geophys Res Lett 28:1107-1110 Ross N, Hazen R (1989) Single crystal X-ray diffraction study of MgSiO3 perovskite from 77 to 400 K. Phys Chem Miner 16:415-420 Ross N, Hazen R (1990) High-pressure crystal chemistry of MgSiO3 perovskite. Phys Chem Miner 17:228237 Sidorin I, Gurnis M, Helmberger D (1999) Evidence for a ubiquitous seismic discontinuity at the base of the mantle. Science 286:1326-1331 Sinelnikov YD, Chen G, Neuville DR, Vaughan MT, Liebermann RC (1998) Ultrasonic shear wave velocities of MgSiO3 perovskite at 8 GPa and 800 K and lower mantle composition. Science 281:677-679 Sinogeikin SV, Bass JD (1999) Single-crystal elasticity of MgO at high pressure. Phys Rev B 59:R14141-R14144 Sinogeikin SV, Bass JD, Katsura T (1998) Single crystal elasticity of MgO at high pressure. Phys Rev B 59: R14141-R14144 Sinogeikin SV, Zhang J, Bass JD (2004) Elasticity of single crystal and polycrystalline MgSiO3 perovskite by Brillouin spectroscopy. Geophys Rev Lett 31:L06620 Spetzler HA (1970) Equation of state of polycrystalline and single-crystal MgO to 8 kbars and 800 K. J Geophys Res 75:2073-2087 Stackhouse S, Brodholt JP, Wookey J, Kendall JM, Price GD (2005) The effect of temperature on acoustic anisotropy of the perovskite and post-perovskite polymorphs of MgSiO3. Earth Planet Sci Lett 230:1-10 Stixrude L, Cohen RE, Yu RC, Krakauer H (1996) Prediction of phase transition in CaSiO3 perovskite and implications for lower mantle structure. Am Mineral 81:1293-1296 Stixrude L, Hemley RJ, Fei Y, Mao HK (1992) Thermoelasticity of silicate perovskite and magnesiowustite and stratification of the Earth mantle. Science 257:1099-1101 Stixrude L, Lithgow-Bertelloni C (2005) Thermodynamics of mantle minerals. I. Physical properties. Geophys J Int 162:610-632 Stixrude L, Lithgow-Bertelloni C (2007) Influence of phase transformations on lateral heterogeneity and dynamics in Earth’s mantle. Earth Planet Sci Lett 263:45-55 Stixrude L, Lithgow-Bertelloni C (2010) Thermodynamics of the Earth’s mantle. Rev Mineral Geochem 71:465484 Stixrude L, Lithgow-Bertelloni C, Kiefer B, Fumagalli P (2007) Phase stability and shear softening in CaSiO3 perovskite at high pressure, Phys Rev B 75:24108 Thomsen L (1972) Fourth-order anharmonic theory—elasticity and stability. J Phys Chem Solids 33:363-378 Troullier N, Martins JL (1991) Efficient pseudopotentials for plane-wave calculations. Phys Rev B 43:19932006 Tsuchiya J, Tsuchiya T, Wentzcovitch RM (2005) Vibrational and thermodynamic properties of MgSiO3 postperovskite. J Geophys Res 110:B02204, doi:10.1029/2004JB003409 Tsuchiya T, Tsuchiya J, Umemoto K,Wentzcovitch RM (2004) Phase transition in MgSiO3 perovskite in the earth’s lower mantle. Earth Planet Sci Lett 224:241–248 Utsumi W, Funamori N, Yagi T, Ito E, Kikegawa T, Shimomura O (1995) Thermal expansivity of MgSiO3 perovskite under high pressures up to 20 GPa. Geophys Res Lett 22:1005-1008 van der Hilst RD, Widiyantoro S, Engdahl ER (1997) Evidence for deep mantle circulation from global tomography. Nature 386:578-584 Vanderbilt D (1990) Soft self-consistent pseudopotentials in a generalized eigenvalue formalism. Phys Rev B 41:7892-7895 Wallace DC (1972) Thermodynamics of Crystals. 1st Edition. John Wiley and Sons, New York Wang Y, Perdew J (1991) Correlation hole of the spin-polarized electron-gas, with exact small-wave-vector and high-density scaling. Phys Rev B 24:13298-13307 Wang Y, Weidner DJ, Liebermann RC, Zhao Y (1994) P-V-T equation of state of (Mg,Fe)SiO3 perovskite: constraints on composition of the lower mantle. Phys Earth Planet Inter 83:13-40 Watt JP (1980) Hashin-Shtrikman bounds on the effective elastic moduli of polycrystals with monoclinic symmetry. J Appl Phys 51:1520-1524 Watt JP, Davies GF, Connell RJO (1976) Elastic properties of composite materials. Rev Geophys Space Phys 14:541-563 Wentzcovitch RM (1991) Invariant molecular dynamics approach to structural phase transitions. Phys Rev B 44:2358-2361

128

Wentzcovitch, Wu, Carrier

Wentzcovitch RM, Justo JF, Wu Z, da Silva CRS, Yuen D, Kohlstedt D (2009) Anomalous compressibility of ferropericlase throughout the iron spin crossover. Proc Nat Acad Sci USA 106:8447-8452 Wentzcovitch RM, Karki BB, Cococcioni M, de Gironcoli S (2004) Thermoelastic properties of MgSiO3perovskite: insights on the nature of the Earth’s lower mantle. Phys Rev Lett 92:018501 Wentzcovitch RM, Karki BB, Karato SI, da Silva CRS (1998) High pressure elastic anisotropy of MgSiO3perovskite and geophysical implications. Earth Planet Sci Lett 164:371-378 Wentzcovitch RM, Martins JL (1991) First principles molecular dynamics of Li: Test of a new algorithm. Solid State Commun 78:831-834 Wentzcovitch RM, Martins JL, Price GD (1993) Ab initio molecular dynamics with variable cell shape: Application to MgSiO3 perovskite. Phys Rev Lett 70:3947-3950 Wentzcovitch RM, Ross N, Price GD (1995) Ab initio investigation of MgSiO3 and CaSiO3-perovskites at lower mantle pressures. Phys Earth Planet Inter 90:101-112 Wentzcovitch RM, Tsuchiya T, Tsuchiya J (2006) MgSiO3 post perovskite at D″ conditions. Proc Nat Acad Sci 103:543-546 Wentzcovitch RM, Yu YG, Wu Z (2010) Thermodynamic properties and phase relations in mantle minerals investigated by first principles quasiharmonic theory. Rev Mineral Geochem 71:59-98 Woodhouse JH, Dziewonski AM (1984) Mapping the upper mantle: three-dimensional modeling of earth structure by inversion of seismic waveforms, J Geophys Res 89:5953-5986 Wu Z, Justo JF, da Silva CRS, de Gironcoli S, Wentzcovitch RM (2009) Anomalous thermodynamic properties in ferropericlase throughout its spin crossover transition. Phys Rev B 80:014409 Wu Z, Wentzcovitch RM (2009) Effective semiempirical ansatz for computing anharmonic free energies. Phys Rev B 79:104304 Wu Z, Wentzcovitch RM, Umemoto K, Li B, Hirose K (2008) P-V-T relations in MgO: an ultrahigh P-T scale for planetary sciences applications. J Geophys Res 113:B06204 Wysession ME, Lay T, Revenaugh J, Williams Q, Garnero EJ, Jeanloz R, Kellogg LH (1998) The D′′ discontinuity and its implications. In: Core-Mantle Boundary Region. Geodynamics Series. Volume 28. Gurnis M, Wysession ME, Knittle E, Buffet B (eds) Am Geophys Union, Washington, DC, p 273-297 Xu WB, Lithgow-Bertelloni C, Stixrude L, Ritsema J (2008) The effect of bulk composition and temperature on mantle seismic structure. Earth Planet Sci Lett 275:70-79 Yamazaki D, Karato S (2002) Fabric development in (Mg,Fe)O during large strain, shear deformation: implications for seismic anisotropy in Earth’s lower mantle. Earth Planet Sci Lett 131:251-267 Yeganeh-Haeri A (1994) Synthesis and re-investigation of the elastic properties of single-crystal magnesium silicate perovskite. Phys Earth Planet Inter 87:111-121 Zha CS, Mao HK, Hemley R (2000) Elasticity of MgO and a primary pressure scale to 55 GPa. Proc Nat Acad Sci USA 97:13494-13499