Large presence of carbonic acid in CO$_2$-rich aqueous fluids under Earth's mantle conditions

TL;DR

This study reveals that carbonic acid is a major carbon species in Earth's deep mantle conditions, significantly impacting the deep carbon cycle and water-rock interactions, contrary to previous assumptions at ambient conditions.

Contribution

The paper demonstrates through ab initio simulations that carbonic acid is abundant in high-pressure, high-temperature aqueous fluids, highlighting its role in Earth's deep carbon cycle.

Findings

Carbonic acid becomes the dominant carbon species at ~10 GPa and 1000 K.

Proton transfer between carbonic acid and bicarbonate ions enhances solution conductivity.

Carbonic acid's buffering capacity influences water-rock interactions in Earth's interior.

Abstract

The chemistry of carbon in aqueous fluids at extreme pressure and temperature conditions is of great importance to Earth's deep carbon cycle, which substantially affects the carbon budget at Earth's surface and global climate change. At ambient conditions, the concentration of carbonic acid in water is negligible, so aqueous carbonic acid was simply ignored in previous geochemical models. However, by applying extensive ab initio molecular dynamics simulations at pressure and temperature conditions similar to those in Earth's upper mantle, we found that carbonic acid can be the most abundant carbon species in aqueous CO$_2$ solutions at ~10 GPa and 1000 K. The mole percent of carbonic acid in total dissolved carbon species increases with increasing pressure along an isotherm, while its mole percent decreases with increasing temperature along an isobar. In CO$_2$-rich solutions, we found significant proton transfer between carbonic acid molecules and bicarbonate ions, which may enhance the conductivity of the solutions. The effects of pH buffering by carbonic acid may play an important role in water-rock interactions in Earth's interior. Our findings suggest that carbonic acid is an important carbon carrier in the deep carbon cycle.

Equations

$$\mathrm{C}{\mathrm{O}}_2(\mathrm{aq})+\mathrm{O}{\mathrm{H}}^{-}$$

$$\mathrm{C}{\mathrm{O}}_2(\mathrm{aq})+{\mathrm{H}}_2\mathrm{O}$$

$$\mathrm{HC}{\mathrm{O}}_3^{-}+{\mathrm{H}}^{+}\rightleftharpoons {\mathrm{H}}_2\mathrm{C}{\mathrm{O}}_3(\mathrm{aq})$$

$${\mathrm{H}}_3\mathrm{C}{\mathrm{O}}_3^{+}\rightleftharpoons {\mathrm{H}}_2\mathrm{C}{\mathrm{O}}_3(\mathrm{aq})+{\mathrm{H}}^{+}$$

Full text

I Introduction

The global carbon cycle is of great importance to Earth’s climate, human’s energy consumption, and sustainable development. Although a tremendous number of studies have focused on the carbon cycle in the atmosphere, oceans, and the shallow crust, less is known about the carbon cycle in Earth’s deep interior hazen2013deep , which may host more than 90% of Earth’s carbon Falkowski2000 and also actively exchanges carbon with Earth’s near-surface reservoirs through volcanism and subduction Manning2014 ; Kelemen2015 . In the deep carbon cycle, fluids containing water as a primary component play a key role in transporting carbon in the deep mantle and back to magmas in volcanoes and eventually to the atmosphere as carbon dioxide (CO2) manning2004chemistry ; Manning2014 . However, many basic questions about carbon-bearing aqueous fluids are poorly known; for example, what are the aqueous carbon species in the deep carbon cycle? A major obstacle to understanding deep carbon transport is the lack of knowledge of carbon reactions in water at the extreme conditions found in Earth’s deep interior manning2013chemistry .

Both pressure (P) and temperature (T) increase with increasing depth inside Earth. For example, in the upper mantle, pressure can reach $\sim$13 GPa and temperature $\sim$1700 K Frost2008 . Most crust and mantle fluids are at vapor-liquid supercritical conditions, and their properties differ fundamentally from those at ambient or near-ambient conditions, which has been shown by many recent studies manning2018 . For oxidized carbon dissolved in water, many geochemical models assume that CO2(aq) is the major carbon species, e.g., Refs Mader1991 ; Duan2006 ; Zhang2009 ; Holland2011 , but recent studies show that carbon-containing ions, HCO${}_{3}^{-}$ and CO${}_{3}^{2-}$, may be more abundant than CO2(aq) in Earth’s mantle Facq2014 ; Pan2016 .

