Structure and Freezing of MgSiO3 Liquid in Earth's Lower Mantle

See allHide authors and affiliations

Science  14 Oct 2005:
Vol. 310, Issue 5746, pp. 297-299
DOI: 10.1126/science.1116952


First-principles molecular-dynamics simulations show that over the pressure regime of Earth's mantle the mean silicon-oxygen coordination number of magnesium metasilicate liquid changes nearly linearly from 4 to 6. The density contrast between liquid and crystal decreases by a factor of nearly 5 over the mantle pressure regime and is 4% at the core-mantle boundary. The ab initio melting curve, obtained by integration of the Clausius-Clapeyron equation, yields a melting temperature at the core-mantle boundary of 5400 ± 600 kelvins.

The structure of silicate liquids over most of the mantle pressure regime is unknown. There is experimental evidence that silicate liquid structure changes substantially on compression and that these changes include an increase in Si-O coordination number from 4 toward 6 (1), although there are no experimental observations within the liquid stability field at lower mantle pressure.

Most simulations of silicate liquids have been based on atomistic models, which permit much faster computation but have the disadvantage of being based on empirical force fields, the forms of which are uncertain (2, 3). First-principles molecular-dynamics (FPMD) simulations are considerably more costly because the electronic structure must be computed at each time step, but they are also more robust because they make no assumptions about the type of bonding or the shape of the charge density.

Here, we present FPMD simulations of MgSiO3 liquid over the pressure-temperature range of Earth's mantle, based on density functional theory in the local density approximation and the plane-wave pseudopotential method (4). The simulations are performed in the canonical ensemble with periodic boundary conditions and a Nosé thermostat (5). The primary cell is cubic, with 80 atoms (16 formula units). The initial condition is a pyroxene structure homogeneously strained to a cubic cell shape and the desired volume. The simulated strained crystal is melted at 6000 K and then cooled isochorically to 4000 K and 3000 K. Melting was confirmed by inspection of the radial distribution function and the mean square displacement. Total run durations are 3 ps, with the last 2.4 ps used for computing equilibrium averages. To ensure that the results are not sensitive to initial conditions, we performed simulations with twice the number of atoms (160 atoms initiated with a cubic majorite structure), twice the duration (6 ps), and different initial configurations (80-atom perovskite super-cell homogeneously strained to a cubic cell shape) and found that the computed properties of the liquid were unchanged within statistical uncertainties. We explore a range of volumes V = VX to V = VX/2, where VX = 38.9 cm3 mol–1 is the experimental value for the liquid at the ambient-pressure melting point (1830 K) (6).

To explore freezing, we also performed simulations of orthorhombic MgSiO3 perovskite, following the same procedures as the liquid simulations except that at each volume the shape of the unit cell was taken to be that found by structural relaxation at static conditions. We found that the unstrained perovskite structure persisted (meta) stably at all volumes V < 0.7VX. We are thus able to compare directly the simulated properties of liquid and perovskite over the entire lower mantle pressure-temperature regime.

Inspection of the equilibrated liquid structure shows that compression has a large influence on the atomic arrangement (Fig. 1). At the experimental zero-pressure volume, the liquid shows four-fold Si-O coordination, as does the stable crystalline phase at this volume, but without the elements of longer ranged order exhibited by the crystal, such as the infinitely long linear chains of tetrahedra that characterize pyroxenes. The liquid also shows a few non-four-fold-coordinated silicon atoms. At high pressure, the liquid structure is much more densely packed.

Fig. 1.

(Top) Mean Si-O coordination number at 3000 K (blue), 4000 K (green), and 6000 K (red) compared with that of the crystalline solidus phases plotted over their approximate range of stability. (Bottom) Distribution of Si-O coordination environments in MgSiO3 liquid at 3000 K. Pressures are those found in our FPMD simulations along the 3000 K isotherm for the liquid. Also shown are snapshots from the equilibrated portion of the simulations at 3000 K and V = VX (right) and V = VX/2 (left). In blue are the Si-O coordination polyhedra and in yellow the Mg ions. The primary simulation cell, which contains 80 atoms, is indicated by the black line.

In the simulations, the average Si-O coordination number increased smoothly and nearly linearly with compression, with no discernible influence of temperature, reaching a value of 6 at a pressure corresponding to the base of Earth's mantle (Fig. 1). The coordination change does not take place over a narrow pressure interval, as has been suggested (7). Such a rapid structural change would be expected to cause changes in the compressibility of the liquid (7) and in the slope of the melting curve (8) over a correspondingly narrow pressure interval, anomalies for which we find no evidence in our simulations.

