Introduction

Aqueous fluids play a critical role in transporting carbon between Earth’s surface and interior1,2,3, which is a substantial part of Earth’s carbon cycle, with great implications for global climate and human energy consumption. It has long been assumed that aqueous carbon solutions under extreme pressure (P) and temperature (T) conditions are made by mixtures of neutral gas molecules4, e.g., H2O, CO2, CH4; however, recent studies showed that important chemical reactions occur between water and carbon species, resulting in significant amounts of ionic products, which may further participate in water-rock interactions and the formation of diamonds in Earth’s interior5,6,7,8,9,10,11. Most of the previous studies focus on the properties of aqueous carbon solutions in the bulk phase. In fact, aqueous solutions in deep Earth are often confined to the nanoscale in pores, grain boundaries, and fractures of Earth’s materials12,13,14, where the physical and chemical properties of solutions may be dramatically different from those of bulk solutions. In addition, in carbon capture and sequestration efforts, CO2 mineralization occurring in water trapped in porous rocks offers an efficient and secure method to permanently store carbon underground with a low risk of return to the atmosphere15. The behavior of aqueous carbon solutions under nanoconfinement at extreme P-T conditions is of great importance to the deep carbon cycle and CO2 storage, but is poorly understood on the molecular scale.

Previous studies reported that nanoconfinement substantially affects properties of water, e.g., equation of state16,17,18, phase behavior19,20,21, dielectric constant22,23,24,25,26, and diffusion27,28,29; as a result, the reactivity of solutes under confinement may be very different from that in bulk solutions30. The dimensional reduction and increased fluid density could enhance reactions between small solutes in nanoconfinement31,32, whereas reactions involving large reactants or intermediates may be sterically hindered33. Further, the increase of the dielectric constant of nanoconfined water parallel to the confining surface leads to the stabilization of aqueous reaction products with charges33, causing the enhanced autodissociation of water23. The solid–liquid interface also greatly affects the properties of confined aqueous solutions34. Preferential adsorption of solutes at the confining interface may shift reaction equilibria. For example, in the production of methane from carbon dioxide at hydrothermal vent conditions (CO2 + 4 H2⇌ CH4 + 2 H2O), hydrophilic pore surfaces adsorb water, favoring the production of methane35.

Nanoconfinement and interface chemistry may both likely change the properties of aqueous carbon solutions, but a molecular understanding is lacking on how chemical speciation and reaction mechanisms are affected. It was experimentally found that magnesite precipitates much faster in nanoscale water films than in bulk water36. Because it is very challenging to study aqueous solutions under nanoconfinement in experiment, atomistic simulations are widely used. Many studies applied classical force fields27,29,34,37, which were usually designed for bulk solutions at ambient conditions; their accuracy at extreme conditions is not well tested. As a comparison, ab initio molecular dynamics (AIMD) simulations do not rely on experimental input or empirical parameters38,39,40. We solve the many-body electronic structure numerically, so the breaking and forming of chemical bonds, electronic polarizability, and charge transfer are all treated at the quantum mechanical level40,41. The AIMD method is widely considered as one of the most reliable methods to make predictions, and many simulation results were later confirmed by experiments40,41.

Here, we performed extensively long AIMD simulations to study CO2(aq) solutions nanoconfined by graphene and stishovite (SiO2) at 10 GPa and 1000 ~ 1400 K. These P-T conditions are typically found in Earth’s upper mantle. We compared the CO2(aq) reactions in nanoconfinement with those in the bulk solutions, and examined how weak and strong interactions between confining walls and confined solutions affect chemical speciation and reaction mechanisms. Although graphene is not found in deep Earth so far, it provides a good comparison with stishovite. In graphene confinement, there are no chemical reactions between graphene and solutions, whereas the dangling atoms in stishovite actively participate in aqueous carbon reactions, so we can compare the effects of spatial confinement with and without interface chemistry. What’s more, thanks to the rapid development in the fabrication and characterization of 2D materials in recent years, experimentalists are now able to delicately measure the properties of aqueous solutions under graphene nanoconfinement30, so we hope our study can also attract many follow-up experiments. Our work is relevant to the carbon transformation in deep Earth, and also helps us to understand atomistic mechanisms of CO2 mineralization in the carbon capture and storage.