Another possible oxidized carbon species is carbonic acid (H2CO3) adamczyk2009real ; loerting2010aqueous . At ambient conditions, dry H2CO3 is very stable with a half-life as long as 180,000 years Loerting2000 , whereas in water it decomposes rapidly and its concentration is only about 0.1% of CO2(aq) concentration manning2013chemistry ; thus, in geochemical modeling H2CO3(aq) has simply been treated as CO2(aq), except in a very recent model huang2019extended . At 2.4 GPa and 553 K, Wang et al. found possible spectroscopic evidence of H2CO3(aq), suggesting it might be an important species in aqueous CO2 solutions under extreme conditions Wang2016 , whereas Abramson et al. commented that the new vibrational band could be also attributed to HCO${}_{3}^{-}$ at high concentration Abramson2017a . Note that these experimental P-T conditions can only be found in the shallow or cold mantle areas. It is very challenging to experimentally detect carbon species at both high pressures (HP) and high temperatures (HT), as found in the deep mantle. In our previous theoretical work, we applied ab initio molecular dynamics simulations to study CO2 in water at $\sim$11 GPa and 1000 K, and found that about 80% of dissolved carbon is HCO${}_{3}^{-}$ and 10$\sim$20% H2CO3(aq) as a minor product Pan2016 . The solution is very dilute with x(CO2) = 0.016. In fact, carbon-bearing fluids in deep Earth may range from water-rich to CO2-rich andersen2001fluid . The carbon speciation in CO2-rich aqueous fluids under Earth’s mantle conditions is still unknown.

Here, by applying extensive ab initio molecular dynamics (AIMD) simulations, we studied the CO2-H2O mixtures at various concentrations from 3 to 11 GPa and between 1000 and 1400 K, P-T conditions similar to those in subduction zones in Earth’s upper mantle. Our total simulation time exceeds 3.8 nanoseconds. We found that at $\sim$10 GPa and 1000 K, H2CO3(aq) can be the predominant carbon species in the aqueous CO2 solutions. With increasing pressure along an isotherm, the mole percent of H2CO3(aq) in total dissolved carbon species increases, while with increasing temperature along an isobar, its mole percent decreases. We show the possible P-T range where H2CO3(aq) may exist. The rich chemistry of H2CO3(aq) may play an important role in the dissolution of carbonate minerals in aqueous fluids in deep Earth. Our findings suggest that H2CO3(aq) is an important carbon carrier in the deep carbon cycle.

II Results and discussion

II.1 Equation of state of CO2-H2O mixtures

First, we calculated the equation of state (EOS) of CO2-H2O mixtures given by density functional theory (DFT) with the PBE exchange-correlation functional. In previous studies we used the same level of theory to calculate many properties of water under HP-HT conditions, such as equation of state Pan2013 , dielectric constant Pan2014 , Raman and IR spectra, and ionic conductivity Rozsa2018 , as well as the reactions of CO2 and carbonates in HP-HT water Pan2016 . In Table SI, we compare the DFT pressures of CO2-H2O systems with those obtained by popular geochemical models Duan2006 ; Zhang2009 and classical force-field MD simulations at various CO2 concentrations. The force field is from Duan and Zhang’s work in 2006 Duan2006 . Since the EOS in Ref. Duan2006 was obtained by fitting the MD simulation data using the force field and the experimental data, it is not surprising the pressures calculated by Duan and Zhang’s EOS in Ref. Duan2006 are in excellent agreement with those calculated by our MD simulations with the same force field. The EOS pressures of Ref. Duan2006 are generally smaller than the results obtained by Zhang and Duan’s EOS Zhang2009 in 2009, which was designed for C-O-H fluids in Earth’s mantle. It was recently reported Zhang2016 that for CO2-H2O binary systems, the EOS in Ref. Duan2006 performs better than the one in Ref. Zhang2009 .

