The fate of carbon dioxide in water-rich fluids under extreme conditions

See allHide authors and affiliations

Science Advances  12 Oct 2016:
Vol. 2, no. 10, e1601278
DOI: 10.1126/sciadv.1601278


Investigating the fate of dissolved carbon dioxide under extreme conditions is critical to understanding the deep carbon cycle in Earth, a process that ultimately influences global climate change. We used first-principles molecular dynamics simulations to study carbonates and carbon dioxide dissolved in water at pressures (P) and temperatures (T) approximating the conditions of Earth’s upper mantle. Contrary to popular geochemical models assuming that molecular CO2(aq) is the major carbon species present in water under deep Earth conditions, we found that at 11 GPa and 1000 K, carbon exists almost entirely in the forms of solvated carbonate (Embedded Image) and bicarbonate (Embedded Image) ions and that even carbonic acid [H2CO3(aq)] is more abundant than CO2(aq). Furthermore, our simulations revealed that ion pairing between Na+ and Embedded Image/Embedded Image is greatly affected by P-T conditions, decreasing with increasing pressure at 800 to 1000 K. Our results suggest that in Earth’s upper mantle, water-rich geofluids transport a majority of carbon in the form of rapidly interconverting Embedded Image and Embedded Image ions, not solvated CO2(aq) molecules.

  • First-principles molecular dynamics
  • deep carbon cycle
  • supercritical water
  • CO2


The threat of global warming and climate change makes the understanding of Earth’s carbon cycle a critical and pressing task. In Earth’s lithosphere, sedimentary carbonates are considered to be the largest carbon hosts (1), with carbon transport between Earth’s surface and interior (2) continuously occurring over millions to billions of years (3). For example, carbon belonging to organic matter and carbonate minerals may go deep into Earth’s mantle within subduction zones, and it may be brought back to the surface by volcanism through CO2 degassing (4). The deep carbon cycle substantially influences the carbon budget near Earth’s surface, which, in turn, affects our energy consumption and global climate change (4). However, whether carbon accumulates in our planet’s interior is still a subject of debate, and the chemical reactions leading to the predominant species of carbon in deep Earth are not yet well understood; hence, the means by which carbon transport occurs within Earth has not yet been established (2).

Carbon may be transported by crustal and mantle fluids in subduction zones (5), where water is a major transport medium (6). However, the forms of dissolved carbon in water under deep Earth conditions are poorly known. In Earth’s upper mantle, the pressure (P) may reach ~13 GPa and the temperature (T) ~1700 K (7). It is challenging to study aqueous solutions of CO2 and carbonates in the laboratory at these conditions, because water becomes highly corrosive (8, 9) and it is difficult to obtain clear vibrational spectroscopic signals, which are commonly used to characterize aqueous solutions under ambient conditions. Additionally, the interpretation of existing Raman spectra measured at high P-T is controversial, because it is yet unclear how to relate the solute concentration to the observed signals under supercritical conditions (10, 11). Experimental data have been obtained either at high P (HP) or at high T (HT), but it has not yet been possible to reach both HP and HT, as shown in Fig. 1. Therefore, our understanding of the chemistry of aqueous carbon is limited to P-T conditions of shallow mantle environments (12). Because of the lack of experimental data, traditional geochemical models have been based on several assumptions, one of which is that CO2(aq) is the major dissolved carbon species in supercritical water under oxidizing conditions [for example, (1315) and references therein].

Fig. 1 High-pressure and high-temperature conditions reached in experiments on dissolved carbon in supercritical water, together with the conditions simulated in the present work: Frantz (10), Caciagli and Manning (30), Martinez et al. (19), Facq et al. (31, 36), and Schmidt (11).

The melting curve of ice is from the study by Datchi et al. (51). The shaded area shows the P-T conditions of the upper mantle (7).

