The wetting of H₂O by CO₂

The boundary between two immiscible fluids represents a unique and fascinating regime. The properties of this biphasic interface, noticeably different from those of bulk fluid, arise from its inhomogeneous composition and the anisotropic distribution of its constituent particles.¹ The resulting asymmetry in forces acting between these particles promotes a series of interesting phenomena, including molecular orientational alignment, density gradients and structuring, phase transfer effects, and variations in reactivity and dielectric properties.^1–4 

Arguably, the most important biphasic interface is that of carbon dioxide and water, CO₂–H₂O. This interface exhibits a diverse range of physiochemical properties, which stem from the distinct phases that CO₂ can form with subtle changes to temperature and pressure. At ambient temperatures, CO₂–H₂O exists as a gas–liquid interface at pressures p less than the critical pressure (p_C = 73.8 bars; T_C = 304.1 K).⁵ At higher pressures, CO₂–H₂O forms a liquid–liquid (or supercritical–liquid) system, suggesting that temperature or pressure can be used as a lever to modulate biphasic properties. CO₂–H₂O finds use in a variety of real-word applications and phenomena, from compressible solvent design⁶ to ocean acidification monitoring⁷ to enhanced oil and gas recovery schemes.⁸ They are also central to various carbon capture and storage (CCS) techniques, where anthropogenic CO₂ is captured at source and injected under pressure into predominately aqueous subsurface storage sites, e.g., saline aquifers or disused coal seams.⁹ In these storage sites, retention of the anthropogenic carbon is mediated by the interactions between injected CO₂ and in situ H₂O, meaning that both storage capacity and the duration are functions of the interfacial behavior.

Realizing these applications for CO₂–H₂O requires a detailed knowledge of its biphasic interface. This includes understanding both its microscopic properties (in terms of the structuring of molecules, their orientations, the nature of their bonding, etc.) and its macroscopic properties, such as phase densities and interfacial tensions (IFTs), γ. The IFT, in particular, is critical to understanding biphasic interfaces. This property encapsulates much of the macroscopic nature of the interfacial region and is used to gauge their relative miscibilities. In terms of CCS and CO₂–H₂O, γ is directly proportional to the sealing capillary pressure, which quantifies the volume of CO₂ that can be stored and thus is a good estimate of the efficacy of a chosen storage site.^10,11

An extensive picture of both micro- and macroscopic properties of CO₂–H₂O has been built up over the past decades. In the experimental literature, researchers have focused mostly on the measurement of γ at differing temperatures and pressures.^12–22 In general, γ is observed to decrease with either increasing pressure or increasing temperature, although the exact nature of this reduction is dependent upon which CO₂ phase is present (see later, Fig. 2). Measurements from experiment, however, are subject to significant variations depending on experimental setup (e.g., thermocouple placement, equilibration times, and density treatment).¹⁰ Therefore, experimental work is often supplemented by theoretical and computational estimates for γ. Statistical perturbation models, such as statistical associating fluid theory (SAFT), do a good job in reproducing experimental values but can sometimes lack sufficient microscopic insight into the origins of these results.^23–25 In contrast, force field molecular dynamics (FF-MD) simulations provide both estimates for γ and full atomistic details on the behavior of CO₂–H₂O. Although the exact description of the CO₂–H₂O interface is sensitive to the choice of the empirical parameters, FFs have provided several important trends and insights. These include the following: CO₂–H₂O interfaces are molecularly sharp; molecules form a layered structure within several angstrom of the interface; water’s dipole orients parallel to the interface; and there exists a strong, lateral network of hydrogen bonds within the first contact layer of water.^10,26–30

Despite the valuable insight uncovered from these studies, a number of open questions remain surrounding CO₂–H₂O. For example, how exactly does γ behave in the vicinity of CO₂’s critical point (p_C, T_C)? How do high pressures impact the value of γ? In addition, what is the nature of interfacial CO₂ and H₂O molecules and how exactly do they interact with one another?