In our simulations, we varied the initial mole fraction of CO2(aq), which is defined as $\mathrm{x(CO_{2})}=\frac{n\mathrm{(CO_{2})}}{n\mathrm{(CO_{2})}+m\mathrm{(H_{2}O)}}$, where $n\mathrm{(CO_{2})}$ and $m\mathrm{(H_{2}O)}$ are the numbers of CO2 and water molecules, respectively. When x(CO2) is smaller than 0.1, our DFT pressures are slightly larger than those obtained by the EOS in Ref. Duan2006 , while for x(CO2) larger than 0.1, the DFT pressures become smaller. Note that in the simulations using the force field from Ref. Duan2006 , no chemical reaction occurs between CO2 and water, whereas in the ab initio simulations the carbon atoms may be converted from sp to sp2 hybridization, and vice versa. The newly generated carbon species may have different molecular volumes, which affect the EOS of the mixtures. Overall, the agreement between our ab initio data and available EOS data is satisfactory, indicating that the PBE functional is not only suitable for water under extreme conditions, but may also be used to study CO2-H2O mixtures at HP and HT.

II.2 Speciation of CO2-H2O mixtures

Fig. 1(A) shows the relative concentrations of carbon species as functions of x(CO2) obtained from AIMD simulations at $\sim$10 GPa and 1000 K. The mole percents of CO2(aq), HCO${}_{3}^{-}$, CO${}_{3}^{2}-$, and H2CO3(aq) as functions of simulation time are shown in Fig. S1. With increasing x(CO2), the mole percent of HCO${}_{3}^{-}$ in total dissolved carbon species decreases, and more CO2(aq) molecules were found in the simulations. The crossover between HCO${}_{3}^{-}$ and CO2(aq) happens when x(CO2) is about 0.25. We found few carbonate ions, and the mole percent of CO${}_{3}^{2-}$ is less than 10% in all the simulations in Fig. 1. Interestingly, when x(CO2) is in the range of 0.14 to 0.31, carbonic acid becomes the predominant carbon species, more abundant than any other carbon-containing molecules or ions. When increasing x(CO2) up to $\sim$0.20, the mole percent of H2CO3(aq) reaches its maximum ($\sim$47%); at higher x(CO2), the mole percent of H2CO3 slowly decreases. During the decrease, H2CO3(aq) is always more abundant than HCO${}_{3}^{-}$. After increasing the temperature from 1000 K to 1400 K at $\sim$10 GPa, the HCO${}_{3}^{-}$-CO2 crossover happens at a lower carbon concentration with x(CO2) = 0.10 (Fig. 1(B)). At 1400 K, the mole percent of H2CO3(aq) is about 20% when x(CO2) is between 0.03 and 0.33, and does not have a clear maximum value.

II.3 Molecular geometry of H2CO3

The carbonic acid molecule has one carbon atom and three oxygen atoms, which form the trigonal planar shape Reisenauer2014 ; bucher2014clarifying . Two hydrogen atoms, which are each bonded to an oxygen atom, appear not always in the plane of the carbon and oxygen atoms in the MD simulations. Fig. 2 shows the distribution of the two O=C-O-H dihedral angles. The H2CO3 conformer with the lowest free energy is the cis-cis conformer, where the two dihedral angles are both 0∘. The conformer with the second lowest free energy is cis-trans, where the two dihedral angles are 0∘ and 180∘. The trans-trans conformer, where the two dihedral angles are both 180∘, has much higher free energy than the former two. We calculated the free energy difference between the $i$th and $j$th conformer by $\Delta F=-k_{B}T\ln\left(\frac{P_{i}}{P_{j}}\right)$, where $k_{B}$ is the Boltzmann constant and $P_{i}$ is the probability of the $i$th conformer. The free energies are calculated with respect to the cis-cis conformer. At $\sim$10 GPa and 1000 K, the free energy of the cis-trans conformer is 0.14 kcal/mol (0.006 eV), while the trans-trans conformer has the free energy of 1.82 kcal/mol (0.079 eV). When increasing the temperature from 1000 to 1400 K along the isobar of $\sim$10 GPa, the free energy of the cis-trans conformer increases to 0.19 kcal/mol (0.008 eV), while the trans-trans one increases to 2.31 kcal/mol (0.100 eV).