Here, by conducting extensive first-principles molecular dynamics (MD) simulations, we studied Na2CO3 and CO2 aqueous solutions at 11 GPa and 1000 K, that is, at pressures similar to those at the bottom of Earth’s upper mantle. We note that ab initio simulations, which do not require any experimental input, have been shown to be powerful tools to study chemical reactions under extreme conditions (16). However, they are computationally demanding and cannot yet be used to carry out a detailed study of the whole P-T range in Earth’s mantle. Instead, the present study focuses on general properties of dissolved carbon in water under extreme conditions; such conditions could not be addressed by simple models, which do not explicitly take into account the molecular nature of water and hence the way it dissociates under pressure. Contrary to the assumption that CO2(aq) is the major carbon species under oxidizing conditions, we found that most of the dissolved carbon is in the form of Embedded Image and Embedded Image, with continuous transformation of the two species into each other on the scale of picoseconds. Surprisingly, even H2CO3(aq) is more abundant than CO2(aq) when CO2 is directly dissolved into water at 11 GPa and 1000 K. We also found that ion pairing occurs between Na+ and Embedded Image or Embedded Image. Our study confirms that the molecular structure of supercritical water is markedly different from that of the liquid under ambient conditions and shows that the fast chemical dynamics in aqueous solutions under extreme conditions is key to understanding dissolved carbon in deep geofluids. Our findings suggest that at the bottom of the upper mantle, solvated CO2(aq) is hardly present in water, with most of the species being carbonate and bicarbonate ions. We expect these highly active ionic species to be involved in the carbon recycling process in subduction zones in deep Earth and to play an important role in the transport of carbon.


Validation of computational strategy

Before carrying out simulations at Earth’s mantle conditions, we validated our computational strategy; in particular, we tested the accuracy of exchange-correlation (xc) functionals used in density functional theory by simulating a carbonate solution at 0.2 GPa and 823 K, for which experimental spectroscopic studies and geochemical models are available (17). We considered 1 Na2CO3 molecule and 62 water molecules (molality, 0.9 m) in a cubic simulation box with a side length of 14.42 Å. According to the equations of state by Zhang and Duan (14, 18), which are fitted to experimental data, the density of water under these conditions corresponds to a pressure of ~0.2 GPa. We chose sodium as a representative cation because Na2CO3 has a large solubility in supercritical water (11, 19, 20); Na is the seventh most abundant element in Earth’s mantle (21).

Figure 2A shows the mole percents of CO2(aq), Embedded Image, Embedded Image, and H2CO3(aq) species as functions of simulation time in first-principles MD simulations carried out with semilocal [PBE (Perdew-Burke-Ernzerhof) (22)] and hybrid [PBE0 (23)] functionals. The overall mole percents in our extensive MD simulations are summarized in Table 1. Using PBE, we found that bicarbonate is the dominating carbon species in the solution, with ~82% (mole percent of total dissolved carbon; see also Methods) bicarbonate ions and ~17% carbonate ions over a ~130-ps MD simulation. The amounts of CO2(aq) and H2CO3(aq) were negligible. The lifetime of bicarbonate ions exceeded 10 ps in our simulations, whereas the lifetime of Embedded Image was much shorter.

Fig. 2 Mole percents of CO2 (orange), Embedded Image (red), Embedded Image (black), and H2CO3 (blue) as functions of simulation time in first-principles MD simulations.

(A) Na2CO3 solution [0.9 m (molality)] at 0.2 GPa and 823 K. (B) Na2CO3 solution (0.9 m) at 11 GPa and 1000 K. (C) CO2 solution (0.9 m) at 11 GPa and 1000 K. Calculations carried out with two xc functionals: PBE (22) and PBE0 (23) are compared.

Table 1 Carbon species in the 0.9 m Na2CO3 and 0.9 m CO2 solutions in first-principles MD simulations.

The concentrations are shown as mole percents of the total dissolved carbon. The mole percents of Embedded Image and Embedded Image include those of ion pairs formed by them, respectively. The experimental results (Expt.) are from the Raman data in the study by Frantz (10) (see text). N/A, not applicable.

View this table:

It is known that generalized gradient approximations (GGAs) may not be accurate in the description of solvated doubly charged anions because of the charge delocalization error (24). For example, in sulfuric acid aqueous solutions under ambient conditions, Wan et al. found that the concentration of Embedded Image was seriously overestimated when using the PBE functional, whereas the hybrid functional PBE0, which better describes the charge localization of solvated anions (25), gave results in qualitative agreement with experiments (2628). Thus, we compared PBE and PBE0 simulations in order to extract robust trends from our results. We started a PBE0 simulation from a ~40-ps-long PBE simulation. As shown in Fig. 2A, Embedded Image ions immediately transformed into Embedded Image anions. In the subsequent simulation, the total mole percent of Embedded Image was 27%, indicating that PBE overestimates the concentration of Embedded Image with respect to PBE0. However, in both PBE and PBE0 simulations, we found that few CO2(aq) or H2CO3(aq) molecules were present in the solution at 0.2 GPa and 823 K. Thus, the absence of these species in the solution was considered to be a robust result of our first-principles MD simulations.