Building on previous work in the field, this paper aims to provide a new perspective for characterizing the CO₂–H₂O interface. Our work focuses on generating machine-learned potentials (MLPs), which we use to perform MD simulations with an ab initio level of accuracy. Unlike previously applied methods, our work breaks away from the empirical fitting to experimental or literature values and instead provides estimates from the ground up. MLPs already have a proven history of success in simulating complex systems,^31–37 and their use in simulating CO₂–H₂O would help address many of the open questions on this system. To date, however, an ab initio or MLP description of biphasic fluid interfaces has remained elusive. Biphasic systems, by definition, require larger simulation boxes than bulk for direct ab initio simulations and for generating training data for MLPs. At the same time, converging the properties of the bulk and interfacial regions requires long simulation time, beyond the reach of standard molecular dynamics.

In this work, we bring together ab initio accuracy and thorough sampling to accurately simulate large and complex mixtures. We utilize a number of distinct methodological features—specifically, MLPs for fast and accurate ab initio total energy predictions, an active learning strategy to build compact and representative training sets, and replica exchange simulations for enhanced sampling of the interface. Our approach allows accurate and converged analysis of both macro- and microscopic properties of the CO₂–H₂O interface. We obtain estimates of the interfacial tension at room temperature as a function of pressure. These are found to be in good agreement with the literature and help discern experimental reference data across the gas–liquid phase transition for CO₂. Improved predictions with respect to previous classical MD results, particularly for low pressure regions, are attributed to the greater accuracy of the underlying theory. Having obtained a correct description of the interfacial tension, we study the structural properties of the CO₂–H₂O interface to obtain key microscopic insights. Notable observations include the buildup of a liquid-like CO₂ layer at low pressures and a reduced structuring in the aqueous phase with increasing pressure. We anticipate the future application of our potential to studying other important state points for CO₂–H₂O as well as characterizing the dynamical system properties, e.g., in the form of diffusion coefficients and interfacial residence times. We believe that our approach will serve as a robust blueprint for investigating other complex biphasic interfaces.

Our approach to modeling the biphasic CO₂–H₂O interface comprises three established components: high-dimensional neural network potentials^38–40 (HD-NNPs) to act as our first-principles force generator; Query by Committee (QbC), an active learning technique for generating and optimizing our structural dataset;⁴⁰ and replica exchange molecular dynamics (REMD) to expedite the statistical convergence of interfacial properties over multiple thermodynamic state points. Combining these features allows us to apply accurate potential energy surfaces for simulating large, complex molecular systems over extended time periods. An overview of our approach is provided in Fig. 1.

Our work follows the committee NNP development procedure outlined in Ref. 40. We employ Behler–Parrinello HD-NNPs trained on DFT-level ab initio total energies and forces estimated with CP2K.^38,41 A set of hand-crafted radial and angular symmetry functions was used.⁴⁰ We employed the BLYP functional^42,43 augmented by Grimme’s D3 corrections.⁴⁴ This setup has been shown to closely reproduce the condensed-phase and structural properties of both pure CO₂ and pure H₂O.^36,45,46 In addition, BLYP-D3 replicates the critical and vapor–liquid coexistence properties for CO₂, making it an appropriate selection for treating the various CO₂ phases exhibited within our chosen pressure range.^36,46 Goedecker–Teter–Hutter (GTH) pseudopotentials were used for the treatment of the core electrons, TZV2P-GTH basis sets were used for the valence electron density, and a plane wave cutoff of 1050 Ry was used.

Training of our NNPs was implemented using QbC over multiple improvement cycles [see Fig. 1(b) for a schematic description]. This was performed using the open-source AML package.⁴⁰ A candidate dataset was provided in the form of chemical structures labeled with energies and atomistic forces. Structures were generated using a mixture of flexible force field potentials (SPC/Fw⁴⁷ + Zhu,⁴⁸ Lorentz–Berthelot combining rules) as well as preliminary NNPs across three active learning cycles. This gave a total of more than 100 000 structures. QbC was used to trim and optimize this dataset through the iterative selection of those structures displaying the largest uncertainties across each committee member. This was repeated for each trajectory, resulting in a refined dataset of some 4000 structures. Using this dataset, a final model was trained using the openly available N2P2 package.⁴⁹ 

