Structure, dynamics and stability of water/scCO2/mineral interfaces from ab initio molecular dynamics simulations

The boundary layer at solid-liquid interfaces is a unique reaction environment that poses significant scientific challenges to characterize and understand by experimentation alone. Using ab initio molecular dynamics (AIMD) methods, we report on the structure and dynamics of boundary layer formation, cation mobilization and carbonation under geologic carbon sequestration scenarios (T = 323 K and P = 90 bar) on a prototypical anorthite (001) surface. At low coverage, water film formation is enthalpically favored, but entropically hindered. Simulated adsorption isotherms show that a water monolayer will form even at the low water concentrations of water-saturated scCO2. Carbonation reactions readily occur at electron-rich terminal Oxygen sites adjacent to cation vacancies that readily form in the presence of a water monolayer. These results point to a carbonation mechanism that does not require prior carbonic acid formation in the bulk liquid. This work also highlights the modern capabilities of theoretical methods to address structure and reactivity at interfaces of high chemical complexity.

where E mol is the energy of gas phase molecule, E tot is the total energy of any combined system, E slab is the energy for the slab, E layer /E liq is the energy of the H 2 O or CO 2 layer, or a box of water or scCO 2 , CE is the cohesion energy, and n is the number of molecules in the simulation. Equation (1) estimates the binding energy of a single, gas-phase molecule to the surface. Equation (2) provides an estimate of the average binding energy to the surface per H 2 O or CO 2 molecule, including both molecule/surface and molecule/ molecule interactions. To de-convolute these energy terms, we introduce equations (3) and (4) the former provides estimates of the average layer-surface interaction energy, while the latter provides the average molecule-molecule interaction within the layer. Equation (5) estimates the cohesion energy in the pure liquids to juxtapose with the binding energy of the layer to the surface. Equation (6), gives BE4, the difference between BE1 and CE, which is a measure of the energetic preference of a molecule being adsorbed on the surface, as opposed to being solvated in the pure liquid phase. Figure 1 describes how the above quantities are related in a thermodynamic cycle. Based on the computed values of adsorption energies, the following conclusions can be drawn regarding water film formation. When a full monolayer of water is formed on the surface, the adsorption energy (− 97 kJ/mol) only mildly decreases compared to the adsorption of a monomer (− 103 kJ/mol). This is a result of strong H 2 O-surface interactions (− 82 kJ/mol) and comparatively weak H 2 O-H 2 O interactions within the layer (− 15 kJ/mol). The surface-H 2 O interaction can be mainly attributed to the 8 coordination bonds between the 4 top-most Ca cations and H 2 O. These interactions are very similar in magnitude to the hydration energy of the Ca (H 2 O) 6 2+ species 36 , which provides an average energy per Ca-water coordination bond of 164 kJ/mol. Using this quantity as guideline, we extrapolate that 8 Ca-H 2 O bonds provide an estimated − 93 kJ/mol each. This is slightly higher than the computed quantity BE2, and suggests that the surface Ca-H 2 O interactions are weaker than those in the fully-solvated Ca 2+ . Similarly, the cohesion energy CE of bulk water (− 57 kJ/mol) arises from approximately 3.5 H-bonds per H 2 O molecule, averaging − 16.3 kJ/mol each.
In the case of CO 2 , the average binding energy of a CO 2 in the monolayer is lower (− 17 kJ/mol) than that of a single CO 2 molecule (− 34 J/mol). More specifically, while the Ca-O CO2 17 , or O surf -C CO2 interactions have a stabilizing effect 37,38 , the OH-O CO2 interactions are only mildly stabilizing 39 . Thus, CO 2 -CO 2 interactions in the monolayer provide only a small stabilization of − 4 kJ/mol and the CO 2 layer is only weakly bound. Overall, there is a stronger enthalpic driving force for H 2 O adsorption on the mineral surface than the appreciably weaker one for the CO 2 , suggesting that in carbon sequestration scenarios, water saturated CO 2 will likely undergo preferential segregation, with H 2 O gravitating towards the surface forming a water film. However, as we will see shortly, entropic considerations may radically influence the speciation at the interface.
Finite Temperature Effects. We assess the impact of the hydration level on mineral dissolution and carbonate formation. For the analysis, we consider the progression of systems ranging from: pure scCO 2 , scCO 2 with one, two and three water layers, and pure water, snapshots shown in Fig. 2. A summary of cell parameters and simulation times for all the systems can be found in Supplementary Table 1. In Fig. 3, we show the molecular density profiles of Ca, H 2 O and CO 2 along the direction of the surface normal. The top two panels show the distribution of CO 2 (sc-CO 2 ) and H 2 O (Water) for the surface/pure liquid systems. For pure scCO 2 , we observe a liquid-like layer (calculated density ~ 0.6 g/ml) that is consistent with a super-critical CO 2 medium. Higher density peaks, at ~2 Å from the surface, are indicative of a denser CO 2 layer, also observed in the case of other minerals 38 . The bottom 4 panels in Fig. 3 show distinct water layer features with or without scCO 2 . At the H 2 O/scCO 2 interface, we still observe a higher density CO 2 layer, however it is less ordered compared to the pure scCO 2 system. In all cases, diffusion in and out of the CO 2 layer is evident on the ps time scale. On the Ca-deficient surface, within 5 ps of simulation, formation of carbonates (CO 3 2− ) is observed from reaction of CO 2 with surface oxygens O 2 − adjacent to Ca 2+ vacancies. Conversely, on the Ca-rich face, neither carbonate formation nor strong Ca 2+ /CO 2 interactions are observed on the simulation time scale. However, as we will discuss later, carbonates can be formed when Ca vacancies occur near the surface.
The water-only system (Water) shows a striation of the molecular density, with a pronounced maximum adjacent to the surface. The water monolayer (1W) is stable at finite temperature and perturbed   only slightly by the presence of the scCO 2 component (1W/CO 2 ). Within this layer, water exhibits a diffusion coefficient D ~ 10 −6 cm 2 /s only an order of magnitude smaller than water self-diffusion in water (22 × 10 −6 cm 2 /s) 40 , indicating a mobile water layer. Although water molecules diffuse within this layer, they do not diffuse out of it. Within the layer, each H 2 O is bound to an average of 1.4/1.4 (w/wo scCO 2 ) waters and 1.0/1.0 surface oxygens via hydrogen bonds. Each Ca 2+ is bound to an average of 1.5 water molecules through direct Ca-O w coordination bonds. Compared to the 1W system, the water coordination environment changes in the 2W and 3W systems. Within the first layer, each H 2 O is bound to an average of 2.0/2.0/2.0 (2W/3W/Water systems) waters and 0.2/0.1/0.2 surface oxygens via hydrogen bonds. Each Ca 2+ is bound to an average of 1.6/1.1/1.1 waters through Ca-O surf bonds. The 2 nd and 3 rd water layers show a distinctly different intermolecular structure, very similar to that of the bulk water. For instance, the 2 nd water layer shows a broader density maxima than the first, even exhibiting a double peak in the 2W system similar to other studies 11 . The 3 rd water layer shows progressively lower density maxima, indicating that a finite percentage of waters, ca 2%, reside in a 4 th layer. Diffusion coefficients and densities for all these systems are summarized in Supplementary Table 2.
The presence of a water layer, adjacent to the surface, induces large disorder in the Ca positions such that these atoms exhibit finite populations up to 1 Å above their normal lattice positions, see Fig. 3. To examine this further, we performed a simulation with a single Ca atom extracted from the surface into the first water layer of the 1W system. The resulting structure ( Fig. 4a and Supplementary Fig. 4) consists of a water-solvated Ca 2+ species within the 1W film as an outer sphere complex. Using the Blue-Moon sampling method 41 , we examined the free energetics of the dissolution process by computing the potential of mean force (PMF) along the reaction coordinate defined by the Ca 2+ z-coordinate, i.e. height above the surface plane. We find that the free-energy barrier for dissolution is 30 kJ/mol indicating that at 323 K this process will be relatively fast. The change of Helmholtz free energy, Δ A, between the two minima is 14 kJ/mol demonstrating that at this temperature there will be an equilibrium between inner and outer sphere Ca 2+ complexes with a K eq = exp(− Δ A/k B T) ~ 10 −3 . Note that this complex shows a positive increase in free energy as it moves out into the scCO 2 layer, as long as the first coordination sphere of Ca 2+ consists only of water. However, complexes with mixed coordination (H 2 O and CO 2 ) may further enhance Ca transport away from the surface 17,42 . Ca 2+ dissolution results in Ca-vacancies where carbonate formation is expected to be spontaneous, in a process similar to sulfite formation 35 . To verify this, we computed a PMF by moving a CO 2 from the scCO 2 liquid across the water layer towards the Ca vacancy to estimate the free energy of carbonate formation. As can be seen in Fig. 4b, this process is effectively barrier-less and exergonic by 17 kJ/mol, implying that carbonation can readily occur and is limited only by diffusion and the available population of Ca-vacancy sites. This emphasizes the paramount importance of a stable water monolayer in establishing the nucleation sites for carbonate formation.
Water adsorption/desorption to/from Anorthite. We now will focus on extracting information about water transport to, from and across these interfaces. This phenomenon is particularly relevant to reactivity under geologic sequestration conditions, as water-saturated scCO 2 can lead to unexpected events, such as steel corrosion in water-lean scCO 2 33,34 or preferential sulfite formation in carbonate minerals 35 . We computed the PMF for moving a single H 2 O from the middle of the liquid slab through the box, such that its immediate environment changes from being solvated in the scCO 2 to being bound on the Ca-rich surface 41 . This process was repeated for the following systems: 1H 2  Supplementary Section S4 provides details on the PMF convergence. The associated changes in the Helmholtz free energy, Δ A, are given in Fig. 5, and the computed coverage dependent adsorption free energy, Δ A(θ), is presented in Fig. 6a.
The free energy landscape for the 1H 2 O/scCO 2 system is essentially flat, up until H 2 O is ca 5 Å away from the anorthite surface, see Fig. 5 and Supplementary Fig. 6. This distance corresponds to the water approaching the periphery of the high density CO 2 layer. The free energy decreases monotonically up to 1.7 Å above the surface where the water is stabilized by forming both Ca-O w and hydrogen bonds with  Supplementary  Fig. 6. Because the constrain is imposed on only one of the waters, the density and free energy profiles are very similar to those for 1 H 2 O, with the other H 2 O molecule(s) moving freely. Adsorption occurs step-wise with the H 2 O dimer/trimer breaking apart along the way and with the constrained water molecule adsorbing first. Subsequently, in both the dimer and trimer simulations, the remaining molecules gravitate towards the surface adsorbed H 2 O: once the additional waters become engulfed in the scCO 2 layer, they adsorb to the surface adjacent to the constrained H 2 O molecule without any additional constraints. These simulations imply that water film growth begins with a single water molecule as a nucleation site followed by the addition of subsequent molecules, to form an island.
To quantify this observation, we focus on the Δ A(θ) of the reactions: The resulting free energy Δ A(θ), as a function of coverage θ (θ = n/14), is plotted in Fig. 5a. We note that in order to directly compare the Δ A(θ) values obtained from the PMF for the n = 2 and 3 simulations (which are relative to solvated clusters) with those of the n = 1, 7, 10 and 14 (which are relative to solvated monomers) we adjust the free energy of the cluster in the solvated phase to be relative to n solvated water monomers at infinite dilution, see Supplementary Section S6 for details. In accord with our above conclusion regarding H 2 O-H 2 O interactions, the computed Δ A(θ) values grow steadily from n = 2-14 by about − 20 kJ/mol. Note that this value is on the order of the molecule-molecule interaction energy (BE3 = − 15 kJ/mol) reported in Table 1. We proceed by de-convoluting Δ A(θ) = Δ E(θ) − TΔ S(θ) into its entropic (TΔ S(θ)) and energetic (Δ E(θ)) components in order to understand the role played by the H 2 O-H 2 O interactions in driving film formation. To assess Δ E(θ), we consider that adsorption of water at the anorthite surface involves replacing adsorbed CO 2 molecules with adsorbing waters at a ratio of 9CO 2 /14H 2 O based on the full coverage of the of CO 2  This process can be subdivided into the sum of two half reactions: The energy Δ E(θ) of (9) and (10) can be easily estimated based on arguments similar to those discussed earlier. As shown in the thermodynamic cycle of Fig. 1, for H 2 O, the energy of (9) can be decomposed utilizing the binding energies of Table 1, while introducing an intermediate state (H 2 O(g)) and the solvation energy of water in scCO 2 (estimated to be − 20 kJ/mol 39 ). Neglecting the molecule-molecule interaction term (BE3) and using only the average surface-molecule adsorption energy (BE2) we obtain: → H 2 0(ads) ΔE= -62kj/mol (9)=(11)+ (12) In an analogous manner, the energy for (10) is found to be Δ E = 6 or 4 kJ/mol, with or without CO 2 -CO 2 interactions. As such, the energy for (8) Fig. 6a. The relationship Δ E(θ) − Δ A(θ) = TΔ S(θ) allows us to estimate the entropy change during adsorption, which agree fairly well with values obtained using quasiharmonic approximation 43 , see Supplementary section S7 for detailed method and results. In all cases, the entropy of adsorption (Fig. 6a) is negative due to a loss of the degrees of freedom of the scCO 2 -solvated H 2 O upon adsorption. The function TΔ S(θ) becomes more negative by − 27 kJ/mol between the adsorption of a monomer and a dimer, then remains relatively flat as the cluster of adsorbed waters increases in size up to a full monolayer. This occurs because addition of water molecules to a surface-bound cluster requires an increase in entropy of adsorption due to the restriction of the number and locality of binding sites. This factor is roughly unchanged for n = 2-14 and hence TΔ S(θ) is relatively constant as the layer forms. As such, Δ A(θ) decreases roughly at the same rate as Δ E(θ), and is proportional to the increase in the H-bonding. Our simulations indicate that the barrier for nucleation of a 1W water layer, via an island growth mechanism, arises due to the entropic considerations and has the strongest influence at n = 2-3 where the number of H-bonds are relatively few.
Finally, our results also suggest that water layer formation will be kinetically hindered due to the lack of a strong thermodynamic driving force at low θ. To assess this quantitatively, we have computed an adsorption isotherm 44 employing the θ -dependent equilibrium constants, K(θ) = exp(− Δ A(θ)/k B T), see Fig. 6b. From this isotherm, we solve the equation to estimate the expected water coverage as a function of water concentration. We find that water concentrations below 10 −5 %mol result in very low surface coverage (θ < 0.1), while for water concentrations above 10 −3% mol, a full coverage water monolayer is expected. In the intermediate concentration range of 10 −5 −10 −3 , consistent with the average 10 −4 %mol concentration of water in wet scCO 2 at P = 90 bar, T = 323 K 45 , the non-monotonic behavior of Δ A(θ) leads to multiple solutions to the equation where both a high and low coverage regions are stable. This is interpreted as the likely occurrence of an equilibrium distribution of high and low coverage patches on the surface. In the context of Ca 2+ extraction and carbonate formation, this implies that the partial water layer formation will impact the population of Ca 2+ vacancies and thus be a major kinetic bottleneck for carbonate formation.

