## Abstract

Carbon dioxide, generally considered as the second most abundant volatile component in silicate magmas, is expected to significantly influence various melt properties. In particular, our knowledge about its dynamical effects is lacking over most of Earth’s mantle pressure regime. Here, we report the first-principles molecular dynamics results on the transport properties of carbonated MgSiO_{3} liquid under conditions of mantle relevance. They show that dissolved CO_{2} systematically enhances the diffusion rates of all elements and lowers the melt viscosity on average by factors of 1.5 to 3 over the pressure range considered. It is remarkable that CO_{2} has very little or no influence on the electrical conductivity of the silicate melt under most conditions. Simulations also predict anomalous dynamical behavior, increasing diffusivity and conductivity and decreasing viscosity with compression in the low-pressure regime. This anomaly and the concomitant increase of pressure and temperature with depth together make these transport coefficients vary modestly over extended portions of the mantle regime. It is possible that the melt electrical conductivity under conditions corresponding to the 410- and 660-km seismic discontinuities is at a detectable level by electromagnetic sounding observation. In addition, the low melt viscosity values of 0.2 to 0.5 Pa⋅s at these depths and near the core-mantle boundary may imply high mobility of possible melts in these regions.

## INTRODUCTION

Carbon is widely considered to be one of the most important volatile elements in Earth’s interior (*1*–*4*). It is possible that the mantle carbon in significant amounts is actually incorporated in partial melts because of the increased solubility of carbon species including CO_{2}, in silicate melts under pressure (*5*–*9*) despite low carbon solubility in mantle silicate minerals (*10*). Dissolved CO_{2} is expected to influence the structural and transport properties of the host melt (*11*). The experimental data on various transport coefficients including diffusivity, viscosity, and electrical conductivity of silicate melt are confined to relatively low pressures over narrow temperature ranges (*12*–*14*). Computational methods were previously used to study carbon diffusion in silicate melts but only at selected low pressures (*7*, *15*). More studies are needed to answer a few key questions: How does carbon diffuse in compressed silicate melt? What is the influence of carbon on the chemical diffusivities and viscosity of silicate melt? How do these dynamical effects compare with those due to the hydrous component dissolved in the melt? It is also possible to link the transport properties of the liquids to their atomic structures. Assuming molecular CO_{2} and species dominance, it has been inferred that the measured carbon diffusivity is independent of the melt composition (*16*). However, the proportions and nature of these species are unknown in the melt stability field because the measurements were carried out for quenched glasses at room temperature.

We have recently performed first-principles molecular dynamics (FPMD) simulations to explore the structure, carbon speciation, and thermodynamics of carbon-bearing MgSiO_{3} liquid (*9*). Here, we extend these simulations to evaluate the transport properties of the MgSiO_{3} liquid containing 5.2, 16.1, and 30.5 weight % (wt %) CO_{2} (thereby allowing us to explore a wide compositional range) at temperatures from 2200 to 5000 K over the whole mantle pressure regime (see Materials and Methods). With knowledge of multiple properties including density, diffusivity, and viscosity of silicate melts, we can estimate useful parameters associated with magma dynamics. One such parameter is the melt mobility, which is a ratio of the melt-solid density contrast to the melt viscosity (*17*). This parameter determines the direction and magnitude of melt migration, that is, whether the melt will sink or float or stagnate in a partially molten mantle or a partially crystallized magma. The electrical conductivity arising from ionic transport in the melt can be related to the mantle electrical conductivity anomalies implied by magnetotelluric sounding data (*18*, *19*) as a possible indicator of partial melting at great depths (*20*). Our computational study of CO_{2}-rich MgSiO_{3} liquids represents a step forward in understanding the dynamical implications of the carbon dioxide component for deep mantle partial melts as well as for a large-scale magma ocean that could have existed in Earth’s early history.

## RESULTS

### Diffusion and viscosity at zero pressure

The self-diffusion coefficient results of a carbonated MgSiO_{3} liquid at ambient pressure show strong temperature dependencies with interspecies differences diminishing with increasing temperature (Fig. 1A). The predicted diffusivity sequence is Mg > C > O > Si, with Mg diffusion being six times faster than Si diffusion at the lowest temperature (2200 K) studied. Silicon always remains as the slowest species, although its diffusion rate varies the most with temperature. The overall temperature variation of carbon diffusivity is intermediate compared to the host melt element diffusion (Fig. 1A). Our results show that the melt viscosity also changes considerably with increasing temperature, decreasing by almost two orders of magnitude over the temperature range considered (Fig. 1B). Compared to the pure MgSiO_{3} liquid, the diffusion coefficients of all species increase, and the viscosity decreases in the presence of carbon dioxide at all temperatures (Fig. 1, A and B), but the differences are small. The overall dynamical effects of CO_{2} incorporation are thus to make the silicate melt more mobile and less viscous to some extent.