This final, refined dataset is shown in Fig. 1(a). A breadth of structures is included here, including gas–liquid, liquid–liquid, crystalline, and pure-phase structures, in order to maximize exploration of the relevant configuration space. A committee of NNPs was trained using N2P2,⁴⁹ and the committee member displaying the lowest errors was selected as the final model. This final NNP was validated against the predictions of ab initio MD in the form of energies and forces (test RMSEs: 3.31 meV/atom and 88.26 meV/Å, of similar performance to that reported in Ref. 50 for pure water) as well as structural predictions of pure-phase properties, which are detailed in Sec. I of the supplementary material. Long-range tail corrections for the chosen NNP cutoff are expected to be on the order of 0.5–1 mN/m based on explicit classical model analysis,⁵¹ which is smaller than our usual statistical error obtained from block averaging. Due to these reasons, we have not applied any explicit tail correction to our data.

Typical system setups for our production run are shown in Fig. 1(c). Each system consisted of 600 water molecules and 200 CO₂ molecules. Lateral dimensions xy were fixed at 20 Å, and periodic boundary conditions were applied along all three axes. Atomic configurations were initialized such that the CO₂–H₂O interface resides parallel to the xy plane. System sizes and compositions were chosen to mitigate periodic error and finite-size effects, which have been shown in FF-MD to impact IFTs and other interfacial properties.^52,53 Results from our system-size analysis in which we varied lateral (xy) dimensions and molecular compositions are shown in Sec. II of the supplementary material. This analysis was performed using SPC/E + EPM2 for l_x = l_y = 12–30 Å and total number of molecules n = 168–3200. Our choice of box size and composition reflects the smallest system for which γ converges to within 0.5 mN/m of the largest box limit.

Fifteen representative pressures were selected spanning the 1–500 bar range. REMD was performed across 15 different pressures, from which we extracted macroscopic measurements (i.e., IFTs) as well as microscopic properties in the form of density profiles, molecular orientations, and radial distribution functions (RDFs). NNP-REMD was performed using the open-source i-pi package⁵⁴ with LAMMPS as the force generator and substituting hydrogen for deuterium atoms.^55,56 For each state point below 40 bars, simulations were run for 3 ns of simulation time, while above 40 bars simulations were run for 10 ns of simulation time corresponding to 3 × 10⁶ and 10⁷ steps, respectively. These were performed under the NP_zAT ensemble, with a fixed number of particles N, a fixed lateral (xy) area A, a fixed temperature (300 K), and P_z set to the target pressure. Collectively, these simulations compile 150 ns of ab initio-quality data, forming a robust basis for our analysis. We have checked the sensitivity of our results on the choice of DFT functional by also training an NNP model to revPBE-D3. Resulting IFTs prove to be very similar overall, and these are shown in Sec. III of the supplementary material.

Inspection of Fig. 2 yields several observations. At low pressures (p < 60 bars), our simulations predict a steep descent in γ with increasing pressure. This phenomenon is associated with the gradual accumulation of CO₂ at the water surface, a process which reduces the spatial anisotropy at the phase boundary and thus also the magnitude of γ (see later, Fig. 3). At intermediate pressures (60 $≤p<$ 100 bars), our NNP results indicate an abrupt change in the behavior of γ as we cross the critical pressure boundary (73.8 bars), going from gas–liquid to liquid–liquid, and finally, at high pressures (p ≥ 100 bars), our model yields a slight reduction in γ with increasing pressure. The comparatively small gradient in γ for this regime is attributed to a saturation in the concentration of CO₂ close to the interface such that changes to the pressure lead to minimal change in the near-surface environment. A minimum of 24 ± 2 mN/m is recorded at 500 bars.