In the gas phase at 0 K, the free energies of the cis-trans and trans-trans conformers are 1.6 and 10.1 kcal/mol respectively Reisenauer2014 , much higher than the respective free energies in the aqueous solutions studied here. This may be attributed to the dielectric screening of water. The three H2CO3 conformers have different electrostatic interactions. The two hydrogen atoms with partial positive charges in the trans-trans conformer are closer than in the other two conformers, so the stronger electrostatic repulsive force causes a higher energy in the trans-trans conformer. Because the dielectric screening of water weakens electrostatic interactions, we see the lower free energies for the trans-trans and cis-trans conformers in water than in the gas phase. When increasing T along an isobar, the dielectric constant of water decreases Pan2013 , so the screening becomes weaker and the free energies become larger again.

II.4 Reaction mechanisms of H2CO3(aq)

In Fig. 3, we analyzed how H2CO3(aq) forms and dissociates at $\sim$10 GPa and 1000 K. In our AIMD simulations, $\sim$90% of the H2CO3(aq) formation reactions follow a two-step mechanism Stirling2010 : (i) the reaction of CO2(aq) with OH- or water to form HCO${}_{3}^{-}$,

[TABLE]

and (ii) the protonation of HCO${}_{3}^{-}$ to form H2CO3(aq):

[TABLE]

Because the solutions are acidic, it is very rare to see free OH-. In reaction (1), the OH- ion is often associated with the H3O+ ion generated by the self-ionization of water. If the proton keeps hopping between the OH- and H3O+ ions, reaction (1) will also go backward and forward frequently. When the proton in H3O+ transfers away from the ion pair along a hydrogen-bond wire via the Grotthuss mechanism Marx2006 , the chemical equilibrium of reaction (1) shifts to the right. Similarly, in reaction (2) where the CO2 molecule reacts with a water molecule, if the released proton is accepted by a third species and cannot transfer back, the generated HCO${}_{3}^{-}$ will not often change back to CO2(aq).

In reaction (3), the HCO${}_{3}^{-}$ ion accepts one proton to become H2CO3(aq). Fig. 3(A) shows that when the carbon concentration is low, most of the protons are from H3O+, whereas with increasing x(CO2), more and more protons come from H2CO3(aq), indicating that protons can transfer between HCO${}_{3}^{-}$ and H2CO3(aq) in CO2-rich solutions (see Fig. 4). We even found H3CO${}_{3}^{+}$ as an intermediate proton donor with lifetime less than 0.15 ps. In fact, about 10% of H2CO3(aq) is generated from the dissociation of H3CO${}_{3}^{+}$:

[TABLE]

The sp2 carbon species such as H2CO3, HCO${}_{3}^{-}$, and H3CO${}_{3}^{+}$ may form a network to conduct protons in water. It was expected that the ionic conductivity of CO2-H2O mixtures decreases significantly with increasing x(CO2), if only the self-ionization of water contributes to ionic conduction manning2018 . This expectation may not be valid, as the generated sp2 carbon species can enhance the ionic conductivity.

At $\sim$10 GPa and 1000 K, about 90% of H2CO3(aq) dissociates to become HCO${}_{3}^{-}$ via backward reaction (3). The lifetime of H2CO3(aq) is typically less than 1 ps (see Fig. S4). When the solutions are dilute, most of protons transfer to water molecules to form H3O+, whereas with increasing x(CO2), the role of HCO${}_{3}^{-}$ or H2CO3(aq) as proton acceptors becomes more important, as shown in Fig. 3(B). Interestingly, only 1$\sim$3% of HCO${}_{3}^{-}$ further decompose into CO2(aq) (see Fig. S5). We found two decomposition pathways: backward reaction (1) and (2). Decomposition via backward reaction (1) is more likely to happen when x(CO2) is low. With increasing x(CO2), decomposition via backward reaction (2) becomes more frequent. In backward reaction (1), if the OH- ion remains in the first solvation shell of CO2, the O2C-OH bond often keeps breaking and forming unless the OH- ion accepts a proton and becomes a H2O molecule. Averagely, CO2(aq) generated through backward reaction (2) has a longer lifetime than that from backward reaction (1). About 10% of H2CO3(aq) does not dissociate or decompose, but receives protons to become the intermediate H3CO${}_{3}^{+}$ through backward reaction (4). In the gas phase, the decomposition of H2CO3 occurs via the concerted mechanism that involves the dehydroxylation of one OH- group and the deprotonation of the other OH- group simultaneously Loerting2000 ; Kumar2007 . Here, we found only the stepwise process through backward reactions (3), and (2) or (1). Galib and Hanna studied the decomposition mechanism of H2CO3(aq) using water clusters with varied sizes at ambient conditions, and found that the hydrogen bonding environment around the H3O+ ion as found in bulk water leads to stepwise decomposition Galib2014 , which indicates that at the supercritical conditions studied here there are still strong water-ion bonding interactions.