Using experimental Raman spectroscopy (10), Frantz studied a 1 m K2CO3 solution under the same P-T conditions as those considered here. The concentration of a given species (Ci) was related to the measured Raman intensity (Ii) by Ii = γiCi, where γi is the Raman scattering cross section. On the basis of the concentrations calculated by geochemical models, Frantz estimated that the ratio Embedded Image decreased by ~75% when increasing T from 523 to 823 K (10). Recently, after reexamining Frantz’s spectroscopic data, Schmidt argued that this ratio should change by less than 10% in that temperature range (11), indicating that the species concentrations obtained by geochemical models in Frantz’s study (10) may be questionable. Our theoretical concentrations are closer to the results obtained by the γ ratio proposed by Schmidt (11) than by Frantz (10). Using Schmidt’s γ ratio and the published Raman spectra (10), we computed the mole percents of oxidized carbon in the solution and found 68% Embedded Image and 32% Embedded Image [no CO2(aq)-related Raman peaks were reported in Frantz’s study]. Therefore, PBE appears to overestimate the amount of Embedded Image by ~20% relative to the value obtained from Frantz’s spectroscopic data (10) and the γ ratio proposed by Schmidt (11); however, the overestimate is not as serious as for sulfuric acid aqueous solutions under ambient conditions, for example (26). This finding is consistent with our previous investigation of water under pressure, where we found that PBE performs better under extreme conditions than under ambient conditions in terms of equation of state of water, as well as in predicting static and electronic dielectric constants (20, 29). On the other hand, our PBE0 simulations yield results much closer to the experimental data. Our validation strategy led us to conclude that the use of the PBE functional is justified at high P and T, where some additional PBE0 simulations were carried out to double-check our PBE results.

Absence of CO2(aq)

The simulations described above led us to question the assumption of CO2(aq) being the major carbon species in supercritical aqueous fluids in deep Earth. We further increased the pressure to ~11 GPa and the temperature to 1000 K, that is, pressures similar to those at the bottom of Earth’s upper mantle, to investigate whether our conclusions on the absence of CO2(aq) in water did hold under more extreme conditions. The density of water at this P-T condition is 1.57 g/cm3 (20). The equation of state of water in this regime has been investigated in our previous study with the PBE functional (20).

Figure 2B shows the mole percents of dissolved carbon species as functions of simulation time obtained in first-principles MD simulations with the PBE and PBE0 functionals. We found 41% Embedded Image, 58% Embedded Image, and 1% H2CO3(aq) when using PBE; PBE0 simulations gave 52% Embedded Image and 48% Embedded Image, as shown in Table 1. Similar to our finding at 0.2 GPa and 823 K, PBE yields a higher concentration of Embedded Image than PBE0, but the results of the two simulations are closer than at lower P and T. Less than 1% CO2(aq) was found in both PBE and PBE0 simulations.

Because of the lack of experimental data, many studies under the upper mantle conditions used geochemical models based on extrapolations of available data. Using extrapolated thermodynamic data, Caciagli and Manning considered CO2(aq) as the dominant dissolved carbon species in aqueous calcite solutions up to 1.6 GPa and 1023 K (30). On the basis of geochemical models at lower P and T, Manning et al. inferred that bicarbonate ions would become less favored relative to CO2(aq) under the mantle conditions (12). However, Facq et al. performed Raman studies of CaCO3 solutions up to 7 GPa and 673 K and found that at P > 4 GPa, Embedded Image is the dominant carbon species (31). Furthermore, on the basis of theoretical models, Facq et al. speculated that at higher temperature, CO2(aq) may become a dominant species, depending on the pressure, though no exact temperature or pressure was reported (31). In the studies by Caciagli and Manning (30) and Facq et al. (31), using solutions at equilibrium with CaCO3 crystals, the concentrations of dissolved carbon were lower than those studied in this work; it is therefore challenging to experimentally detect the amount of CO2(aq) and definitely rule out the existence of molecular CO2. Our first-principles MD simulations clearly showed that up to 11 GPa and 1000 K, few Embedded Image ions are converted to CO2(aq).