Across Fig. 2, we then observe close agreement between experimental and our simulation results, within the range of pressures for which experimental data are available. This is true at both low pressures (p < 60 bars), where γ_NNP replicates the steep descent in gas–liquid IFTs, and high pressures (p ≥ 100 bars), where γ_NNP follows the gradual decline in liquid–liquid IFTs as reported by Hebach et al.¹³ and Bachu and Bennion¹⁶ A notable aspect of this work is the extension of IFT measurements beyond the currently available experimental pressure range for ∼ 300 K. Moreover, γ_NNP values recorded at intermediate pressures (60 $≤p<$ 100 bars) provide additional clarity on the nature of γ in the vicinity of CO₂’s phase transition pressure p_T (73.8 bars). This region has previously been a point of contention, with experimental results differing by up to 20 mN/m for pressures close to p_T, as seen by the selection of results displayed in Fig. 2. Our results clearly support the presence of a break in the slope between gas–liquid and liquid–liquid regimes, as reported by a number of authors.^13,16,17 There is no evidence of the “dip” or “cusp” in γ near p_T as predicted by some experiments.^12,14 Such observations are now thought to be experimental artifacts that stem either from differences in thermocouple placement or from the inherent uncertainty in near-phase-transition properties.^10,13

Our results also qualitatively corroborate the results of previous FF simulations. These are encompassed within the shaded portion of Fig. 2, the limits of which represent the extent of γ predictions for varying combinations of FF models (e.g., SPC/E + EPM2, SPC/E + PPL, and TIP4P/2005 + PPL).³⁰ While the exact value of γ is sensitive to the choice of models and parameterization, the overall range of FF predictions is in qualitative agreement with our NNP predictions, with a notable exception of the low pressure regime.³⁰ At these pressures, classical FFs systematically underpredict γ by 15–30 mN/m. We attribute the improved performance of our NNPs to the more accurate level of underlying theory compared to FFs and to the fact that our NNP is trained explicitly to replicate the ab initio-level interactions between CO₂ and H₂O. In comparison, bicomponent FFs are constructed using geometric/mathematical mixing rules, the choice of which is often ad hoc and can lead to substantial variations in the resulting thermodynamic properties.^58–60 This is exemplified by the large spread in γ values shown for FF-MD in Fig. 2.

To better understand the interfacial tension shown in Fig. 2, we plot a series of microscopic profiles detailing the variation in density ρ with distance d from the instantaneous interface. This interface is calculated at each time step using the Willard–Chandler formalism,⁶¹ allowing for a full resolution of the dynamic nature of the interface. Resulting profiles are plotted in Fig. 3 for both H₂O (red) and CO₂ (blue). An inspection of Fig. 3 reveals a stark contrast in the phase behavior of water and CO₂. Focusing first on the water phase, we observe a layered structuring in ρ across all pressures, in close agreement with the behavior of the air–water interface.³⁵ The shape of these profiles is relatively unvarying with changes to the pressure, although we do note an enhanced structuring of molecules at p ≤ 20 bars approaching atmospheric pressure. In contrast to this behavior, the CO₂ profile exhibits large changes in both magnitude and shape with changes to the pressure. At p < 10 bars, density profiles register a sharp peak at roughly 1 Å with a density that tails off exponentially with distance from the interface. Physically, this represents a film of monolayer CO₂ adsorbed at the water surface and a gaseous CO₂ phase extending beyond to large distance. At 20 $≤p<$ 60 bars, a second peak is observed at 4 Å, suggesting the formation of a bilayer of CO₂ at the interface. Beyond 60 bars, the exponential decay in density is replaced by liquid-like layering that tends toward bulk density for large distances. For the sampled pressure range, the fraction of CO₂ located within the aqueous component increases with increasing pressure. This is a product of both a greater CO₂ accumulation at the aqueous surface and a reduction in the interfacial sharpness at higher pressures. The calculated CO₂ solubilities are shown in Sec. I of the supplementary material and found to replicate the general trend in experimental results.