Results and discussion

Graphene nanoconfinement

We first studied CO2(aq) solutions confined by two graphene sheets at ~ 10 GPa, and 1000 ~ 1400 K (Fig. 1a). The graphene sheet separation was 9.0 and 9.2 Å at 1000 and 1400 K, respectively. We modeled the graphene sheets using a distance-dependent potential acting on the carbon and oxygen atoms, which was fitted to the interaction energies calculated using diffusion quantum Monte Carlo42 and van der Waals density functional theory43 (see Supplementary Methods). We calculated the pressure of confined solutions parallel to the graphene sheets, which is ~10 GPa (see Supplementary Methods). In addition, we also used atom number density profiles to calculate actual volumes that aqueous carbon solutions occupy, and then applied the equation of state of CO2 and water mixtures to obtain the pressure44.

Fig. 1: Snapshots of ab initio molecular dynamics simulations under confinement.
figure 1

a CO2(aq) confined by graphene sheets. b CO2 confined by cleaved stishovite (SiO2) (100) slabs. Red balls represent oxygen atoms in solutions, and pink balls are oxygen atoms in SiO2. Gray, white, and yellow balls are carbon, hydrogen, and silicon atoms, respectively.

Full size image

We directly dissolved CO2 molecules in the supercritical water, and the initial mole fraction of CO2(aq) is 0.185. The CO2 molecules reacted frequently with water, and we performed long AIMD simulations until the concentrations of carbon species reached equilibria (see Supplementary Fig. 3). Initially, the reaction between CO2(aq) and H2O produces bicarbonate ions (HCO\({}_{3}^{-}\)):

$${{{{{{{{\rm{CO}}}}}}}}}_{2}({{{{{{{\rm{aq}}}}}}}})+2{{{{{{{{\rm{H}}}}}}}}}_{2}{{{{{{{\rm{O}}}}}}}}\rightleftharpoons {{{{{{{{\rm{HCO}}}}}}}}}_{3}^{-}+{{{{{{{{\rm{H}}}}}}}}}_{3}{{{{{{{{\rm{O}}}}}}}}}^{+}.$$
(1)

This reaction in some cases occurs in one step, or involves the dissociation of water so that OH− can react with CO2(aq) to form HCO\({}_{3}^{-}\)(aq). The generated bicarbonate ion may further accept a proton to become a carbonic acid molecule (H2CO3(aq)), or may lose a proton to become a carbonate ion (CO\({}_{3}^{2-}\)(aq)). The major carbon species in the solutions are CO2(aq), CO\({}_{3}^{2-}\), HCO\({}_{3}^{-}\), and H2CO3.

We compared the chemical speciation of the solutions under nanoconfinement and in the bulk phase at the same P-T conditions10. Figure 2 shows that at 1000 K, the mole percent of CO2(aq) in total dissolved carbon under nanoconfinement is 1.3 ± 0.9%, whereas it is 15.2 ± 2.0% in the bulk solution. The mole percent of HCO\({}_{3}^{-}\) under nanoconfinement (50.0 ± 1.0%) is higher than that in the bulk solution (35.9 ± 0.7%), and the concentrations of H2CO3(aq) are similar (42.7 ± 1.7% vs. 46.8 ± 1.5%). With increasing temperature from 1000 K to 1400 K, the mole percents of CO2(aq) under nanoconfinement and in the bulk solution increase to 14.5 ± 3.2% and 58.8 ± 2.0%, respectively. The equilibrium concentrations of CO2(aq) in the nanoconfined solutions are lower than those in the bulk solutions at the two temperature conditions studied here, suggesting that nanoconfinement promotes the CO2(aq) reactions. When increasing temperature along an isobar, due to thermal entropy effects, small molecules like CO2(aq) are more favored. We did not see obvious difference in reaction rates between nanoconfined and bulk solutions.