The reaction mechanisms of H2CO3(aq) formation and dissociation help us understand why H2CO3(aq) can be the predominant carbon species in supercritical water at $\sim$10 GPa and 1000 K. When the CO2-H2O solution is dilute, most CO2 reacts and becomes HCO${}_{3}^{-}$ via reaction (1) or (2) Pan2016 . With increasing x(CO2), the solution becomes more acidic and HCO${}_{3}^{-}$ stays closer to H3O+, so reaction (3) is more likely to happen. Thus, when increasing x(CO2) up to $\sim$0.20, the concentration of H2CO3(aq) increases to its maximum. However, when further increasing x(CO2), a significant amount of CO2(aq) does not react and the mole percent of HCO${}_{3}^{-}$ in total dissolved carbon species keeps decreasing. As a result, the mole percent of H2CO3(aq) also decreases. When x(CO2) $>$ 0.14, H2CO3(aq) is more abundant than HCO${}_{3}^{-}$, indicating that H2CO3(aq) is energetically favored over HCO${}_{3}^{-}$ in CO2-rich solutions under extreme conditions.

II.5 P-T range of H2CO3(aq)

When increasing the temperature from 1000 K to 1400 K at $\sim$10 GPa, the mole percent of H2CO3(aq) becomes smaller. The concentration of CO2(aq) at thermal equilibrium at 1400 K is higher than that at 1000 K when the initial CO2 concentration is the same. As a result, the H2CO3 formation reactions are less likely to happen. Big molecules or ions tend to decompose at high temperatures, so at 1400 K, 9$\sim$10% of HCO${}_{3}^{-}$ decomposes into CO2(aq), almost one order of magnitude more than that at 1000 K.

Not only temperature but also pressure affects the presence of H2CO3(aq). We lowered the pressure from $\sim$10 GPa to $\sim$8.0 GPa, $\sim$6.1 GPa, and $\sim$3.8 GPa in the AIMD simulations at 1000 K and x(CO2) = 0.016 (see Fig. 1(C)). Because the solutions are dilute, H2CO3(aq) is not the predominant carbon species, but we can estimate the formation of H2CO3(aq) with increasing x(CO2) based on the reaction mechanisms discussed above. At $\sim$8.0 GPa, most of the CO2 molecules still convert into HCO${}_{3}^{-}$, but the concentration of HCO${}_{3}^{-}$ is lower than that at $\sim$10 GPa; thus, with increasing x(CO2), we also expect to have a significant amount of H2CO3(aq), despite a lower concentration than that at $\sim$10 GPa. Fig. 1(C) shows that when the pressure is lower than $\sim$6.3 GPa, CO2(aq) becomes predominant. Because CO2(aq) reacts slowly at low pressures, the transition pressure of $\sim$6.3 GPa should be an upper limit. Generally, with decreasing pressure at 1000 K, the mole percent of H2CO3(aq) decreases.