Relating these microscopic observations to the IFT profile of Fig. 2, it is clear that the variation in γ is a product of the considerable changes observed in CO₂ with changing pressure. Unlike water, CO₂ is a nonpolar molecule that exhibits relatively weak intermolecular interactions. Its critical pressure p_C resides within our range of sampling, so we end up simulating two distinct phases of CO₂ and a phase transition region. In comparison, at 300 K, water is firmly within the liquid phase, and its high natural surface tension—a product of its extensive hydrogen bonding network—ensures preservation of its structural integrity and prevents extensive mixing with the CO₂ phase. The most significant change in the structuring of water occurs at low pressures (p ≤ 20 bars). It is interesting to note that this change coincides with the emergence of the monolayer in CO₂. The combination of these observations possibly suggests that structuring of the aqueous phase is mediated by the extent of monolayer coverage of the adsorbed phase and that subsequent adsorbed layers have a minimal impact.

Our ab initio-level profiles share many similarities with previous force field modeling.^26,27 This includes the formation of monolayer and bilayer CO₂ and the recovery of liquid-like properties at high pressures. Interestingly, however, the pressure at which our model predicts the emergence of a CO₂ monolayer (and, therefore, also the starting of a second CO₂ layer) is much lower compared to conventional FF predictions (e.g., those of SPC/E + EPM2). In our results, a fully saturated monolayer forms at 20 bars; in FFMD studies, this occurs at higher pressures (around 60 bars for SPC/E + EPM2), thereby suggesting a lower wettability for CO₂.²⁷ These observations corroborate the ab initio work of Morishita and Shiga,⁶² who have suggested that differences in the interfacial behavior predicted by ab initio-level modeling and FF-MD stem from the weaker attractive CO₂–H₂O forces given by the latter. Similarly, we conclude that differences in the density profiles arise from the improved CO₂–H₂O descriptions provided by NNPs, which are trained to account for both polarization and charge transfer effects. On the other hand, the classical FFs used so far to study this interface are limited by the nature of their mixing rules and their rigid-body formulation, which we would not expect to fully represent the interfacial regime.

In Fig. 4(b), we show the distribution of angles between the instantaneous interface and the CO₂ orientation. At all pressures, our results suggest a clear preference for CO₂ to lie flat at the water surface. This corroborates the previous computational work, which suggests that interactions between water and CO₂ are maximized when the latter molecule lies parallel at the surface.²⁶ We note that the distribution is more pronounced at 90° at p < 20 bars, where the adsorbed CO₂ is more exposed and has fewer molecules to interact with. The addition of more adsorbed molecules at the aqueous surface reduces the strength of water–CO₂ interactions, allowing for greater freedom in the latter molecule’s orientational alignment with increasing pressure.

Biphasic interfaces represent complex and dynamic regimes. In this paper, we have highlighted an approach for analyzing these regions in a way that combines ab initio accuracy with converged statistics and computational tractability. Our methodology allowed for nanosecond treatment of large CO₂–H₂O systems, which has yielded a converged IFT profile at an ab initio level of accuracy. Our results show good agreement with the experimental literature for low- and high-pressure regimes, and a microscopic insight into this behavior has been provided. The reproduction of results in the vicinity of a phase transition (gas to liquid) is notable, given the difficulty associated with treating this highly dynamic regime. We observe the formation of a saturated CO₂ monolayer at low pressures (20 bars) with structural properties akin to those of bulk CO₂. The emergence of this monolayer coincides with reduced structuring for near-interface water molecules, suggesting that interactions between the first contact layers of water and CO₂ are critical for aqueous phase structuring. We envisage that such insights will be important for realizing CO₂–H₂O’s many applications as well as shedding light on other significant biphasic systems, for example, in biological membranes⁶³ and for liquid–liquid interfaces facilitating nanoparticle assembly.^64,65