The calculated temperature variations of all self-diffusion (*D*_{α}) coefficients and viscosity (η) can be accurately described by the respective Arrhenius equations(1)where *E*_{α} and *E*_{η} are the activation energies for diffusion and viscosity, respectively, and *D*_{α0} and η_{0} are the corresponding preexponential factors (Table 1). Mg diffusion has a much smaller activation energy than the other three species. An interesting finding is that the higher the diffusivity, the lower the activation energy. Thus, all diffusivities increasingly converge to a similar value at some very high temperature—a tendency toward the compensation law of isokinetic behavior. The viscous activation energy takes an intermediate value compared to the diffusion activation energies. This means that all ions, including the network formers and modifiers, play a role in stress relaxations in silicate melt. The diffusion and viscosity activation energies of both the carbonated and pure MgSiO_{3} liquids shown in Table 1 are generally smaller than the experimentally inferred values for different melt compositions (*21*–*24*).

### Diffusion and viscosity under pressure

The general effects of pressure and temperature are suppression and enhancement of the melt mobility, respectively. However, these effects are neither monotonous nor complementary for the calculated diffusion coefficients of both the pure and carbonated MgSiO_{3} liquids. The log*D* versus *P* representations show significant nonlinearity, which is also sensitive to temperature (Fig. 2). The predicted complex *P*-*T* behavior of the diffusivity cannot be described by the standard Arrhenius relation assuming constant activation energy and volume values. The activation volume of species diffusion estimated as *V*_{α} = (*dD*_{α}/*dP*)_{T} appears to be dependent on both pressure and temperature. At and below 3000 K, the Si and O diffusivities tend to increase or remain unchanged initially with increasing pressure (Fig. 2). This represents anomalous diffusion because pressure is expected to systematically suppress the mobility of diffusing ions. A normal decreasing behavior is resumed after 10 GPa where the activation volume also becomes positive. On the other hand, carbon diffusion (shown in Fig. 2, lower right), similar to Mg diffusion, behaves normally, displaying a linear log*D*-*P* relationship over most of the pressure regime considered. The C diffusivity takes intermediate values compared to the host element diffusivities at all pressures. The behavior is quite different at higher temperatures where all species diffusivities decrease more rapidly initially with compression and gradually decrease thereafter. The diffusivity order and differences among the atomic species change with compression. Mg is generally the fastest species except at very high pressure, where O becomes the fastest species. Si remains as the least mobile element under all conditions.

Similar to the pure MgSiO_{3} liquid, the calculated viscosity of the carbonated MgSiO_{3} liquid is positively correlated with pressure over most of the pressure range considered (Fig. 3). At 4000 K, the viscosity increases monotonically with pressure by more than two orders of magnitude. At lower temperatures, the viscosity remains steady or even decreases initially with increasing pressure before it starts to increase rapidly on further densification. Anomalous viscous behavior can be clearly seen at 2200 K. This anomaly is consistent with the experimental observations of decreasing viscosity with pressure in highly polymerized silicate melts at low pressures (*25*). The calculated viscosity-pressure profiles mostly show a linear trend (in the log representation) except for a narrow low-pressure regime where the activation volume defined as *V*_{η} = (*d*η/*dP*)_{T} is sensitive to temperature. The activation volume appears to be nearly zero or even negative at ≤3000 K, but it takes large, positive values at higher temperatures. At pressures above 10 GPa, the viscosity-pressure slope decreases with increasing temperature, but the activation volume remains positive (Fig. 3). Temperature systematically suppresses the effects of pressure on the melt viscosity.

The addition of the CO_{2} component in the MgSiO_{3} liquid systematically enhances the diffusion rates of all host elements (Fig. 2) and lowers the melt viscosity (Fig. 3). However, it does not greatly affect the dynamic anomalies predicted in the low-pressure regime. The initial increases in the Si and O diffusivities and the initial decrease in the melt viscosity with compression are usually attributed to the presence of fivefold oxygen-coordinated silicon species acting as transition states (*26*, *27*). The carbon-induced changes in ionic diffusivities and viscosity tend to remain roughly the same under most conditions except for very high pressures where they become notably large (Figs. 2 and 3). For instance, the viscosity of the melt containing 16.1 wt % CO_{2} is almost half of the pure melt viscosity above 100 GPa along the 4000 K isotherm. The comparison of the transport coefficients between 5.2, 16.1, and 30.5 wt % CO_{2} concentrations in the melt hints toward an overall positive (negative) correlation between the diffusivity (viscosity) and the CO_{2} content (Fig. 4). Our finding that the effects on the transport properties are stronger at higher CO_{2} concentration is generally consistent with the experimental data on carbon-bearing silicate melts (*13*).

### Electrical conductivity

Mobile ions carrying finite charges are expected to contribute to the electrical conductivity of silicate melt. The estimated electrical conductivity (σ) behaves qualitatively similar to the diffusivities with respect to temperature and pressure (Fig. 5). At ambient pressure, σ is positively correlated with temperature. The melt electrical conductivity also shows nonmonotonic pressure evolution at all temperatures—increasing with compression initially, reaching broad maximum, and then decreasing gradually on further compression (Fig. 5). As such, the overall pressure variations of σ are smaller than those of the self-diffusion and viscosity coefficients. It is interesting to note that the calculated electrical conductivity of the carbonated MgSiO_{3} liquid is similar (within the computational uncertainty) to that of the pure phase under most conditions. This prediction is at odds with the general notion that the volatile components including CO_{2} and H_{2}O accelerate the dynamics, thereby making the melt more conductive. For instance, recent experimental data on hydrous carbonated basaltic melts show that the addition of 10 to 25 wt % CO_{2} can increase the electrical conductivity by an order of magnitude (*28*). Small effects of CO_{2} on the electrical conductivity have also been predicted by the previous computational study of carbonated basaltic and kimberlitic melts in the CaO-MgO-Al_{2}O_{3}-SiO_{2} system (*15*). Only at very high pressures does the dissolved CO_{2} tend to increase the melt conductivity.

The first-principles electrical conductivity values of the pure and carbonated MgSiO_{3} melts at the lowest temperature (2200 K) studied are larger than the available low-pressure experimental data in the temperature range 1673 to 2073 K for basaltic and ultramafic silicate melts with or without the carbon dioxide component (*28*–*32*). The discrepancies can be attributed to higher temperatures in simulations than in experiments, but the compositional effects might be significant as well. The increased electrical conductivity values for basaltic melt in the experiments (*28*) actually represent the combined effects of CO_{2} and H_{2}O. On an independent basis, the effects of CO_{2} on electrical conductivity of dry basalt melt are presently unknown, although H_{2}O is experimentally shown to significantly increase electrical conductivity (*31*). Our calculated values of the hydrous MgSiO_{3} at 0 GPa are systematically larger than the corresponding pure melt values (Fig. 5). It may be possible that the addition of CO_{2} in silicate melts facilitates transport of alkali cations Na and K, which can markedly raise the conductivity (*31*). Note that alkali cations are absent in our simulations unlike experimental compositions. Our calculated results also show that the electrical conductivity of carbonated MgSiO_{3} liquid does not show systematic dependence on CO_{2} content (Fig. 5), consistent with the experimental inferences for a low-concentration regime (*14*, *16*). The previously calculated conductivities of basaltic and kimberlitic melts containing 20 wt % CO_{2} at 2073 K and 12 GPa are lower than our results for the carbonated MgSiO_{3} liquid at 2200 K. The effective ionic charges used by those studies (*15*, *33*) are also smaller than our values derived from first-principles simulations (see Materials and Methods).

## DISCUSSION

A recent structural analysis of the carbonated MgSiO_{3} liquid (*9*) shows that, at low pressures, carbonate () groups prevail over CO_{2} molecules, as experimentally inferred for other silicate compositions (*34*). The molecular species–to–carbonate species ratio increases as temperature is raised. The CO_{2} molecules are free in the sense that they are exclusively linked to Mg via the C–O–Mg connections. On the other hand, the carbonate groups are also bound to the silicate polyhedral network via the C–O–Si connections. Some carbonates exist as free units; that is, none of their three oxygen atoms is bonded to any silicon. Visualization of the temporal behavior of these species confirms that carbon diffuses as free CO_{2} molecules and free carbonate groups, with polyhedral carbonates mostly serving as the transition states at low pressure and low temperature (Fig. 6). Carbon self-diffusion thus mostly occurs as a multi-atom movement, which involves Mg–O bond (breaking/formation) events much more frequently than C–O or Si–O bond events. This finding appears to contradict the general assumption that the carbonate groups are strongly bound to the silicate network and remain immobile (*14*, *16*, *35*). The degree of polyhedral connection considerably increases with compression. At high pressure, CO_{2} is incorporated as polyhedral bound carbonates, tetrahedral coordination (CO_{4}), and polymerized di-carbonates. These species are polymerized with the silicate framework (*9*) and are no longer able to move as stable units. All bond events including C–O and Si–O bond events are involved in self-diffusion processes.

The melt electrical conductivity (σ) can be related to the self-diffusion coefficients of *n*_{α} species via the Nernst-Einstein equation: , where ρ_{i} is the number density of species *i* carrying charge *z*_{i}. The Haven ratio *H*_{R} is a correlation factor for the ionic motions. For *H*_{R} = 1, ionic movements are independent of each other. We find that the values of σ_{NE} estimated with *H*_{R} = 1 and the effective ionic charges used are much larger than the directly calculated conductivity σ at ambient pressure. For instance, the σ_{NE}/σ ratio is about 3 at 0 GPa and 3000 K (implying that *H*_{R} ≈ 3) and becomes smaller at lower temperatures. This suggests that the motions of ions are not cooperative in causing the net charge transfer and that the cation-anion ionic dissociation is not very effective. Charge neutral CO_{2} molecules and doubly negative charged carbonate groups (also serving as transition states) likely move as stable units, as discussed in the previous paragraph. This also means that the concentration of charge carriers is reduced in the carbonated silicate melt, which tend to provide similar conductivity to the pure melt despite their higher mobility (Fig. 5). The carbonate () groups whose abundance increases initially with increasing pressure (*9*) appear to contribute to the electrical conductivity more as the melt is compressed. The difference between two conductivities decreases rapidly with increasing pressure, and the σ_{NE}/σ ratio eventually becomes smaller than 1 (implying *H*_{R} ≤ 1) above 50 GPa. All atoms diffuse at similar rates (Fig. 2), perhaps independent of each other, such that the diffusive motions of the charge carriers at the individual level are mostly operative in a highly compressed melt. Thus, our analysis implies that CO_{2} speciation plays a crucial role in the electrical conductance of the silicate melt.

Our results suggest that the presence of CO_{2} accelerates silicate melt dynamics to some extent by increasing the diffusion coefficients of the host ions (Fig. 2) and decreasing the melt viscosity (Fig. 3). However, these effects are smaller than those due to the hydrous component (Fig. 4). The addition of 16.1 wt % CO_{2} increases the Si/O diffusivity and decreases the melt viscosity by factors of 1.5 to 3 depending on the *P-T* condition, whereas the addition of 10 wt % H_{2}O changes these transport coefficients by factors of 3 to 10 (*26*, *36*). In addition, H diffusion rate is larger than the diffusion rates of all species including carbon by one or two orders of magnitude. The large differences between the carbonated and hydrous silicate melts can be linked to the speciation of the CO_{2} and H_{2}O components and how they influence the melt structure. Highly mobile protons can contribute substantially to the electrical conductivity of silicate melt (*27*, *37*), as shown in Fig. 5, which is consistent with the experimental data (*28*). Nevertheless, our first-principles simulations predict that the electrical conductivity of anhydrous silicate melt with or without CO_{2} may remain largely undiminished, taking high values over the extended part of the mantle regime [for example, lying above 10 S/m at pressures up to 25 GPa along the 2200 K isotherm and around 200 S/m near the core-mantle boundary (CMB) pressure at 4000 K]. Thus, it is possible that the incipient melting of the mantle at a depth causes high electrical conductivities in the range 0.1 to 2.0 S/m, as found by the electromagnetic sounding observations of the upper mantle and transition zone (*18*–*20*).

To assess the implications of the calculated density and viscosity results for magma dynamics, we can estimate the melt mobility parameter defined as Δρ/η = (ρ_{M} − ρ_{S})/η. Here, ρ_{M} is the density of the melt for which we use the recently estimated density-pressure profiles of the (Mg_{1–x},Fe_{x})SiO_{3} liquid containing 5 wt % CO_{2} with *x* = 0.22 at 2200 K and *x* = 0.27 at 4000 K (*9*). For the density ρ_{S} of the surrounding solid rock, we take the mantle as representative (*38*). For the melt viscosity η, we use the first-principles results. Corresponding to the 660-km depth, we have ρ_{M} = 4.18 g cm^{−3} and η = 0.5 Pa⋅s at 2200 K and 23 GPa and ρ_{S} = 4.38 g cm^{−3}. The mobility value is −0.4, whose magnitude rapidly decreases with increasing depth. This means that any melt in the uppermost part of the lower mantle rapidly floats toward the 660-km depth where the melt becomes gravitationally stable. Our estimated mobility ratio takes a positive value (0.2) in the lowermost mantle (using ρ_{M} = 5.59 g cm^{−3}, ρ_{S} = 5.55 g cm^{−3}, and η = 0.25 Pa⋅s at 4000 K and 135 GPa), suggesting that the melt there sinks toward the CMB. Thus, melts tend to rapidly migrate out of their source regions, perhaps forming narrow concentrated regions, such as melt pockets or layers, at the 660-km depth and CMB. The possibility of this excessive melt accumulation at the boundary between the lithosphere and asthenosphere (80- to 100-km depths) was previously suggested (*17*). Our analysis hints toward the possibility of dense silicate melts at depths corresponding to the top and bottom of the lower mantle.

## MATERIALS AND METHODS

FPMD simulations of the liquid systems under consideration were performed in a canonical NVT ensemble using the local density approximation and projector augmented wave potentials, as implemented in VASP (Vienna Ab initio Simulation Package) (*39*). The cubic supercells used to sample the pure and 16.1 wt % CO_{2}-bearing compositions contained 32MgSiO_{3} and 32MgSiO_{3} + 14CO_{2}, respectively, as in our recent study (*9*). Our simulations covered a wide range of volume (*V*/*V*_{X} = 1.5 to 0.5, where *V*_{X} = 2067.1 Å^{3} for the melt containing 16.1 wt % CO_{2}), corresponding to the whole mantle pressure regime (0 to 140 GPa). At each volume, initial melting and equilibration were performed at 8000 K so that no reminiscence memory of solid phase would be retained by the simulated liquid phase. The system was then quenched isochorically to lower temperatures 5000, 4000, 3000, 2500, and 2200 K. We also simulated 32MgSiO_{3} + 4CO_{2} (5.2 wt %) and 16MgSiO_{3} + 16CO_{2} (30.5 wt %) systems under selected conditions to explore the effects of the carbon content. A plane wave cutoff energy of 400 eV with gamma-point Brillouin zone sampling was considered. The Pulay stress was estimated as the difference in the calculated pressures using this cutoff and a larger cutoff of 600 eV at each volume. Its value increases from 2 to 7 GPa almost linearly with decreasing *V*/*V*_{X} over the entire compression range considered. The run durations with a time step of 1 fs ranged from 50 to 300 ps, depending on the thermodynamic points considered to achieve acceptable convergence.

The self-diffusion coefficient *D*_{α} of species α (Mg, Si, O, and C) was calculated from the mean square displacement (MSD) using the Einstein relation(2)where ’s represent the atomic trajectories and 〈 〉 represents an average at time *t* to be taken with different time origins (*t*_{0}). Diffusion coefficient was obtained in the diffusive regime, where MSD in Eq. 2 varies linearly with time (*40*, *41*). We evaluated the melt electrical conductivity from the simulated atomic position time series data using the following relation(3)Here, *z*_{i} represents the charge on each ion for which we use the corresponding formal charges of *z*_{Mg} = +2e, *z*_{Si} = +4e, *z*_{O} = −2e, and *z*_{C} = +4e. We can actually estimate the effective ionic charges from the Bader analysis of the simulated electronic charge distribution on a three-dimensional grid (*42*, *43*). In the melt, the charge on each ionic species varies from site to site (within the same snapshot), as well as with time (between snapshots), depending on their local surroundings. Because these variations are small, we can take the appropriate average for each ion type. The effective charges thus obtained from the first-principles simulations are *z*_{Mg} = +1.6e, *z*_{Si} = +3e, *z*_{O} = −1.4e, *z*_{C} = +1.9e, and *z*_{H} = +0.6e.

Similarly, the shear viscosity (η) was calculated from the time integral of the stress tensor autocorrelation function (ACF) using the Green-Kubo relation(4)where *P*_{αβ} is the deviatoric stress tensor that represents the symmetrized traceless part of the stress tensor computed at every FPMD step. Averaging was performed by calculating the integrand in Eq. 4 at time *t* with different time origins (*t*_{0}). The viscosity was estimated from the integral computed for long time, that is, beyond the point where the ACF decays to zero (*26*, *44*). The typical uncertainties in the evaluation of the self-diffusivities, conductivity, and viscosity vary between ±10 and ±30% with larger errors at higher pressure and lower temperature due to slow dynamics (*41*), as well as for electrical conductivity and viscosity due to their collective nature.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**High-performance computing resources were provided by Louisiana State University (www.hpc.lsu.edu) and the Louisiana Optical Network Institute.

**Funding:**This work was funded by NSF (EAR-1426530).

**Author contributions:**B.B.K. conceived and supervised the research and D.B.G. performed most simulations. Both D.B.G. and B.B.K. analyzed the results, designed the figures, and wrote the manuscript.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper. Additional data related to this paper may be requested from the authors.

- Copyright © 2017 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).