The liquid equation of state can be accurately described by the Mie-Gr”neisen form Embedded Image(1) where Pc is the reference isotherm, taken to be at T0 = 3000 K, γ is the Gr”neisen parameter, and CV is the isochoric heat capacity. The latter quantities are computed from the simulations by noting that the internal energy E and pressure P vary linearly with temperature along isochores to within our uncertainty. The value of γ is then calculated at each volume from Embedded Image(2) and CV = (E/∂T)V. We find that CV decreases by about 10% from 4.12 ± 0.11 Nk at V = VX [the experimental value is 4.17 ± 0.18 Nk (9)] to 3.54 ± 0.11 Nk at V = VX /2, where N is the number of atoms and k is the Boltzmann constant. γ increases by a factor of three over the same range of volume, reflecting the increase in thermal pressure. All known mantle crystalline phases show the opposite behavior: γ decreases on compression (10).

The unusual behavior of γ in the liquid can be understood on the basis of the pressure-induced change in liquid structure (Fig. 2). Experimental data on crystalline MgSiO3 phases show that, although compression reduces the value of γ in each individual phase, polymorphic phase transformations have an even larger and opposing effect, with the higher coordinated phase showing the larger value of γ. As the coordination number of the liquid increases on compression, the liquid adopts values of γ characteristic of the more highly coordinated state. The behavior of the liquid is consistent with theoretical analyses of the influence of pressure-induced coordination changes on γ (11). Previous theoretical and experimental studies have also found that γ increases on isothermal compression in liquids (1214).

Fig. 2.

Equation of state of MgSiO3 liquid. FPMD results (open symbols) and Mie-Grüneisen equation of state (lines) fit to the FPMD results at 6000 K (red), 4000 K (green), and 3000 K (blue). Uncertainties in pressure, estimated on the basis of the appropriate statistics (35), are 1 to 2 GPa and are smaller than the size of the symbols on this scale. (Inset) Grüneisen parameter of liquid (open circles) and perovskite (open squares) from FPMD simulations, and experimental value (filled circle) for the liquid near its ambient melting point, calculated from the identity γ = αKT V/CV. Also shown are the Grüneisen parameter of MgSiO3 solidus crystalline phases at ambient conditions (filled squares) and under compression (dashed lines) (10). Experimental values of the thermal expansivity α, isothermal bulk modulus KT, and volume V are from (6), and isochoric heat capacity CV is from (9). The Mie-Grüneisen equation of state of the liquid is given by Eq. 1, with a third-order Eulerian finite strain expression for Pc and T0 = 3000 K, V0 = 1.13VX, K0 = 10.1 GPa, K0′ = 7.6, CV = 3.7Nk, and γ = 0.21 – 1.75(VV0)/VX.

The volume of the liquid closely approaches, but does not fall below, that of perovskite over the mantle pressure regime (Fig. 3). The volume difference is 4% at 6000 K and 140 GPa, as compared with 18% at the ambient-pressure melting point (6). Studies based on dynamic compression also find that the volume of fusion is markedly reduced at deep lower mantle conditions (15). The enthalpy difference between liquid and perovskite increases nearly linearly with increasing temperature and changes little with compression over the pressure regime of perovskite stability (25 < P < 140 GPa) (Fig. 3). The quantity ΔH/T varies little throughout this regime and has a value near 1.5 Nk, significantly greater than the value Nk ln2 often assumed for the entropy of fusion in high-pressure melting studies, and based primarily on one-component, close-packed materials. The larger value in our simulations can be understood on the basis of liquid structure: Throughout the lower mantle pressure regime, the liquid shows a wide range of different coordination environments, which adds a configurational component to the entropy that is absent in structurally simpler liquids (16). A recent FPMD study of MgO also found values of the entropy of melting ∼1.5 Nk at high pressure (17).

Fig. 3.

(Top) Volume and (Bottom) enthalpy of liquid (open circles and solid lines) and perovskite (filled circles and dashed lines) at 3000 K (blue), 4000 K (green), and 6000 K (red). For clarity, values at 3000 K and 4000 K are shifted downward in volume by 0.1 and 0.2, and in enthalpy by 500 and 1000 kJ mol–1, respectively. Statistical uncertainties are smaller than the size of the symbols. Lines with shading indicating estimated uncertainties are referred to the right-hand axes and show the liquid-perovskite volume and enthalpy differences. The volume difference is divided by the volume of the solid, and the enthalpy difference is divided by NkT. Lines are computed from the Mie-Grüneisen equation of state previously described (Fig. 2) and with H(P = 0,T0) = –3295 kJ mol–1 for the liquid. Parameters for the solid equation of state are T0 = 3000 K, V0 = 0.696VX, K0 = 162 GPa, K0′ = 4.3, CV = 3.1Nk, γ = 1.62 + 2.00(VV0)/VX, and H(P = 0,T0) = –3372 kJ mol–1. Previous results for perovskite are shown as plusses (32), crosses (33), and error bars (36). Open symbols show the volume and entropy of melting according to experiments at the ambient-pressure melting point (6, 9) and monatomic model liquids interacting through inverse power law repulsive pair potentials with the power indicated (37).

