Viscosity of MgSiO3 Liquid at Earth’s Mantle Conditions: Implications for an Early Magma Ocean

See allHide authors and affiliations

Science  07 May 2010:
Vol. 328, Issue 5979, pp. 740-742
DOI: 10.1126/science.1188327


Understanding the chemical and thermal evolution of Earth requires knowledge of transport properties of silicate melts at high pressure and high temperature. Here, first-principles molecular dynamics simulations show that the viscosity of MgSiO3 liquid varies by two orders of magnitude over the mantle pressure regime. Addition of water systematically lowers the viscosity, consistent with enhanced structural depolymerization. The combined effects of pressure and temperature along model geotherms lead to a 10-fold increase in viscosity with depth from the surface to the base of the mantle. Based on these calculations, efficient heat flux from a deep magma ocean may have exceeded the incoming solar flux early in Earth’s history.

Silicate liquids likely played a crucial role in terrestrial mass and heat transport in Earth’s history. Molten silicates would have controlled the dynamics of the predicted magma ocean [a largely or completely molten mantle that is expected during Earth’s earliest stages (1)] and continue to influence the transport of modern magmas at the present. If such a magma ocean existed, the rates of initial thermal evolution (via convection) and chemical evolution (via crystal settling and melt percolation) of Earth’s interior would be primarily controlled by the melt viscosity (2). The ability of melts to carry xenoliths from great depths in the mantle (3) also depends on the melt viscosity, in addition to melt composition. Moreover, melts are considered to be responsible for the ultralow velocity zone (ULVZ) in the deep mantle detected by seismology (4, 5).

Despite their importance, transport properties, including the viscosity of molten silicates, are unknown over almost the entire mantle pressure regime, which reaches 136 GPa at the core-mantle boundary. Because of experimental difficulties, the viscosity of MgSiO3 liquid, the dominant composition of Earth’s mantle, has only been measured at ambient pressures (6). In fact, viscosity measurements of any silicate melts have been limited to relatively low pressures (<13 GPa) (712). In many silicate liquids, the viscosity depends non-monotonically on pressure over the range that has been measured, making extrapolations highly uncertain. Theoretical computations serve as a complementary approach. Previous calculations were primarily based on atomistic models (1315), which permit much faster computation but have the disadvantage of being based on empirical force fields, the forms of which are uncertain. On the other hand, the first-principles approach is more robust because it makes no assumptions about the nature of bonding or the shape of the charge density and is thus in principle equally applicable to the study of a wide variety of materials problems, including liquids. We previously calculated the structure and thermodynamic properties of MgSiO3 and MgSiO3-H2O liquids from first principles (16, 17), finding good agreement with extant experimental data over the entire mantle pressure-temperature regime. Unlike these equilibrium properties, the transport properties such as viscosity require much longer simulation (18).

Here, we determine the viscosity of two key liquids over the entire mantle pressure regime from density functional theory (18). MgSiO3 serves as an analog composition for a magma ocean, whereas MgSiO3-H2O liquid allows us to explore the role of melt composition, focusing on H2O as the component that is known to have the largest influence on the viscosity at low pressure (19). The shear viscosity (η) was calculated by using the Green-Kubo relationη=V3kBT0i<jσij(t+t0).σij(t0)dt (1)where σij (i and j = x, y, z) is the stress tensor, which is computed directly at every time step of the simulation, V is volume, kB is the Boltzmann constant, T is temperature, t is time, and t0 represents the time origin. The shear-stress autocorrelation function (the integrand of Eq. 1) decays to zero more slowly at lower temperature and higher pressure, requiring longer simulation runs (Fig. 1). We find that the integral values, and hence the computed viscosity, converge over time intervals much shorter than the total simulation durations (Fig. 1, inset). The fact that the shear stress autocorrelation function decays to zero within the time scale of our simulations means that the Maxwell relaxation time of silicate liquids (20) remains much shorter than seismic periods over the entire mantle regime and that seismic wave propagation through melts that may exist in the ULVZ will occur in the relaxed limit. We further confirm that the simulated system is in the liquid state at each pressure-temperature condition by examining the mean-square displacements (fig. S1) and radial distribution functions (fig. S2). Our approach is expected to be more robust than the commonly used indirect approach of estimating the melt viscosity from the self-diffusion coefficient via the classic Erying relation (20, 21). The validity of the Eyring relation as applied to silicate liquids has been questioned on the basis of experiments (11).