In Fig. 5, we show the possible P-T range for H2CO3(aq) as a significant solute in water based on our AIMD simulations and experimental findings. With increasing P along an isotherm, the concentration of H2CO3(aq) increases, while with increasing T along an isobar, less H2CO3(aq) forms. There can be a significant amount of H2CO3(aq) at P $\sim$10 GPa and T $\leq$ 1400 K. Wang et al. showed possible spectroscopic evidence of H2CO3(aq) at 2.4 GPa and 553 K in a mixture of CO2 and H2O Wang2016 . For aqueous CO2 solutions, Abramson et al. found that with increasing P from 4.1 GPa to 6.2 GPa at 523 K, the Raman signal of CO2(aq) gradually disappears and a new peak at $\sim$1040 cm*-1* might be attributed to high concentration HCO${}_{3}^{-}$ or H2CO3(aq) Abramson2017a . No evidence for the formation of H2CO3(aq) below 2.4 GPa is reported experimentally, so in Fig. 5 we use 2.4 GPa and 523 K as the P-T boundary of H2CO3(aq). In Fig. 5, we also show the P-T conditions of water in Earth’s upper mantle Thompson1992 and subducting slab surfaces Syracuse2010 , which largely overlap with the P-T range of the formation of H2CO3(aq). This indicates that there can be a significant amount of H2CO3(aq) in aqueous geofluids in subduction zones in the upper mantle. Dissolved CO2 in aqueous geofluids may be gradually converted into H2CO3(aq) and HCO${}_{3}^{-}$/CO${}_{3}^{2-}$ with increasing depth. The sp2-carbon-bearing fluids may also be brought back to Earth’s surface by volcanism and the CO2 degassing would occur in shallow areas.

In the DFT calculations, we used the semilocal exchange-correlation functional PBE Perdew1996 . It has been reported that generalized gradient approximations (GGAs) may not be able to describe doubly charged anions in water accurately at ambient conditions wan2014electronic . However, PBE performs better at HP-HT conditions than at ambient conditions, as shown in many of our previous studies Pan2013 ; Pan2014 ; Pan2016 ; Rozsa2018 . Our previous finding that CO2(aq) is absent at $\sim$11 GPa and 1000 K Pan2016 was further extended to lower pressures and temperatures by experiment Abramson2017a , indicating that our computational method is capable to predict the properties of aqueous carbon solutions under extreme conditions. To validate our results, we compared the PBE simulations with the simulations using the hybrid functional PBE0 adamo1999toward ; Pan2016 , which often gives results in better agreement with experimental ones due to the smaller charge delocalization error Cohen2008 . When we dissolved CO2 in water (x(CO2) = 0.016) at $\sim$11 GPa and 1000 K , we found that in both PBE and PBE0 simulations, HCO${}_{3}^{-}$ is the predominant carbon species and its mole percents are very similar: 79.8% vs. 75.0% Pan2016 . PBE0 predicts even more H2CO3(aq) than PBE. Thus, our finding that H2CO3(aq) can be the most abundant carbon species under extreme conditions would not be affected by the PBE error.

In addition, van der Waals interactions are not well captured by the PBE and PBE0 functionals. However, previous studies showed that PBE and van der Waals functionals give closer results at high pressure than at ambient pressure for the equilibrium volume and dielectric constant of ice murray2012dispersion , indicating that for water and ice under pressure, the dispersion interactions are not so important as those at ambient conditions. Besides, the breaking and forming of covalent bonds are little affected by weak van der Waals interactions, so the lack of dispersion interactions should not change our main results here.

The chemical reactions involving H2CO3(aq) greatly affect the acidity of aqueous carbon solutions, which may have an important impact on water-rock interactions in deep Earth manning2013thermodynamic . It has been estimated that about half of carbon in subducting slabs or even more is brought back to Earth’s surface Kelemen2015 . The previous geochemical studies on the devolatilization in subduction zones considered only CO2(aq) as the major carbon species in water Manning2014 . If so, the amounts of carbon removed from subducting slabs are too small to account for the large volumes of carbon found in magmas Manning2014 . It was recently reported that significant quantities of carbonate minerals are dissolved by the infiltrating fluids, and are brought back to magmas, which may potentially balance CO2 emissions from volcanoes Ague2014 . However, the chemical mechanism of this dissolution process is not yet known. At HP-HT conditions, the solubilities of carbonate minerals in pure water increase, but are still not very large Facq2014 ; Pan2013 . If we consider only molecular CO2 in the solutions, the solubilities should become even lower manning2013chemistry . Our current work suggests that the effects of pH buffering by H2CO3(aq) may play an important role in enhancing the dissolution of carbonate minerals manning2013thermodynamic .

III Conclusion