Fig. 2: Mole percents of carbon species in the CO2(aq) solutions in bulk and nanoconfined by graphene and stishovite (SiO2) at chemical equilibria.
figure 2

The initial mole fraction of CO2(aq) is 0.185. The pressure is ~10 GPa. The temperatures are a 1000 K and b 1400 K. The data of bulk solutions in (a) are from ref. 10, and the bulk data in (b) were interpolated using the simulation results in ref. 10. Uncertainties were obtained using the blocking method68.

Full size image

To understand why nanoconfinement promotes CO2(aq) reactions, we analyzed the water structure in Fig. 3. In the graphene-confined solutions, there are two sharp density peaks for oxygen atoms, corresponding to two water layers (Fig. 3a, b). We found that the carbon-containing ions and molecules in these two layers tend to align parallel to the graphene sheets (Fig. 4), and CO2(aq) mostly reacts with water molecules in the same layer (Supplementary Figs. 6 and 7). Nanoconfinement increases the probability of reactive encounters between CO2(aq) and solvent molecules, as diffusion is restricted to two dimensions33. It has been reported that the dielectric constant of nanoconfined water in the direction parallel to the confining surfaces (ϵ∥) increases significantly compared to the bulk value (ϵ0). In water, the Coulomb interaction between two ions is \(F=\frac{{q}_{1}{q}_{2}}{{\epsilon }_{0}r}\), where q1 and q2 are the charges of the two ions, and r is their distance. With increasing the dielectric constant, the magnitude of F decreases, so it is easier to separate a cation from an anion. Consistent with this, water molecules dissociate more easily under nanoconfinement23,25. In CO2(aq) solutions, the produced OH− ions from the water self-ionization are subsequently available to react with CO2(aq) (reaction (1)). The enhancement of ϵ∥ also further stabilizes HCO\({}_{3}^{-}\) and CO\({}_{3}^{2-}\) ions generated in the reaction between CO2(aq) and H2O or OH−. As a result, more CO2(aq) molecules react under nanoconfinement than in bulk. When the interlayer distance between graphene sheets increases beyond ~1.5 nm, the bulk behavior of water is recovered in the center of the slit pore and the effects of nanoconfinement become less obvious25.

Fig. 3: Number density profiles of oxygen atoms and carbon atoms along the z axis, normal to the confining surfaces.
figure 3

a, b The solutions under graphene confinement at 1000 and 1400 K, respectively. c, d The solutions under stishovite confinement at 1000 and 1400 K, respectively. Osol refers to the oxygen atoms in solutions (black lines), O\({}_{{{{{{{{{\rm{SiO}}}}}}}}}_{{{{{{{{\rm{2}}}}}}}}}}\) refers to the oxygen atoms in stishovite (orange lines), and C refers to carbon atoms (blue lines). The initial mole fraction of CO2(aq) is 0.185. The pressure is ~10 GPa. The center of confined fluids is set at z = 0, and the density distributions have been symmetrized. The horizontal black dashed lines represent the oxygen density in the bulk solutions at the corresponding P-T conditions44. e–g The H2O molecule, the OH−, and H+ ions bonded to the stishovite (100) surface, respectively.

Full size image
Fig. 4: Orientation distribution of sp2 carbon species in CO2(aq) solutions.
figure 4

The dihedral angle α is between the confinement interface and the plane defined by the three oxygen atoms in sp2 carbon species. a, b The solutions under graphene confinement at 1000 and 1400 K, respectively. c, d The solutions under stishovite confinement at 1000 and 1400  K, respectively. The pressure in all solutions are ~10 GPa. The initial mole fraction of CO2(aq) is 0.185. The center of confined fluids is set at z = 0, and the angle distributions have been symmetrized. e, f The CO\({}_{3}^{2-}\) and HCO\({}_{3}^{-}\) ions adsorbed on the stishovite (100) surface, respectively.