Fig. 1

Time convergences of the calculated stress autocorrelation function (ACF) and viscosity (inset) of MgSiO3 melt (without water) at different conditions. The run durations are 18 ps (VX, 6000 K, 7.5 GPa), 60 ps (VX, 3000 K, 1.8 GPa), 72 ps (0.7VX, 3000 K, 25 GPa), and 172 ps (0.5VX, 4000 K, 135 GPa), where VX is the reference volume (38.9 cm3 mol–1).

Over most of the pressure range of our investigation, viscosity increases with increasing pressure (Fig. 2A). The calculated viscosity increases by a factor of ~140 for anhydrous silicate melt over the entire mantle pressure regime at 4000 K. The activation volume Vη=(dlnη/dP)T varies systematically over most of this range, tending to decrease with increasing pressure and with increasing temperature. At the lowest temperature and pressure, we find that the viscosity behaves anomalously: decreasing with increasing pressure initially, reaching a minimum value near 5 GPa at 3000 K, and then increasing on further compression. Low-pressure experimental studies have found viscosity decreasing with increasing pressure in highly polymerized silicate melts (8, 11, 23). We attribute the initial decrease in viscosity with increasing pressure to the presence of fivefold coordinated silicon, which acts as a transition state accommodating viscous flow (16, 24). The variation of viscosity in this anomalous regime is small compared with the total variation in viscosity over the mantle pressure-temperature range. The calculated viscosities show large and systematic deviations from Arrhenian behavior (Fig. 2B). The activation energy decreases with increasing temperature, consistent with the behavior of moderately fragile liquids (25).

Fig. 2

Calculated viscosity (η) of anhydrous (anhy, solid symbols) and hydrous (hy, open symbols) MgSiO3 melts. Our results are compared with experimental (Expt) data (6) at lower temperatures and ambient pressure for anhydrous MgSiO3 liquid. Error bars indicate the statistical uncertainties. (A) Pressure variations along 3000 K (circles), 4000 K (squares), and 6000 K (diamonds) isotherms. The anhydrous results can be represented by the modified VFT (Vogel-Fulcher-Tammann) equation (35): Embedded Image, with T0 = 1000 K. Also shown are the results at 3000 K (14) and 4000 K (15) from previous molecular dynamics (MD) studies of anhydrous liquid based on semi-empirical pair potentials. (B) Temperature variations at the reference volume (VX) together with experimental data (6) represented by the VFT equation, Embedded Image, where A = 0.00033 Pa s, B = 6400 K for anhydrous liquid and A = 0.00024 Pa s and B = 4600 K for hydrous liquid.

Silicate melt with 10 weight percent H2O is two to four times less viscous than the anhydrous melt at all pressure-temperature conditions studied (Fig. 2). We have previously shown that the self-diffusion coefficients of the hydrous liquid are systematically higher than those of the anhydrous liquid (26). The region of anomalous pressure dependence of the viscosity and diffusion are weak or absent in the case of hydrous silicate liquid. Our first-principles results confirm that the dynamical enhancement (smaller viscosity and larger diffusivity) occurs in hydrous silicate liquid because water systematically depolymerizes the melt structure (27). The mean O-Si and Si-Si coordination numbers decrease in the presence of water: The hydrous values vary from 1.1 to 1.7 on compression (compared with anhydrous value of 1.4 to 2) and 1.8 to 4.7 (compared with anhydrous value of 2.5 to 5.5), respectively, over the compression range studied (17, 26).