We used our simulations to predict the melting curve of MgSiO3 by integrating the Clausius-Clapeyron equation (18) Embedded Image(3) The quantities on the right-hand side are computed from our simulations over the entire lower mantle regime. The integration constant is set by assuming a single fixed point along the melting curve. We choose this point as 25 GPa, 2900 K near the low-pressure limit of perovskite melting, where many different experiments are mutually consistent.

Our predicted melting curve lies between previous experimental determinations, which diverge widely at pressures above 30 GPa (Fig. 4). The differences between theory and experiment can be understood in detail. Because we find ΔV > 0, our melting slope and melting temperature must be greater than experimental results with a vanishingly small melting slope, which requires ΔV = 0 (19, 20). Zerr and Boehler (21) argued that their much higher melting slope is consistent with the Lindemann law, a simple scaling argument that is valid if neither solid nor liquid undergo considerable structural change on compression (22). This picture of melting is incompatible with the large structural changes that we find in MgSiO3 liquid. Moreover, because the structural changes we see lead to more efficient packing, and therefore to reduced liquid-solid volume contrast, our melting slope and melting temperature should be correspondingly reduced. Some of the discrepancy among experiments may be due to differences in composition—Fe-free (21) versus ∼10% Fe (8, 19, 20)—but the wide divergence among diamond-anvil cell studies more likely reflects difficulties in detecting melt at high pressure (23). In this context, we note that a recent determination based on dynamic compression (15) agrees well with our predicted melting curve.

Fig. 4.

FPMD melting curve of MgSiO3 perovskite (red line with shading indicating uncertainties) compared with experimental determinations: HJ (19), KJ (8), IK (38), SH (20), ZB (21), SL (39), and ALAA (15). The fixed point assumed for the Clausius-Clapeyron integration is indicated. The dashed curve is the extrapolation by Zerr and Boehler (21) of their experimental results according to the Lindemann law. Thin long-dashed line shows the perovskite (pv) to post-perovskite (ppv) phase transformation (40). The cross is the pv-ppv-melt triple point computed using a classical pair-potential model “corrected” by comparing the location of lattice instabilities in quantum and classical molecular dynamics simulations (3).

The marked decrease in the volume of fusion of MgSiO3 that we find means that liquids are likely to be denser than coexisting solids in the deep mantle. In a multicomponent system, the density contrast between liquid and coexisting solid is a function of the volume contrast and the difference in chemical composition. In the shallow mantle, the volume contrast dominates, and liquids are less dense despite being enriched in heavy major elements such as Fe and Ca. In the deep mantle, where the volume contrast is small, the difference in composition is likely to be the determining factor. Experiments show that Fe and Ca continue to be preferentially incorporated in the liquid phase at least up to 25 GPa (24). For example, an apparent Fe-Mg partition coefficient between solid and liquid K = [(XFesol XMgliq)/(XFeliq XMgsol) ≈ 0.4, where XAα is the mole fraction of element A in phase or phase assemblage α, would be sufficient to overcome the volume difference between MgSiO3 liquid and perovskite and make the liquid denser.

The mantle solidus has been estimated to lie 800 to 1300 K below that of pure MgSiO3 (25). Combining this estimate with our ab initio value of the melting point at 136 GPa of 5400 ± 600 K yields a mantle solidus temperature of 3500 to 5200 K, which overlaps a range of previous estimates of the temperature at the core-mantle boundary (26, 27). Melting of normal mantle at its base cannot be ruled out and appears to be plausible.

The presence of neutrally or negatively buoyant melt in the deep mantle would have fundamental implications for our understanding of the Hadean Earth (28) and of the structure and chemical evolution of the coremantle boundary. If melting extended to the base of the mantle during accretion, perovskite may have floated in coexisting liquid as it crystallized out of the magma ocean, providing a possible means of large-scale chemical differentiation. In the present-day Earth, dense melt may provide a natural explanation for the ultralow velocity zone at the base of the mantle (29).

References and Notes

View Abstract

Navigate This Article