Computational Mineral Physics and the Physical Properties of Mantle Minerals

1. Computational Mineral Physics and the Physical Properties of Mantle Minerals John Brodholt UCL

2. Outline of the talk • Different computational techniques: Molecular dynamics Ab initio (quantum mechanics) calculations Ab initio molecular dynamics. • Ab initio calculations on elastic constants of perovskite at lower mantle P and T and using this to make an estimate of temperature variations on the lower mantle. • The effect of aluminium on the physical properties of perovskite. • Future directions?

3. Why use computational mineral physics? • P and T at centre of Earth ~ 3.5 Mbar and ~ 6000 K • P and T at core/mantle ~ 1.3 Mbar and ~ 3000 to 4000 K • Piston cylinder ~ 40 kbar • Multi-anvil ~ 400 kbar • Diamond-anvil ~ 2 Mbar (temperatures uncertain, gradients high) • Shock guns ~ 2 Mbar (temperatures extremely uncertain) • High P and T experiments are hard and can have large uncertainties. Use computational minerals physics to predict the properties of Earth materials at mantle and core conditions. Also to understand properties at an atomistic level.

4. Predict: • Volumes, bulk modulus • Elastic constants (seismic velocities) • Heat capacities • Free energies of reactions (phase diagrams) • Defects • Diffusion • Viscosities • Element partitioning • etc

5. Molecular Dynamics • Simulate the properties of minerals, fluids, melts, etc at an atomistic (or sub-atomic) level. • If you know the forces (or potential energy) between the constituent atoms, you can find the total force on the atoms (or ions) and solve the equations of motions and move them through time • Obtain averages over time that give you the properties of the material of interest. • i.e. the kinetic energy gives the temperature.

6. In classical simulations, the potentials are parameterised empirically; fit to some experimental data.

8. Simulation of Na diffusing in Quartz

9. Matsui, P.E.P.I., 2000: (Uses breathing shell for the oxygen)

12. Pros and Cons of Empirical Potentials Pros: • Large systems (millions of atoms) • Long runs (millions of time steps ~ ns) Cons • Very difficult to assess how transferable the potentials are: will work in one phase but not another. So, solution is to use ab initio quantum mechanical techniques to get the forces.

13. HY = EY Alternatively, “I think I can safely say nobody understands quantum mechanics” Richard Feynman

14. Ab Initio Techniques • Numerically solving Schrodinger’s equation • One major approximation; the effect on any one electron from all the other electrons (a very serious many body problem) is wrapped up into a term called the exchange-correlation. • The exchange correlation functional can be either the Local Density Approximation (LDA) or the Generalised Gradient Approximation (GGA). • Also use pseudopotentials for the interactions between valence electrons and the tightly bound core electrons (not to be confused with empirical pair potentials). • Run on national supercomputers – T3E. • Until recently has been used only at 0K.

15. LDA vs GGA • LDA tends to underestimate volume (1-3%) and overestimate the bulk modulus (by a few %). • GGA tends to do the opposite – overestimate the volume and underestimate the bulk modulus. • GGA is better for reaction energies. • GGA is better for hydrogen bonded systems. This is just a fact of life.

20. Masters et al., 2000

21. 2740 km depth ~ D’’ Forte and Mitrovica, Nature, 26 April 2001

23. Oganov, Brodholt and Price, EPSL, 2001

25. ATHERMAL (0 K) ELASTIC PROPERTIES OF PEROVSKITE Volume is constrained to experimental zero pressure V. Calc. Exp. K(GPa) 267 264 G 179 177 C11(GPa) 493 482 C22 546 537 C33 470 485 C12 142 144 C13 146 147 C23 160 146 C44 212 204 C55 186 186 C66 149 147 According to Artem Oganov, “These are simply brilliant”

26. Ab Initio Molecular Dynamics • VASP (Vienna Ab-Initio Simulation Package) • Generalised Gradient Approximation (GGA) • Ultrasoft pseudopotential for O • Norm-conserving pseudopotential with partial core corrections for Mg and Si • Plane-wave cutoff of 500 eV • 1 kpt (gamma point) • 80 atom cell • NVT • Run lengths of about 2 ps with a 1 fs timestep • Run on CSAR T3E (which NERC has ¼).