The viscosity of silicate melts increases modestly along temperature profiles characteristic of Earth’s interior because of the competing effects of pressure and temperature (Fig. 3). For example, along a slightly super-liquidus magma ocean isentrope (28), the viscosity of anhydrous melt increases by a factor of 10 from the surface to the core-mantle boundary. The variations are similar in size but non-monotonic along the estimated mantle solidus and liquidus (28). The non-monotonic variation of the viscosity along these curves is due to the rapid increase in temperature with increasing pressure at low pressure. The viscosity profiles of the hydrous melt show similar variations but are systematically shifted downward (Fig. 3).

Fig. 3

Predicted viscosity for magma ocean based on anhydrous silicate liquid results (using the modified VFT relation shown in Fig. 2) along different temperature profiles (inset) from (28): magma ocean isentrope, mantle liquidus, mantle solidus, and melting curve of pure MgSiO3 perovskite (dashed line). The hydrous result is shown only along the mantle liquidus (gray line). The symbol represents a viscosity value at mid-mantle condition.

Our results provide a fundamental basis for any dynamical model of magma ocean evolution. To illustrate, we use our viscosity results and those of previous ab initio simulations of silicate melts to estimate critical dynamical parameters. For a completely molten mantle such as one that may have occurred early in Earth’s history (1), the estimated Rayleigh and Prandtl numbers (29) are ~6 × 1030 and ~60, respectively, for the viscosity value of 0.048(10) Pa s for anhydrous MgSiO3 liquid at mid-mantle condition (70 GPa and 4000 K along the magma ocean isentrope, Fig. 3). The Rayleigh number lies in the regime of turbulent convection: The presence of turbulence may substantially influence the settling of crystals as they form upon cooling. The surface heat flux, F ~ 6 × 104 W m−2, estimated from mixing length theory far exceeds the incoming solar flux (30) and suggests that the surface temperature was set by heat exchange of the magma ocean with a dense silicate atmosphere rather than by solar radiation balance (2). This value of F implies a cooling time for the magma ocean ~20 ky (30). In fact, a number of processes are likely to increase the cooling time of the magma ocean substantially, including crystallization, which is predicted to initiate in the mid-mantle (28) and to separate the magma ocean into upper and basal layers (31). The evolution at this stage also depends strongly on the viscosity, which will set the time scale for buoyancy-driven motion of crystals and liquid that can lead to chemical differentiation. The direction of motion will be set by the crystal-liquid density contrast, the sign of which varies with pressure and temperature. Indeed, crystals are expected to float near the base of the mantle (16), producing a buoyantly stable basal magma layer that may be long-lived (31).

Supporting Online Material

Materials and Methods

Figs. S1 and S2

References and Notes

  1. Computations were performed by using the VASP software (32) with the local density approximation and ultrasoft pseudopotentials as before (16, 17, 26). Methods are available as supporting material on Science Online.
  2. The Rayleigh and Prandtl numbers are Ra=[αρg(TMTS)L3]/(κη) and Pr=η/(ρκ). We have adopted the density (ρ = 4410 kg m−3) and thermal expansivity (α = 2.6 × 10−5) (16) and assumed the depth scale L = 3000 km, acceleration due to gravity g = 10 m s−2, thermal conductivity k = 1.2 W m−1 K−1 (33), thermal diffusivity κ = k/(ρcP) where the specific heat cP = 1660 J kg−1 K−1 (16), a mantle potential temperature TM = 2500 K, the lowest temperature at which the mantle will be completely molten, and a surface temperature TS = 1000 K set by a dense atmosphere (2).
  3. The surface heat flux is F=0.22k(TMTS)Ra2/7Pr1/7L1 (34) and the cooling time is τcool=TMcPMM/4πR2F, where MM is the mass of the mantle and R is the radius of Earth.
  4. This work was supported by NSF (EAR-0809489) and the UK National Environmental Research Council (NE/F01787/1). Computing facilities were provided by the Center of Computation and Technology at Louisiana State University. The authors thank J. Brodholt, M. Ghiorso, and S. Karato for useful comments and suggestions.
View Abstract

Stay Connected to Science

Navigate This Article