To further verify whether CO2(aq) may be stable in water under extreme conditions, we directly dissolved CO2 into water under the same P-T condition of the simulation previously described (the carbonate ion was replaced by one CO2 molecule and one water molecule in the same simulation box at ~11 GPa and 1000 K) and found again that only a low percentage of CO2(aq) remained after the solution was equilibrated. Bicarbonate ions accounted for approximately 80% of dissolved carbon; this result was obtained using both xc functionals, as shown in Fig. 2C and Table 1.

We also found that H2CO3(aq) was present in the solution with a mole percent between 11 and 23%, depending on the functionals; note that this concentration is about 20 times larger than that of CO2(aq), as shown in Table 1. The formation of H2CO3(aq) via CO2(aq) hydration may follow one of two pathways (32): (i) CO2(aq) directly reacts with a water molecule to generate H2CO3(aq) or (ii) an intermediate Embedded Image is formed before H2CO3(aq). At 11 GPa and 1000 K, in both PBE and PBE0 simulations, we found that Embedded Image formed first, followed by H2CO3(aq); however, frequent proton hopping events were detected between H2CO3(aq) and water. The lifetime of H2CO3(aq) is on the subpicosecond scale. The existence of H2CO3 was also found in a 55.5 m CO2 water solution at ~32 GPa and 2000 K using first-principles MD with the GGA functional, but no concentration was given because the simulation was only 10 ps long (33).

In aqueous solutions under ambient conditions, H2CO3(aq) exists only in a very low concentration and rapidly reverts to CO2(aq) and water (12)Embedded Image(1)

Because it is difficult to detect H2CO3(aq) at ambient conditions, current geochemical models cannot distinguish between H2CO3(aq) and CO2(aq), and H2CO3(aq) and CO2(aq) are modeled as the same species (12). However, our findings point out that this conventional treatment may ignore the important presence of H2CO3(aq) as a neutral carbon species in deep Earth.

Under ambient conditions, the dissociation reaction of CO2(aq) in pure waterEmbedded Image(2)is rather slow, with an activation free energy of 21.7 kcal/mol (0.94 eV) (32, 34); more than 99% of CO2 molecules remain intact when the concentration of CO2(aq) is above 0.1 m. However, when the hydroxide ion (OH) is present, the activation free energy of the dissociation of CO2(aq) is reduced to 12.0 kcal/mol (0.52 eV). As a result, the dissociation reaction is accelerated markedly; it is about 107 faster than that in pure water (34). Under HP-HT conditions, water molecules are more easily dissociated, producing hydroxide ions, which then react with hydrated CO2(aq), thus enhancing its dissociation. The enhanced bending motion of C–O bonds at high P-T also facilitates the change of carbon hybridization from sp to sp2. In our simulations, the exact concentrations of charged ions depend on the xc functionals used, but both the GGA and hybrid functionals give the same result on CO2(aq) concentration: The latter is negligible in aqueous solutions at 11 GPa and 1000 K. Hence, we conclude that our calculations provide strong and robust evidence of the absence of CO2(aq) in water under this P-T condition.

Ion pairing

Another important question regarding the aqueous solutions studied here concerns the interactions between (bi)carbonate ions and metal ions. In Fig. 3, we plot the radial distribution functions (RDFs) and distributions of the C–Na distances obtained from our MD trajectories using the PBE and PBE0 functionals. We first discuss the results given by the two functionals, and we then extract robust features of both simulations.

Fig. 3 The solvation structures of carbon species and ion pairing distances in the 0.9 m Na2CO3 solutions.

(A) RDFs of carbon atoms (C) versus oxygen atoms of water (Ow). (B) Probability distributions of the distances between carbon atoms and sodium ions. The MD simulations on the 0.9 m Na2CO3 solutions were conducted at 0.2 GPa, 823 K and 11 GPa, 1000 K. Two xc functionals were compared: PBE and PBE0.

We examined the solution structures of (bi)carbonate ions in the solutions. In Fig. 3A, the C–Ow RDFs (where Ow denotes an O belonging to water molecules) obtained by PBE and PBE0 are very similar at 0.2 GPa, 823 K or 11 GPa, 1000 K. At 11 GPa and 1000 K, the first local minimum of the C–Ow RDFs is at ~4.54 Å, corresponding to the separation of the first and second solvation shells of (bi)carbonate ions. Although this minimum is not easily defined at 0.2 GPa and 823 K (because the first peak of the C–Ow RDFs also appears to vanish at ~4.54 Å), we used the same radial cutoff of 4.54 Å for the first solvation shells of (bi)carbonate ions under two P-T conditions to have a consistent comparison.