29. Elastic constants of MgSiO3Perovskite • Oganov, Brodholt & Price, Nature 2001 – QMD on perovskite at P and T. • Remains orthorhombic in mantle P-T space. • Elastic constants obtained from finite deformation of box and stress/strain relations:

30. At 38 GPa dlnVs/dlnVp = 1.48 At 88 GPa dlnVs/dlnVp = 1.91

31. High Temperature Vp and Vs from Ab Initio Molecular Dynamics • R = 1.48 to 1.91 • If we extrapolate, R ~ 2.3 at core-mantle boundary • This is in accordance with other minerals. • Observed R is about 1.7 to perhaps greater that 3. • So, this supports the suggestion that the deeper lower mantle is either chemically heterogeneous or it is caused by anelasticity. • Anelasticity could perhaps take R up to about 2.7 (Karato, 2001). • Need to know the dependence of elastic constants on chemistry.

32. High Temperature Vp and Vs from Ab Inito Molecular Dynamics • VF = Sqrt(K/r) is relatively independent of composition and anelasticity. • Therefore, use VF to make an estimate of temperature variations in the lower mantle. • Masters et al, (2000) give a maximum δVF /VF of about 1.4% at 1225 and 2195 km. • Our results give: dlnVF/dT = -1.7×10-5 K -1 at 1000 km dlnVF/dT = -0.9×10-5 K -1 at 2000 km • Maximum temperature contrast is ~800 K at 1000 km and about ~1500 K at 2000 km. • Possibly > 2000 K at core-mantle boundary • Using Kennett et al, (1998) we get smaller temperature variations: 500 and 1000K.

33. 2740 km depth ~ top of D’’ Forte and Mitrovica, Nature, 26 April 2001

35. Summary • Experiments on MgSiO3 perovskites with only 5% Al2O3 have unexpectedly high compressibilities (mantle has ~ 4 wt%). • This is only consistent with a mechanism where Al3+ substitutes for Si4+, charge balanced by oxygen vacancies. • But, this mechanism is only favoured in the shallow part of the lower mantle. In the deeper part Al is incorporated be a vacancy free mechanism. • This means that alumina will only effect the properties of perovskite (compressibility, creep strength, diffusion, etc) in the shallow lower mantle; the deeper lower mantle will have properties in accord with alumina-free perovskite.

37. Incorporation of Al3+ in Perovskite from Empirical Potential Calculations • At 0 GPa Al3+ goes into Si site only, charge balanced by oxygen vacancies. • But P > 15 GPa, Al substitutes for both Si and Mg with no vacancies. Richmond and Brodholt (1999)

42. Homologous series of compounds along the oxygen deficient axis. Smyth et al. 1996 A2B2O5 A3B3O8 A4B4O11 AnBnO3n-1 ABO3

43. Computational Details • Used commercially available DFT code CASTEP. • Ultrasoft pseudopotentials • Default planewave cutoff (380 eV). Energy difference between 380 eV and 900 eV was 0.04 eV. • K-point grid 0.1 1/Å. Increased sampling gave negligible difference on equation-of-state and reaction enthalpies. • All atomic positions were fully relaxed.

45. (Mg3Al)(Si3Al)O12 • K = 224 GPa • V = 163.31 Ang**3 Adjust for 5% Al2O3and assume linear mixing • V = <1.001*V(MgSiO3) • Vexp = 1.005*V(MgSiO3) • K = < 1% smaller • Kexp = 10% smaller

46. Brownmillerite • K = 155 GPa • V = 345.16 Ang**3 Adjust for 5% Al2O3 and assume linear mixing • V = 1.006*V(MgSiO3) • Vexp = 1.005*V(MgSiO3) • K = 6% smaller • Kexp = 10% smaller

47. Oxygen Vacancies • Enhanced diffusion • Weaker – faster creep rates • Oxygen ionic conductor – implications for electrical conductivity • Seismic velocities might be lower • Source for water: Vo + O2- + H2O = 2OH-

49. Oxygen Vacancies • Might be an explanation for the mid-mantle discontinuity that is sometimes reported (900 to 1100 km) • Might be a zone of mantle below the lower mantle which is weaker than that below.

50. Cizkova et al., GRL 1999