Discussions
We have demonstrated that in the case of a prototypical mineral surface, here anorthite, in contact to water-saturated scCO 2 , a liquid, boundary layer of water will form and facilitate the extraction of Ca 2+ cations in the form of stable, hydrated, outer-sphere complexes. This process creates cation vacancies at the top-most layers, which can serve as nucleation sites for carbonate formation limited only by the water layer formation and CO 2 diffusion. Water film formation was found to be energetically driven by both Ca 2+ -H 2 O and H 2 O-H 2 O interactions, but entropic contributions largely override the stabilizing factors. The entropic component arises from H 2 O adsorption on the surface as well as the availability of binding sites. Our results highlight that the existence of a boundary water layer, significantly different in composition to the bulk liquid phase, can occur at even relatively low concentrations, as in the case of water-bearing scCO 2 . In the present case, it also points to a novel mechanism for carbonation, distinct from the traditional view that carbonation must begin with CO 2 dissolution in the water phase. We have shown that a molecular level understanding of enthalpic and entropic contributions is now achievable through a combination of simple and intuitive chemical arguments and large-scale AIMD simulations. With this capability, we can now access an unprecedented level of qualitative and quantitative molecular scale understanding of structure, dynamics and reactivity at a chemically complex solid-liquid interface which may be difficult to obtain by direct experimental observation alone.
Scientific RepoRts | 5:14857 | DOi: 10.1038/srep14857 Methods DFT parameters. Density-functional-theory(DFT) based AIMD simulations were carried out with periodic boundary conditions (3D PBC) with the exchange correlation functional of Perdew, Burke and Ernzerhoff (PBE) 46 as implemented in the CP2K package 25 . Dispersion corrections were accounted for by Grimme's third generation corrections DFT-D3 47 that have been evaluated extensively in the literature for both scCO 2 and H 2 O and shown to adequately represent weak interactions in these systems 16,17,38,[48][49][50] . The core electrons were described by the norm-conserving pseudopotentials 51 , while the valence wavefunctions were expanded in terms of double-zeta quality basis sets optimized for condensed systems to minimize linear dependencies and superposition errors 52 . An additional auxiliary plane wave basis set with a 350 Ry cutoff was used for the calculation of the electrostatic terms. We used a time-step of 1.0 fs and the Nosé-Hoover thermostat for NVT ensemble calculations at 323 K.
Computational models. Anorthite, with the chemical composition of CaAl 2 Si 2 O 8 , is the calcium end-member of the plagioclase feldspars. A coupled substitution of two Al 3+ for two of Si 4+ , allows the incorporation of one Ca 2+ and overall neutrality 53  The mineral-liquid interface was modeled with 2 different simulations of the anorthite slab with 32 CO 2 or 84 H 2 O molecules in the interlayer space, determined by optimization of cell parameters at 90 bar. MD simulations within the NVT statistical ensemble were performed at 323 K, resulting in a pronounced peak near the surface in the density profile plot indicating formation of strong first H 2 O/CO 2 layer near the interface. Based on this analysis, the first layer of H 2 O/CO 2 was determined to contain 14H 2 O/9CO 2 molecules respectively and was used to conduct simulations for stand-alone monolayer properties.  Table 1 summarizes the optimized simulation cell parameters and simulation times for each system. Estimates of the diffusion coefficients and densities in the liquid slab of each simulation are summarized in Supplementary Table 2. Overall, it was found that the scCO 2 layer in our model adequately reproduces the structure, density and dynamics of an scCO 2 liquid phase at ~90 bar and T = 323 K.