In conclusion, we conducted ab initio molecular dynamics simulations to study dissolved CO2 in both water-rich and CO2-rich aqueous solutions from 3 to 11 GPa and at 1000 and 1400 K, P-T conditions similar to those found in subduction zones in Earth’s upper mantle. We found that H2CO3(aq), whose concentration is negligible in water at ambient conditions, is the most abundant carbon species in the aqueous CO2 solutions at $\sim$10 GPa and 1000 K, when the initial mole fraction of CO2(aq) is between 0.14 and 0.31. With increasing pressure along an isotherm, the mole percent of H2CO3(aq) in total dissolved carbons species increases, while with increasing temperature along an isobar, its mole percent decreases. Contrary to popular geochemical models assuming carbon mainly exists in the form of molecular CO2 in aqueous geofluids, our study suggests that big carbon-containing molecules like H2CO3(aq) cannot be ignored under deep Earth conditions. We found significant proton transfer between H2CO3(aq) and HCO${}_{3}^{-}$ in CO2-rich solutions, which may enhance the ionic conductivity of the fluids. The rich chemistry of H2CO3(aq) may substantially affect water-rock interactions in deep Earth. For example, the effects of pH buffering by H2CO3(aq) may play an important role in enhancing the dissolution of carbonate minerals. Our findings suggest that H2CO3(aq) is an important carbon carrier in the deep carbon cycle. Not only inside Earth, it may also largely exist in other water-rich planets Marounina2019 .

IV Methods

IV.1 Ab initio molecular dynamics (AIMD)

AIMD simulations were carried out with the Qbox code (version 1.63.8, http://qboxcode.org/) gygi2008 in the Born-Oppenheimer approximation. We used plane wave basis sets, the PBE exchange-correlation functional Perdew1996 , and norm-conserving pseudopotentials with a kinetic energy cutoff of 85 Ry (Pseudopotential Table, http://fpmd.ucdavis.edu/potentials) Hamann1979 ; Vanderbilt1985 . In the pressure calculations, the cutoff was increased to 145 Ry. We sampled the Brillouin zone of the simulation box at the $\Gamma$ point only. We used deuterium instead of hydrogen to have a larger time step (0.24 fs) for computational convenience. The temperature was controlled by the Bussi-Donadio-Parrinello thermostat ($\tau$ = 24.2 fs) Bussi2007 . The numbers of CO2 and water molecules in simulation boxes are listed in Table SI. In MD simulations, the equilibration time is 20 ps. Each production run is between 180 and 580 ps long, as shown in Fig. S1-S3.

We used the atomic trajectories obtained from AIMD simulations to determine the chemical form of carbon species. To determine whether a carbon atom is sp or sp2 hybridized, we sorted the three smallest C-O distances near each carbon atom. If the difference between the second and third C-O distances is more than 0.4 Å, the species was defined as CO2; otherwise, it was considered as a CO${}_{3}^{2-}$ anion. For each hydrogen atom near carbon species, we looked for the closest oxygen atom to form the O-H bond. The bicarbonate ion and carbonic acid have one and two O-H bonds, respectively.

IV.2 Classical molecular dynamics

Classical MD simulations were performed using simulation boxes containing 320$\sim$512 molecules. We used the Gromacs package version 5.1.2 Abraham2015 . The time step is 1 fs. The force fields for CO2 and water are from Ref. Duan2006 . The temperature was maintained by the Bussi-Donadio-Parrinello thermostat ($\tau$ = 0.1 ps). Electrostatic interactions were calculated using the fast smooth particle-mesh Ewald method with a grid spacing of 1 Å and cubic interpolation. A cutoff radius of 10 Å was used for van der Waals interactions. Production runs are $\sim$2 ns after 1 ns equilibration.

V Acknowledgements

We thank Giulia Galli and Viktor Rozsa for their helpful discussions. N.S. acknowledges the Hong Kong PhD Fellowship Scheme. D.P. acknowledges support from the Croucher Foundation through the Croucher Innovation Award, Hong Kong Research Grants Council (Projects ECS-26305017, GRF-16307618), National Natural Science Foundation of China (Project 11774072), and the Alfred P. Sloan Foundation through the Deep Carbon Observatory.

VI Contributions

D.P. designed the research. N.S. performed the ab initio and classical simulations. All authors contributed to the analysis and discussion of the data and the writing of the manuscript.