Full size image

Stishovite nanoconfinement

After studying the effects of graphene nanoconfinement, we turned to the confinement by a realistic mineral in deep Earth, stishovite, which is a stable phase of SiO2 (space group: P42/mnm) at the P-T conditions studied here45,46 and a major component of subducted oceanic crust47, playing a substantial role in transporting water into Earth’s mantle48. We exposed the cleaved stishovite (100) face, one of the low-energy surfaces49, to the carbon solutions as shown in Fig. 1b. We carried out constant-pressure (NPT) simulations to keep the pressure perpendicular to the solid–liquid interface at ~10 GPa, and then we found that the distance between the outermost oxygen atoms in two stishovite (100) surfaces is ~7 Å (Supplementary Table 1). In deep Earth, aqueous solutions released from subducting materials in devolatilization processes tend to locate at grain boundaries14. The size of the pores at grain boundaries is typically between 0.4 and 1.2 nm13, and our confinement width is within this range.

Figure 2 shows the chemical speciation of aqueous carbon solutions under stishovite confinement at ~10 GPa and 1000 ~ 1400 K. We found that at chemical equilibrium, 48.5 ± 2.1% (1000 K) and 19.1 ± 3.1% (1400 K) of carbon species, mostly HCO\({}_{3}^{-}\) and CO\({}_{3}^{2-}\), are bonded to the stishovite surfaces, unlike in the graphene-confined solutions. In the atomic density profiles shown in Fig. 3c, d, there are oxygen and carbon density peaks near the SiO2 surfaces, where the oxygen and carbon atoms come from the solutions, indicating that the solid–liquid interface plays an important role.

In bulk stishovite crystals, silicon atoms are octahedrally coordinated, and oxygen atoms are trigonally coordinated, whereas at the cleaved stishovite (100) surface, silicon atoms form bonds with five oxygen atoms, and each oxygen atom bridges two undercoordinated silicon atoms. Water molecules can directly bond to the stishovite surface, or dissociate under the influence of the surface. The hydroxide ion (OH−) from water dissociation can bond to an undercoordinated silicon atom to form a silanol (Si-OH) group, and the extra proton can bond with the surface oxygen atom to become a Si-(OH+)-Si bridge (Fig. 3e–g). Similar hydroxylation occurs at the quartz (1000) surface50,51. In our simulations, we found reactions between hydroxyl groups at the SiO2 surface and CO2(aq) in the solutions forming HCO\({}_{3}^{-}\). Figure 5a shows the reaction snapshots at 1000 K. The surface hydroxyl groups or the undercoordinated oxygen atoms also accept protons released in reaction (1), driving the reaction forward (Fig. 5b).

Fig. 5: Reactions of CO2(aq) catalyzed by the stishovite (100) surface.
figure 5

a The formation of HCO\({}_{3}^{-}\) at the interface. b The proton released from the reaction between CO2(aq) and water is accepted by the silanol (Si-OH) surface group.

Full size image

We analyzed the spatial orientation of the sp2 carbon species such as CO\({}_{3}^{2-}\)(aq), HCO\({}_{3}^{-}\)(aq), and H2CO3(aq) in the stishovite-confined solutions in Fig. 4c, d. We found that the molecular plane of the sp2 carbon species bonded to the stishovite surface tends to form an angle of ~60∘ or ~90∘ with the solid–liquid interface plane, dramatically different from the orientation of carbon species in the graphene-confined solutions. When the angle is ~60∘, the carbon species has two Si–O bonds by straddling two silicone atoms (Fig. 4e), while with the angle of ~90∘ the carbon species forms only one Si–O bond (Fig. 4f). The strong solid–liquid interaction substantially affects the molecular structure of the confined carbon solutions.