In terms of carbon capture and sequestration, our NNP provides a robust tool for extending IFT measurements and providing benchmark figures for higher pressure regimes. In this way, coupled with additional geological measurements, one might use our model to provide estimates for a particular storage site of known temperature and pressure, thereby also providing an estimate of that site’s suitability. Compared with statistical perturbation models, our work has the advantage of providing an additional microscopic insight to supplement predictions of γ. Future work could look to exploit this microscopic insight toward estimating other important storage parameters, for example, the contact angle between CO₂ and water. It would also be interesting to investigate ideas related to thermodiffusion in underground storage sites. Previous classical work has suggested that CO₂ aggregates in areas of low temperature,^66,67 i.e., at the top of storage sites, as a result of the Soret effect. Investigating how these predictions change with NNP treatment would allow further comparison between classical and ab initio-level predictions and shed more light on this important phenomenon. A further degree of realism can also be added by incorporating Na⁺ and Cl⁻ ions into our training data and simulations. In doing this, we would better replicate the saline conditions of underground storage reservoirs and thus enable more accurate predictions of a site’s storage capacity.

Our work also unearths several differences between classical and ab initio-level modeling of the CO₂–H₂O interface. Compared with previous classical results, we suggest that previous analyses of low-pressure CO₂ coverage in aqueous systems may require reanalyzing. This could have profound implications for our understanding of processes such as CO₂ adsorption and its role in ocean acidification. Higher levels of CO₂ adsorption at the ocean surface imply a greater CO₂ dissolution and, therefore, a higher overall pH. Extending our current model through the addition of other relevant species, e.g., carbonate, bicarbonate, and carbonic acid, it would be interesting to investigate this acidification process and understand how the adsorption, dissolution, and subsequent reaction of CO₂ proceeds under differing conditions.

In conclusion, we provide new insight into the nature of biphasic interfaces and demonstrate the significance of interactions occurring between molecules adsorbed in the first contact layer. We anticipate application of the combined set of techniques used here to other important biphasic systems. The immediate findings of this study about preferred CO₂ layering at H₂O are expected to be of direct relevance in geoscience, climate research, and materials science.

See the supplementary material for NNP validation tests, system-size convergence tests, and a comparison of IFT predictions for NNP and classical MD.

SGHB was supported by the Syntech CDT and funded by EPSRC (Grant No. EP/S024220/1). C.S. acknowledges financial support from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) Project No. 500244608. V.K. acknowledges support from the Ernest Oppenheimer Early Career Fellowship and the Sydney Harvey Junior Research Fellowship, Churchill College, University of Cambridge. A.M. acknowledges support from the European Union under the “n-AQUA” European Research Council Project (Grant No. 101071937). We are grateful for computational support and resources from the UK Materials and Molecular Modeling Hub, which is partially funded by EPSRC (Grant Nos. EP/P020194/1 and EP/T022213/1). We are also grateful for computational support and resources from the UK national high-performance computing service, Advanced Research Computing High End Resource (ARCHER2), and the Swiss National Supercomputing Centre under project s1209. Access to both the UK Materials and Molecular Modeling Hub and ARCHER2 was obtained via the UK Car-Parrinello consortium, funded by EPSRC Grant Reference EP/P022561/1. Access to CSD3 was obtained through a University of Cambridge EPSRC Core Equipment Award No. (EP/X034712/1).

The authors have no conflicts to disclose.

Samuel G. H. Brookes: Conceptualization (equal); Investigation (lead); Formal analysis (equal); Writing - original draft (Lead); Writing – review and editing (equal). Venkat Kapil: Conceptualization (equal); Methodology (equal); Formal analysis (equal); Writing – review and editing (equal); Supervision (equal). Christoph Schran: Conceptualization (equal); Methodology (equal); Formal analysis (equal); Writing – review and editing (equal); Supervision (equal). Angelos Michaelides: Conceptualization (equal); Formal analysis (equal); Writing – review and editing (equal); Supervision (equal).

All data required to reproduce the findings of this study are available at https://github.com/water-ice-group/co2-on-h2o.