AIMD simulations.
We have employed the Blue-Moon ensemble method 41 to calculate the potential of mean force (PMF) required to move a H 2 O molecule from the middle of the liquid region towards the Ca-rich face of the anorthite (001) surface. This method is broadly used to compute the free energy profile along a reaction coordinate direction, which is not likely explored with MD simulation due to potentially high barriers. For the initial and final points, we also performed unconstrained MD simulations to locate the thermodynamic minima. The average forces subject to geometric constraints were calculated as ensemble-averages, and converted to a PMF based on the Blue-Moon ensemble method 41 : the z-distance between the surface and the water Oxygen atoms was chosen as the geometric constraint, and all other parameters were allowed to vary. This approach allowed the H 2 O molecule to move laterally over the mineral surface through the liquid boundary layer. Integrating the average force along the reaction coordinate generated the free energy profile. Increments in dz were chosen to be 0.2-0.5 Å apart depending on the steepness of the PMF to better estimate the free-energetics via thermodynamic integration, which resulted in a total of 11-18 simulations per PMF. The force on each point in the PMF was converged to within ± 0.0004 a.u. See Supplementary section S4 for details.
Each AIMD trajectory, both constrained and unconstrained, was equilibrated for approximately 5-10 ps and data was collected for analysis from at least 10 ps of well-equilibrated trajectory. The sum total of the overall compiled AIMD trajectories accounts for approximately 2.0 ns of total simulation time (system-specific simulation times are recorded in Supplementary Table 1). This data set represents over 10 6 DFT configurations for systems of total size of ~500 atoms.