After analyzing the carbon species at the solid–liquid interface, we investigated the carbon species not bonded to the stishovite surface, i.e., fully dissolved in the stishovite-confined solutions. Figure 2 shows that at 1000 K, the mole percent of dissolved CO2(aq) is 10.5 ± 2.3%, which is larger than 1.3 ± 0.9% in the graphene-confined solution, and slightly smaller than 15.1 ± 2.0% in the bulk solution. At 1400 K, the mole percent of dissolved CO2(aq) in the stishovite-confined solution (39.8 ± 3.6%) is also between those in the graphene-confined (14.5 ± 3.2%) and bulk (58.8 ± 2.0%) solutions.

Both hydroxide ions and protons can be chemically adsorbed on the SiO2 surface, which affects the acidity of carbon solutions. Considering that the pH value of neutral water is no longer 7 under extreme P-T conditions, we calculated the difference between pH and pOH to quantify the acidity of solutions52:

$$f={{{{{{{\rm{pH}}}}}}}}-{{{{{{{\rm{pOH}}}}}}}}=-{\log }_{10}\left(\frac{[{{{{{{{{\rm{H}}}}}}}}}_{3}{{{{{{{{\rm{O}}}}}}}}}^{+}]}{[{{{{{{{{\rm{OH}}}}}}}}}^{-}]}\right),$$
(2)

where [H3O+] and [OH−] are the concentrations of hydronium and hydroxide ions, respectively. Because CO2(aq) reacts with water to generate H3O+, the solutions studied here are all acidic, i.e., f < 0, as shown in Table 1. The interesting finding is that the f value of stishovite-confined solutions is more negative than that of graphene-confined solutions at the same P-T conditions, which means that the former is more acidic than the latter, even though less CO2(aq) reacts in the stishovite-confined solutions. We have discussed that the stishovite surface adsorbs the hydroxide ions and protons from solutions. In addition, our AIMD trajectories show that the SiO2 surface favors the adsorption of OH− over that of H+ (see Supplementary Fig. 8); as a result, the stishovite-confined solutions are more acidic than the graphene-confined ones at the same P-T conditions. Increasing the concentration of H3O+ shifts the equilibrium of reaction (1) towards the left, so there is more CO2(aq) in the solutions.

Table 1 Acidity of aqueous carbon solutions: f = pH - pOH
Full size table

The nanoconfinement enhances ϵ∥, which stabilized charged ions, so in both graphene- and stishovite-confined solutions, more CO2(aq) reacts than in the bulk solutions. However, it has been reported that ϵ∥ near the hydrophobic surface increases more than near the hydrophilic surface, because the motion of water molecules is more hindered at the hydrophilic surface22. Considering that the stishovite surface is more hydrophilic than graphene, charged ions are less stabilized, so we found more CO2(aq) in the stishovite-confined solutions than in the graphene-confined solutions. This comes in addition to the reactions between SiO2 and solvent molecules, which make the fluids more acidic and thereby lead to destabilization of HCO\({}_{3}^{-}\) in favor of CO2(aq). Therefore, the CO2 concentration increase in the stishovite-confined solutions is a combined result of the hydrophilic confinement and the adsorption preference of OH− on the stishovite (100) surface.

In our simulations, we used the semilocal Perdew-Burke-Ernzerhof (PBE) exchange-correlation (xc) functional53, which was reported insufficient to describe aqueous systems at ambient conditions54; however, our previous studies showed that PBE performed better for the equation of state and dielectric properties of water55,56 and the carbon speciation in water7 at extreme P-T conditions than at ambient conditions. Particularly, we compared the simulations using PBE and a hybrid xc functional, PBE057. For an aqueous carbon solution at ~11 GPa and 1000 K, whose initial mole fraction of CO2(aq) is 0.016, both PBE and PBE0 suggest that HCO\({}_{3}^{-}\) is the dominant carbon species, and its mole percents are 79.8% and 75.0%, respectively7. Both PBE and PBE0 lack van der Waals (vdW) interactions, so we performed an additional simulation using the RPBE xc functional58 with Grimme’s D3 vdW corrections and the Becke-Johnson damping (RPBE-D3)59. For the solution confined by graphene at 10 GPa and 1000 K, we found that the mole percents of carbon species change by <6% (see Supplementary Fig. 4 and Supplementary Table IV). Particularly, the RPBE-D3 simulation gives that the concentration of CO2(aq) is 0%, while it is 1.3 ± 0.9% at the PBE level, indicating that our main conclusion that nanoconfinement enhances reactivity of CO2 is not affected by the neglect of the vdW corrections. VdW interactions do not play a major role in breaking and forming of covalent bonds, so do not much affect the chemical speciation studied here.