The ion pairing strength in an aqueous solution greatly depends on the static dielectric constant of water, ε0: the electrostatic interaction between ions in water is Embedded Image, where q1 and q2 are the electrical charges of the ions, and r is the distance between them. We note that PBE appears to overestimate ε0 under ambient conditions, whereas PBE0 yields a result in better agreement with experiment (35); hence, at the PBE level of theory, the overestimate of electrostatic screening leads to weaker ion pairing, as shown by lower PBE peaks at ~3.2 Å in Fig. 3B. Interestingly, under upper mantle pressure conditions, the discrepancy between ε0 evaluated within PBE and the experimental value is rather small (20), and we expect PBE to yield more realistic ion associations under those conditions than at ambient conditions.

At 0.2 GPa and 823 K, both functionals show marked peaks in the distribution of the C–Na distance at ~3.2 Å, indicating strong ion pairing between Na+ and Embedded Image or Embedded Image. The first peak of the C–Ow RDF is at ~3.5 Å, so the Na+-Embedded Image/Embedded Image pairs tend to reside inside the first solvation shell of (bi)carbonate ions. With increasing P and T to 11 GPa and 1000 K, respectively, the main peak of the distribution function of the C–Na distance decreases in intensity, indicating that ion pairing becomes weaker in denser water.

Table 2 shows the mole percents of the ion pairs in the first solvation shells of (bi)carbonate ions. As shown in Tables 1 and 2, we found that more than 99% of the Embedded Image anions were accompanied by one or two Na+ cations in their first solvation shells at 0.2 GPa and 823 K, whereas at 11 GPa and 1000 K, less than 63% of the Embedded Image ions had counterions in the first solvation shells. This may be attributed to the change of ε0: Our previous study showed that with increasing pressure, ε0 increases along an isotherm (20), and thus, the dielectric screening is enhanced and the electrostatic attraction between Na+ and Embedded Image or Embedded Image is weakened. The large dielectric screening also enhances the autoionization of water, producing a higher concentration of OH ions, whose presence accelerates the dissociation of CO2.

Table 2 The ion pairing in the first solvation shells of the carbon species in the 0.9 m Na2CO3 solutions.

The concentrations are shown as mole percents of the total dissolved carbon.

View this table:

The pairing of carbonate ions with alkali cations at 0.2 GPa and 823 K was also reported in a 1 m K2CO3 solution by using geochemical models (10), but no spectroscopic evidence was available because of the detection limit (10). Recently, Schmidt reported the Raman spectra, which were interpreted as showing Na+-Embedded Image pairing, for a 1.6 m Na2CO3 solution interfaced with quartz at T above 673 K. In the CaCO3-water-NaCl solutions, Facq et al. interpreted the measured Raman data by assuming the existence of the Embedded Image complex (36). In our study, despite the different ion pairing strengths given by the two functionals, we found that ion pairing between Na+ and Embedded Image or Embedded Image occurs at 0.2 GPa and 823 K, with its strength decreasing at 11 GPa and 1000 K.

Kinetics of reactions involving Embedded Image and Embedded Image

Bicarbonate and carbonate ions are the major species in the solutions studied here, with frequent proton transfer over picosecond time scales. Their chemical balance in the Na2CO3 solution isEmbedded Image(3)

The kinetics of the reaction is key to understanding the chemistry of oxidized carbon in supercritical water.

Figure 4 shows the distribution of protons hopping between a carbonate ion and its nearest water molecule in the Na2CO3 solution at 11 GPa and 1000 K. From Fig. 4, we calculated the effective free energy barrier of proton transfer to the right in reaction (3) as Embedded Image, where Pb and Pr are the probabilities of being on the top of the barrier and in the reactant state, respectively, and kB is the Boltzmann constant. At the PBE0 level of theory, the free energy barrier is estimated to be 4.6 kcal/mol (0.20 eV), whereas at the PBE level, it decreases to ~2.8 kcal/mol (0.12 eV). According to the Arrhenius equation, the rate constant is Embedded Image, where Ea is the activation energy. If we assume that PBE and PBE0 yield a similar prefactor A and approximate the activation energy by the free energy barrier, the rate constant of forward reaction (3) using PBE is about three times as large as that found using PBE0. This difference may be explained by the difference in O–H bond strengths obtained using the two functionals. We computed the vibrational density of states of the solution at 11 GPa and 1000 K in fig. S1 and found that the stretching frequency of O–H bonds given by PBE is about 150 cm−1 lower than that obtained with PBE0. As a result, at the PBE level of theory, the O–H bond is weaker and easier to break than at the PBE0 level, and the proton transfer rate is enhanced. It has been reported that at ambient conditions, the vibrational frequencies of the OH stretching mode are also underestimated using PBE compared with experimental data, whereas PBE0 yields results much closer to experiment (37). The underestimate decreases with increasing temperature (37). These findings suggest that it is necessary to use hybrid functionals to get accurate reaction rates in aqueous solutions, especially under ambient conditions. Under extreme conditions, the error of PBE becomes smaller, which is again consistent with the results of our previous studies (20, 29).

Fig. 4 Probability distributions of positions of protons hopping between Embedded Image and H2O in the Na2CO3 solution at 11 GPa and 1000 K.

The unit is Å−2. The reaction coordinate RO−O is the distance between the two neighboring oxygen atoms, Oc and Ow, in carbonate ions and water molecules, respectively, and δ is the proton displacement ROc−HRH−Ow. Two xc functionals were compared: PBE and PBE0.

At T = 1000 K, the thermal energy kBT is 2.0 kcal/mol (0.086 eV), comparable to the free energy barriers given by PBE and PBE0. Because of the low free energy barrier and enhanced thermal effects, both functionals yield a rate constant of forward reaction (3) about six orders of magnitude faster than that at ambient conditions [the rate constant of forward reaction (3) is 3.06 × 105 s−1 in seawater under ambient conditions (38)], indicating highly active ionic interactions in our solutions at 11 GPa and 1000 K.

It is worth noting that Fig. 4 shows only part of reaction (3); the proton transfer not only is limited to the distance between a carbonate ion and a neighboring water molecule but also may involve multiple water molecules connected by a hydrogen bond (HB) wire, via a Grotthuss mechanism (39, 40). In our solutions, protons move back and forth along HB wires over short time scales (<0.1 ps), leading to short-lived (bi)carbonate ions. However, if the thermal fluctuation breaks HB wires, protons cannot hop back and then bicarbonate ions persist for a few picoseconds. The forming and breaking of the HB wires are critical to the lifetime of carbon species in the calculations reported here.

In Raman experiments, the symmetric stretching mode of C–O/C–OH at ~1000 cm−1 is used to detect carbonate or bicarbonate ions. The vibrational period of this mode is 33 fs. For the carbon species, whose lifetime is comparable to 33 fs or even shorter, the Raman peaks are broadened considerably and may be difficult to identify experimentally.

In our MD simulations, the forces acting on nuclei were computed quantum mechanically, but protons were still treated as classical particles. If quantum nuclear effects are taken into account, we expect the proton transfer rate to be enhanced, though at such high temperatures, quantum effects may not be as important as under ambient conditions (41).

With both PBE and PBE0 functionals, van der Waals interactions are described inaccurately; however, in our study, we focused on the breaking and forming of covalent bonds, and the lack of dispersion forces is not expected to affect our main conclusions. In addition, we note that the results for the equilibrium volume and dielectric constant of ice obtained with PBE and van der Waals functionals agree much better at high pressure (for example, for ice VIII) than at low pressure (for example, for ice Ih and Ic) (42, 43), indicating that the lack of dispersion forces may affect the results for water and ice more under ambient conditions than under pressure.

Widely used geochemical models based on the Born function, for example, the Helgeson-Kirkham-Flowers model as implemented in SUPCRT92 (17), consider water as a continuum medium with no molecular structure taken into account, when calculating thermodynamic properties of electrolytes (44). In our study, we found that the microscopic properties of water at the molecular scale, for example, proton transfer mechanisms, play a substantial role in the chemical reactions of oxidized carbon, particularly under the high P-T conditions studied here.