In summary, we performed extensively long AIMD simulations to study the chemical reactions and speciation of aqueous carbon solutions nanoconfined by graphene and stishovite at 10 GPa and 1000 ~ 1400 K. We found that the graphene nanoconfinement promotes the CO2(aq) reactions. When graphene is replaced by stishovite, less CO2(aq) reacts, but still more than in the bulk solutions. We found that contacting the stishovite (100) surface makes the solutions more acidic, which shifts the chemical equilibria, though the stishovite surface also catalyzes the CO2(aq) reactions by adsorbing HCO\({}_{3}^{-}\) and H+.

The enhanced reactivity of CO2(aq) in nanoconfinement has important implications for carbon transport and fluid-rock interactions in deep Earth. Aqueous fluids located at grain boundaries in minerals can either exist in isolated fluid-filled pores, or form a connected network of channels along grains facilitating fluid transport13. It is known that adding molecular CO2(aq) to water increases the rock-fluid-rock dihedral angle θ, which inhibits fluid flow60,61. However, our study shows that CO2(aq) reacts with water under nanoconfinement and also reacts with the solid interface, which may decrease θ and promote the interconnectivity of fluids14. Our study also sheds light on atomistic mechanisms of CO2 storage through mineral carbonation. CO2 reacts more in nanoconfined water, which benefits CO2 mineralization. If we choose minerals with larger points of zero charge than that of SiO2, such as forsterite62 and magnesium oxide63, the CO2 reactivity may be further enhanced.

Methods

We carried out Born-Oppenheimer ab initio molecular dynamics using the Qbox package64. We used periodic boundary conditions and employed plane-wave basis sets and norm-conserving pseudopotentials65,66, with a plane-wave cutoff of 85 Ry. The cutoff was increased to 145 Ry for pressure calculations. We applied density functional theory and the PBE exchange-correlation functional53. We sampled the Brillouin zone at the Γ point. We performed AIMD simulations in the canonical, i.e., NVT, ensemble. Stochastic velocity rescaling was used to control the temperature67, with a damping factor of 24.2 fs. We replaced hydrogen by deuterium to use a large time step of 0.24 fs in the simulations, but still referred to these atoms as hydrogen atoms.

We ran simulations for 180–480 ps after 20 ps equilibration to reach chemical equilibria (see Supplementary Table I for the simulation details). We analyzed the AIMD trajectories to determine the nature of carbon-containing molecules. For each carbon atom, we searched for the three nearest oxygen atoms, and sorted the C-O distances in increasing order. If the difference between the third and second C-O distance is <0.4 Å, the carbon species is a CO\({}_{3}^{2-}\) ion; otherwise, it is CO2. Hydrogen atoms were considered being bonded to their nearest-neighbor oxygen atoms. For the solutions confined by stishovite, the oxygen atoms were considered bonded to silicon atoms when the interatomic distance fell within the first peak of the Si-Oaq radial distribution function (RDF), i.e., 2.6 Å, as shown in Supplementary Fig. 5. We also varied the cutoff distances (0.4 and 2.6 Å) by ± 10%, and found that the changes of species concentrations are within the statistical fluctuations of our AIMD simulations (see Supplementary Tables V and VI).