In conclusion, we carried out first-principles MD simulations to investigate dissolved carbonate and CO2 in water at high pressure and high temperature, up to the conditions of Earth’s upper mantle. Contrary to popular geochemical models assuming that CO2(aq) is the major carbon species present in water, we found that most of the dissolved carbon at 11 GPa and 1000 K is in the form of solvated Embedded Image and Embedded Image anions. The two anions exchange protons with water on a picosecond time scale. The proton transfer along hydrogen bond wires driven by thermal fluctuations is responsible for the fast dynamics observed in these solutions. Under these extreme conditions, even H2CO3(aq) is more abundant than CO2(aq), a marked difference with respect to ambient conditions. Although it is well known that the autoionization of water is greatly enhanced under extreme conditions, for example, under the P-T conditions of Earth’s mantle, and that the pH of neutral water may be well below 7 (6), it has been long unclear how the notable properties of water under pressure affect the solvation of oxidized carbon. Using atomistic ab initio simulations, we showed that the coordination number of carbon changes in water under Earth’s upper mantle conditions. We also found that ion pairing between Na+ and Embedded Image or Embedded Image is greatly affected by P-T conditions, decreasing with pressure at 800 to 1000 K. Our results suggest that under extreme conditions, water transports carbon mostly through highly active ions, not CO2(aq) molecules.

The results reported here give important insights into deep carbon science. Our calculations predicted that the carbon dissolved in water-rich geofluids is in ionic forms in the bottom of Earth’s upper mantle, indicating that CO2 degassing may mostly occur close to Earth’s surface.


First-principles MD simulations were performed in the Born-Oppenheimer approximation using the Qbox code ( (45). Two xc functionals were used: the semilocal functional PBE and the hybrid functional PBE0. To improve the efficiency of our calculations using PBE0, the recursive subspace bisection method was used with a threshold of 0.02 (46, 47). The Brillouin zone of the supercell was sampled with the Γ point only. We used norm-conserving pseudopotentials (Pseudopotential Table, (48, 49) with a kinetic cutoff of 85 Ry, which was increased to 220 Ry for pressure calculations. The MD time step was 0.24 fs. Deuterium was used instead of hydrogen, so as to use a larger time step for computational convenience. Note that the density given in the main text was computed for light water. The temperature was maintained constant by the Bussi-Donadio-Parrinello thermostat (τ = 24.2 fs) (50).

To determine whether CO2 or Embedded Image was present in our simulations, we searched the atomic trajectories generated in our simulations for the three closest O atoms to a given C atom and ordered these O atoms according to their respective C–O distances. If the C–O distance of the third closest O atom was larger than that of the second one by more than 0.4 Å, the species was considered to be CO2; otherwise, it was defined as a Embedded Image anion. The C–O bond length in all sets of our MD runs was found to be 1.3 ± 0.2 Å. For H atoms around the carbonate ion, we sought the nearest O atom to each H atom, and if this O atom belonged to the carbonate ion, we concluded that a bicarbonate ion, Embedded Image, was formed. When two H atoms were bonded to Embedded Image simultaneously, the solute then became carbonic acid.

Using the criteria described above, we determined the type of solute molecules present in our MD simulations at each time step. Note that under HP-HT conditions, covalent bonds are frequently broken and reformed; hence, to establish how the concentration of the solute varied along our MD trajectories, we calculated the mole percents of carbon species as functions of simulation time. The mole percent of the ith species at time t was obtained asEmbedded Image(4)where ni(t) is the number of the snapshots containing the ith species between the time (t − τw) and t, and Nw is the total number of snapshots in this time interval. Here, τw was set to 1 ps.


Supplementary material for this article is available at

fig. S1. The vibrational density of states of the 0.9 m Na2CO3 solution at ~11 GPa and 1000 K.

fig. S2. Probability distributions of positions of protons hopping between Embedded Image and H2O in the Na2CO3 solution at 0.2 GPa and 823 K.

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.


Acknowledgments: We thank F. Giberti and D. Gygi for many useful discussions and D. A. Sverjensky for the critical reading of the manuscript and many discussions. Funding: This work was supported by the Alfred P. Sloan Foundation through the Deep Carbon Observatory (DCO) and by U.S. Department of Energy (DOE)–Basic Energy Sciences grant no. DE-SC0008938. Part of this work was carried out using the DCO Computer Cluster and resources of the Argonne Leadership Computing Facility (ALCF), provided by an award of computer time by the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program. ALCF is a DOE Office of Science User Facility supported under contract DE-AC02-06CH11357. We also acknowledge the University of Chicago Research Computing Center for support of this work. Author contributions: D.P. and G.G. designed the research. Calculations were performed by D.P. All authors contributed to the analysis and discussion of the data and the writing of 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 and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.
View Abstract

Stay Connected to Science Advances

Navigate This Article