Carbonation Reaction Mechanisms of Portlandite Predicted from Enhanced Ab Initio Molecular Dynamics Simulations

: Geological carbon capture and sequestration (CCS) is a promising technology for curbing the global warming crisis by reduction of the overall carbon footprint. Degradation of cement wellbore casings due to carbonation reactions in the underground CO 2 storage environment is one of the central issues in assessing the long ‐ term success of the CCS operations. However, the complexity of hydrated cement coupled with extreme subsurface environmental conditions makes it difficult to understand the carbonation reaction mechanisms leading to the loss of well integrity. In this work, we use biased ab initio molecular dynamics (AIMD) simulations to explore the reactivity of supercritical CO 2 with the basal and edge surfaces of a model hydrated cement phase— portlandite—in dry scCO 2 and water ‐ rich conditions. Our simulations show that in dry scCO 2 conditions, the undercoordinated edge surfaces of portlandite experience a fast barrierless reaction with CO 2 , while the fully hydroxylated basal surfaces suppress the formation of carbonate ions, resulting in a higher reactivity barrier. We deduce that the rate ‐ limiting step in scCO 2 conditions is the formation of the surface carbonate barrier which controls the diffusion of CO 2 through the layer. The presence of water hinders direct interaction of CO 2 with portlandite as H 2 O molecules form well ‐ structured surface layers. In the water ‐ rich environment, CO 2 undergoes a concerted reaction with H 2 O and surface hydroxyl groups to form bicarbonate complexes. We relate the variation of the free ‐ energy barriers in the formation of the bicarbonate complexes to the structure of the water layer at the interface which is, in turn, dictated by the surface chemistry and the degree of nanoconfinement. the density profiles characteristic of the crystal atoms, shown the (001) surface. The C and O c peaks strongly overlap with the surface O h and H atom peaks. We attribute the C and O c peaks on the surface region to the reacted (bi)carbonate


Introduction
The development of effective and sustainable approaches for the reduction of global carbon footprint is crucial for combating climate change and its negative effects. The technology of geological carbon capture and sequestration (CCS) has demonstrated an enormous potential in this respect [1][2][3][4][5]. This technology, however, poses several scientific questions which are important for its long-term success. One of the major concerns is the degradation of cement wellbore casings which are used to prevent migration of fluids to the surface through leakage pathways [6][7][8][9][10][11][12]. Although leakage pathways can occur in wells due to faulty construction and other mechanical defects, geochemical reactions induced by the injected fluids could cause cement degradation, resulting in damage of wells and the development of leaks [6,7,[13][14][15]. Therefore, to make accurate and reliable assessments of the near-and long-term performance of wellbore cements when exposed to injected CO2, it is important to understand geochemical interactions between CO2 and cement at relevant subsurface operation conditions. The carbonation reactions of the injected CO2 with wellbore cement pose serious well integrity issues because of the acid attack on cement [16]. CO2-saturated fluids are thermodynamically incompatible with cement, since cement is highly alkaline (pH > 12. 5) and not in chemical equilibrium with the acidic CO2-rich fluids (pH < 6) [17]. Depending on the subsurface conditions (e.g., temperature, pressure, salinity), studies report different degrees of cement degradation due to CO2-rich fluids [11,13,[18][19][20][21]. For instance, in dry supercritical CO2 (scCO2) conditions, carbonation of cement results in decrease of permeability and a slight increase in hydrated cement strength [19,20,22]. On the other hand, in the presence of water and/or brine, cement undergoes a progressive decalcification process, which leads to a drastic loss of mechanical strength and an increase of permeability [19,20,[22][23][24][25]. Nevertheless, due to the complexity of the subsurface environmental conditions, it is difficult to experimentally disentangle various factors controlling the carbonation of wellbore cement. In this regard, atomistic computer simulation approaches have proven to be very useful because of their ability to provide atomic-scale quantitative information about the physical and chemical processes occurring during carbonation [26][27][28][29].
The degradation of cement due to carbonation reactions mainly affects portlandite (Ca(OH)2) and calcium silicate hydrate (C-S-H) phases, which both constitute about 70-85% of hydrated cement [17]. In this work, we have focused on the reactivity of CO2 with the portlandite phase. Portlandite is a hydrous mineral which is very reactive with CO2, making it a good model for understanding carbonation reaction mechanisms in cement and in other hydrous minerals with similar structure (e.g., Mg(OH)2). The carbonation reaction mechanism of portlandite has been a subject of many recent studies using experimental [30][31][32][33][34][35][36][37][38] and atomistic simulation [28,39] techniques. It is proposed that carbonation of portlandite is based on either a solid-state mechanism [34,37,40] or via coupled dissolution-precipitation reaction mechanism in a water-rich environment [31]. Although the reaction products are well known, the rate at which these two reaction mechanisms proceed critically depends on the presence of water, pressure and temperature [40]. Galan et. al [37] highlight that at ambient conditions, gas-solid carbonation reaction is self-arresting due to the formation of a protective carbonate layer around the portlandite crystal. The same observation is reported by Longo et. al [29] on the carbonation of wollastonite (CaSiO3). In contrast, Vance et. al [30] suggests that so long as suitable reaction conditions are provisioned, the direct carbonation of Ca(OH)2 and ordinary Portland cement (OPC) can progress substantially due to the volume change accompanying the reaction and the release of water which favors advancement of the reaction.
To provide a better understanding of the carbonation reaction mechanism of portlandite, this work examines and compares the reactivity of basal and edge surfaces of portlandite with scCO2 and H2O/CO2 mixture using enhanced ab initio molecular dynamics (AIMD) simulations. In AIMD simulations, the forces are calculated on-the-fly from accurate electronic structure calculations, and thus the technique is a reliable approach to predict reaction pathways and the energetics of CO2-portlandite interactions. However, AIMD calculations are computationally very demanding, which limits the model system size to only a few hundred atoms and the timescales to only a few tens of picoseconds. Since the picosecond timescale is appropriate to simulate only barrierless processes, the dynamics of the rare reaction events was accelerated by means of the metadynamics approach [41,42]. This method allows us to compute the free-energy landscape for the carbonation reaction of CO2 with the portlandite surfaces. We show that in dry scCO2 conditions, the edge surfaces of portlandite undergo a rapid barrierless carbonation reaction resulting in the formation of a carbonate surface layer. In contrast, the presence of water hinders direct interaction of CO2 with portlandite surface, with the free-energy barrier being governed by the thickness of the interfacial water layer.

Computational Details
To model the reaction mechanisms of CO2 with portlandite crystals, we used the ab initio molecular dynamics (AIMD) simulation technique [43][44][45]. AIMD is a versatile tool which combines both first principles and molecular dynamics concepts where the interatomic forces are calculated based on quantum mechanics (electronic structure theory) and the motion of the nuclei is described classically according to Newton's mechanical laws. This allows realistic simulations to be performed without adjustable parameters. Nevertheless, an accurate description of the processes is highly dependent on the different approximations employed in the electronic structure calculations, such as the exchange-correlation functional approximation, the description of the dispersive (van der Waals) interactions as corrections, and the basis set limitations.
The characteristic times accessible by AIMD simulations are small-on the order of picoseconds-and thus are not often sufficient to capture chemical reactions which occur as are rare events at longer timescales. To overcome this potential limitation, the metadynamics enhanced sampling method [41,42,[46][47][48] was used to facilitate the reaction events and evaluate the free energy of the CO2 carbonation reactions. Metadynamics enables the acceleration of conformational transitions by iteratively filling the potential energy surface by a sum of Gaussians centered along the trajectory of a set of collective variables [48]. Thus, its effectiveness depends on the appropriate choice of collective variables defined to describe the processes of interest.
The AIMD simulations were performed within the density functional theory (DFT) using the Gaussian and plane-wave basis approach, as implemented in the CP2K simulation code [49]. The core electrons were described using the Goedecker-Teter-Hutter (GTH) [50] pseudopotentials and the valence density was developed on a double-zeta valence polarized (DZVP) basis set along with an auxiliary plane-wave basis set with cutoff energy of 400 Ry obtained after a series of convergence tests (see Figure S1 in the supporting information). The generalized gradient approximation (GGA) parametrized by Perdew, Burke and Ernzerholf (PBE) [51] was used for the exchange-correlation terms with the dispersion interactions included using the Grimme D3 [52] corrections. The nuclei dynamics was treated within the Born-Oppenheimer approximation and the convergence on forces was chosen to be 10 −6 a.u. The timestep for the integration of the equations of motion was set to 0.5 fs. All calculations were performed in the NVT statistical ensemble (constant temperature T, volume V, and the number of particles N) at a temperature of 323 K which is maintained using the Canonical Sampling through Velocity Rescaling (CSVR) thermostat [53]. The systems were pre-equilibrated for 5 ps before a production run of 30 ps was performed for further statistical property analysis.
The metadynamics technique was used to calculate the reaction free energies using the PLUMED code [54]. To predict the reactivity of CO2 molecules with the portlandite surfaces in pure supercritical state and in a CO2/H2O mixture, we chose the carbon atom (for all the CO2 molecules) coordination number (CN) to oxygen atoms of both the surface hydroxyls (Oh) and of the neighboring H2O molecules (Ow) as the collective variable (CV). This choice of the collective variable ensures that the reaction barrier for sp3 hybridization of the carbon atom is estimated and does not include any proton transfers in the general reaction coordinates. The calculation of the atom coordination numbers is defined by the switching function, Cij(rij): where rij is the instantaneous distance between atoms i and j, r0 is the distance beyond which the bond is considered to be broken and d0 is the central value of the switching function (i.e., Cij(d0) = 1). m and n are parameters of the switching function chosen in such a way that Cij(rij) tends to zero beyond r0. Gaussian hills with the width of 0.02 and an initial height of 3.3 kJ/mol were added every 50 fs. The metadynamics simulations were ran for 30 ps, during which only forward reactions were observed. However, the free energy surface (FES) curves were constructed till the first barrier crossings and therefore the free-energy barrier is estimated from the bias potential accumulated after a single barrier-crossing event. This implies that while the reactant regions of the paths are reproduced properly by the bias potential, the product regions do not feature their real depths. This strategy of cutting the FES at the moment of crossing commonly used in DFTbased metadynamics simulations is aimed at dealing with hysteresis (where the barrier keeps going down once the system is in the product well) [55]. The information on the variation of the CVs with time is provided in the supporting information ( Figure S2). The free energy (F) within the different states in the FES is obtained by integrating the (FES) F(s), with respect to the collective variables (s) as follows [41,56]: where , is the Boltzmann constant, T the temperature of the system and F(s) defines the free-energy surface.
To determine the portlandite surfaces exposed to the environment, the surface energies () of various low-index surfaces were calculated. The surface energies were calculated from the energies obtained in the DFT calculations of the bulk portlandite crystal and the surface slabs using the following relation: where Eslab is the total energy of the slab and Ebulk is the bulk total energy, n is the number of bulk unit cells used to construct the slab, A is the area of the generated surfaces and the factor of 2 accounts for the two equivalent surfaces in the slab. Based on the Wulff theorem [57], the equilibrium morphology of the portlandite crystal was obtained from a polar plot of the surface energy versus the orientation of the surface normal vector using the WulffMaker software [58]. The relative contribution of a crystal facet to the morphology is dependent not only on the surface energy but also on its orientation [57,59].

Atomistic Models
Portlandite, Ca(OH)2, is a layered mineral, with the layers consisting of distorted edge-sharing CaO6 polyhedra. Each hydroxyl group is connected to three calcium atoms in its layer and surrounded by three other hydroxyl groups belonging to the adjacent layer (see Figure 1a). The portlandite crystal model used in this work is characterized by the hexagonal crystal structure with lattice parameters a = b = 3.589 Å, c = 4.911 Å,  =  = 90° and  = 120° [60]. The lattice parameters obtained after optimization are a = b = 3.696 Å, c = 5.147 Å,  =  = 90° and  =120°, within 5% of the experimental parameters.
To build the surface models, supercells representing the crystal were created by multiplying the unit cell by 4 × 4 × 3 along the a, b, and c crystallographic directions, respectively. The supercells were then cleaved along the planes (001), (100), 11 0 , (110), (101) and 11 1 , exposing various crystal facets to investigate the relative surface stability and crystal morphology in vacuum conditions. The created surfaces were separated by a vacuum space of 10 Å along the surface normal.
According to Equation (2), the calculated surface energies are 0.042, 0.336, 0.338, 0.544, 0.796, and 0.781 J/m 2 for the (001), (100), 11 0 , (110), (101) and 11 1 surfaces, respectively. We note that the (100) and 11 0 cleavages expose the same surface atoms hence they have very similar values for the surface energies. The (001) surface which results from cleaving the crystal through the weakly interacting hydroxyl layers ( Figure  1b) has the lowest surface energy. No chemical bonds were broken during the cleavage of the (001) surface and the calcium atoms remain fully coordinated. The 0.042 J/m 2 value obtained in our calculations compares well with the 0.051 J/m 2 value obtained from other DFT simulations [61], 0.073 J/m 2 by molecular mechanics calculations [62] and 0.107 J/m 2 from theoretical studies [63]. Indeed, experiments predict a perfect cleavage along the (001) plane [64] for a naturally occurring portlandite mineral. For the second most stable surface; the (100), the surface calcium atoms are five-fold coordinated, while they are fourfold coordinated for the (110) surface. It is important to note that the trend of the relative surface energies follows the degree of under-coordination of the calcium atoms at the surface and the presence of kinks and steps. The position of the cleavage planes for the low surface energy planes is illustrated in Figure 1b,c. The equilibrium morphology of portlandite based on the calculated surface energies and Wulff construction is a hexagonal platelet (Figure 1d). The basal (001) surface is dominant with the (100) surfaces constituting the thin edges. The platelet morphology of portlandite is typical for cementitious materials hydrated in the presence of sulfates [65][66][67] and naturally occurring portlandite minerals [64]. The (001), (100) and (110) faces are chosen for our carbonation analysis due to their high surface stability and variation of surface chemistry.
To develop a clear picture of the carbonation reaction of CO2 with portlandite surfaces, we analyzed two scenarios: portlandite nanopores filled with pure scCO2 and a H2O/CO2 mixture. The portlandite/scCO2 and portlandite/H2O/CO2 atomistic models were built by inserting the CO2 and H2O/CO2 molecules within the 10 Å vacuum space of the (001), (100) and (110) pores. The choice of the pore size ensures that the systems are within the range of few hundred atoms, thus limiting the computational cost of the DFT calculations, while it is already large enough for the simulations to realistically model not only the first coordination sphere of the reacting interfacial fluids molecules, but also their second coordination sphere. The number of molecules is selected such that the densities of the pore solution are 0.78 g/cm 3 and 0.95 g/cm 3 for the pure scCO2 and H2O/CO2 mixture, respectively (Table 1). To ensure that the pore fluids are at supercritical conditions, the simulations were performed at a temperature and pressure values of 323 K and 90 bar, respectively. Consequently, all the systems were first pre-equilibrated by performing classical MD simulations in the isothermal-isobaric (NPT) statistical ensemble for 2 ns to acquire the desired pressure conditions. The classical MD simulations were performed using the LAMMPS [68] package and the atomic interaction parameters were described by the ClayFF [69] force field. The force field parameters for the CO2 molecules were taken from the work by Cygan et. al [70]. The timestep was set to 0.5 fs and the Coulomb interactions are modelled using the PPPM method [71] with a precision of 10 −5 and a cutoff distance of 10 Å for pair interactions. Periodic boundary conditions were applied in all the three dimensions. The final snapshots of the equilibrated trajectories were then taken as input models for the AIMD simulations.

The Reactivity of Pure ScCO2 with Portlandite Surfaces
Due to the absence of water within the portlandite pores, the carbonation reaction of portlandite described in this section corresponds to a solid-state reaction [32,33,40], which thermodynamically proceeds as shown in equation (3). On exposing a portlandite crystal to scCO2, various crystal facets react at different reaction rates. Our results show that the (100) and the (110) surfaces undergo a rapid barrierless reaction to form (bi)carbonate complexes on the surfaces within the first 5 ps of equilibration (see Figure 2a,b for the snapshots of the reacted surfaces). The spontaneous reaction is ascribed to the defects in the local coordination environment of the two surfaces as the reactivity of atoms surrounding a vacant site is reported to increase due to their reduced coordination state. [72]. The carbonation reaction occurs through a concerted C-O bond formation with the surface oxygen atoms (sp3 hybridization of the C atom) and a proton transfer to form an H2O molecule. Some protons hop to the adjacent layers in the bulk region, facilitating formation of H2O molecules within the layers. The water produced from the reaction has been shown to impact the carbonation reaction by facilitating the dissolution of the calcium atoms, thus allowing further carbonation reactions [39,[73][74][75][76]: After a few ps of AIMD simulations, the (100) and the (110) surfaces are saturated with the (bi)carbonate complexes, which occupy the positions formerly occupied by the surface hydroxyl groups. The surface (bi)carbonate complexes form a carbonate layer which acts as a barrier to further intake of CO2. Consequently, the carbonation reaction ceases, in accordance with previous studies [37,73,74,77,78]. There is no penetration of CO2 molecules into the bulk crystal, and thus a surface-only reaction is observed. Therefore, the rate-limiting step of the solid-state carbonation reaction on the (100) and (110) surfaces is interpreted as the formation of a carbonate surface barrier. On the other hand, the (001) basal surface does not react spontaneously with scCO2 like the edge surfaces, and no reaction was observed during the equilibration time. This implies that fully hydroxylated surfaces suppress the formation of carbonate ions, hence, a higher reactivity barrier for CO2 with the (001) surface in comparison to the (100) and (110) surfaces. Indeed, the probability of carbonation of the (001) surface in solid-state form has been reported to be very low under ambient conditions [28] and experimental results have indicated that almost no carbonation occurs in the absence of water [73,79].
To enhance the sampling of the configurational space and estimate the free-energy barrier of the carbonation reaction within the (001) pore, a metadynamics calculation was set up with the coordination number of C atoms to all surface O atoms as the collective variable. As shown in Figure 2c, the metadynamics calculations activate the reaction of CO2 with the (001) surface to form a bicarbonate ion. The free-energy reaction barrier for the reaction is calculated as 11.07 kcal/mol (Figure 2d). This result is in good agreement with previous metadynamics simulations [80,81] and the 11.5 kcal/mol value observed experimentally [82]. The resultant bicarbonate ion formed from the reaction is stable and remains coordinated to the surface. No transformation to carbonate ion was observed during the entire simulation time. This could be attributed to the high energy cost required to dehydroxylate the (001) surface.

Structural and Dynamical Analysis of ScCO2 Within Portlandite Pores
To further analyze the reaction products from the carbonation reactions, the probability distribution of the angle O-C-O was plotted, as shown in Figure 3a. For all the surfaces, most of the CO2 molecules remain unreacted as demonstrated by the high peaks at ~180°. The reacted (bi)carbonate ions produced on carbonation present a peak at ~120°. Additionally, vibrational modes of the formed (bi)carbonate ions were calculated from the power spectra, obtained by Fourier transformation of the velocity-velocity auto correlation function of the CO2 molecules. The vibrational modes of interest here are the modes associated with the carbonate or bicarbonate ions and are highlighted with arrows in Figure 3b. In the region between 800 and 900 cm −1 , carbonates exhibit out-of-plane and in-plane bending absorptions ( ) [83,84]. These bending modes are present in all the three systems at 802 cm −1 , 818 cm −1 and 823 cm −1 for the (001), (100) and (110) surfaces, respectively. The weaker bands referred to as the totally symmetric stretching frequency, ( -, tot sym) are IR active only when the D3h symmetry of the carbonate anion is broken, and it is hydrogen bonded or directly coordinated to the surface or to protons [84]. In our calculations, the symmetry of the formed carbonate anions is expected to degrade from D3h to C1 asymmetries in the local hydrogen bonding environment, and this is demonstrated by the occurrence of thetot sym bands at 1047 cm −1 and 1052 cm −1 for the (100) and (110) surfaces, respectively. vibrational mode is shifted to ~1195 cm −1 -an indication that the adsorbed bicarbonate ions form a bridged structure on the surface [85]. The bridged structure has previously been proposed to best represent the bonding mode of bicarbonate anion adsorbed on aluminum oxide and iron oxide surfaces, and other hydroxylated metal oxides that have similar vibrational spectra [85]. The presence of distinct -, and , vibrational bands associated with the carbonate anion is an indication that the D3h symmetry of the anion is broken by coordination to cations [84]. For the (100) and (110) surfaces, the -, band is situated at ~1332 and 1342 cm −1 and the -, band at ~1480 and 1500 cm -1 . This implies that the formed carbonate anions are strongly bound to the portlandite surfaces through interaction with the surface calcium cations.
The atomic density profiles of the CO2 carbon (C), CO2 oxygen (Oc), hydroxyl oxygen (Oh), hydroxyl hydrogen (H) and calcium (Ca) atoms along the z direction are shown in Figure 4 for the (001), (100) and (110) pores. For the (001) pore, both the oxygen and hydrogen atoms of the surface hydroxyls show sharp peaks (Figure 4a), implying that they mostly remain undistorted on interaction with pure scCO2 and their density profile corresponds to their crystallographic positions. Similarly, the CO2 molecules are highly ordered as indicated by the formation of one distinct peak near the surface and a small peak in the middle of the pore. The first CO2 molecule peak is located at ~4 Å from the surface (the surface is defined by the average position of the surface calcium atoms) and both C and Oc atomic density peaks are coinciding, which suggests that the CO2 molecules are lying almost parallel to the surface. The small peaks of the C and Oc atoms, whose positions coincide with the hydroxyl groups, represent the reacted bicarbonate ion.
On the other hand, for the (100) and (110) pores, the atomic density profile of the hydroxyl groups forms peaks (see Figure 4b,c) that are less sharp in comparison to the (001) case. This feature highlights a distorted local structure, in contrast to the sharply defined maxima of the density profiles characteristic of the crystal atoms, shown here for the (001) surface. The C and Oc peaks strongly overlap with the surface Oh and H atom peaks. We attribute the C and Oc peaks on the surface region to the reacted (bi)carbonate ions while the peaks beyond ~2 Å from the surface to the unreacted CO2 molecules. Based on the density profiles, we deduce that the (bi)carbonate ions formed by the carbonation reaction do not diffuse away from the surface. Furthermore, no dissolution of the calcium atoms is observed during our simulations. The slight distortion of the calcium atoms positions in the (100) and (110) surfaces is attributed to their rearrangement due to the saturation of the (bi)carbonate ions on these surfaces. The diffusivity of CO2 molecules and CO anions formed after the carbonation reaction is analyzed by plotting the mean square displacement (MSD) as a function of time for all simulated systems (Figure 5a-c). The results show that the CO2 molecules diffuse mainly within the xy plane as indicated by the linear relationship between the MSD and time. The slow diffusivity along the z direction is a typical feature of fluids in nanoconfinement [87][88][89][90][91][92]. Of the three systems considered, CO2 diffusion within the (001) pore is faster in comparison to the (100) and (110) pores. This behavior could be explained by the weak interactions of CO2 molecules with the (001) surface, as well as their location ~4 Å away from the surface, as shown in the density profiles, Figure 4a. The implication of this finding to fluid transport of CO2-rich fluids within cement is that CO2 will diffuse faster within the pores of this nature, without necessarily readily reacting. The MSD of the CO ions within the (100) and (110) systems shown in Figure 5b,c is flat both in xy and z directions. This indicates that the diffusion of the CO ions is very slow, and it can be concluded that these ions are static, at least within the time scales of our simulations.

Density Profiles and H2O Molecular Orientation
The inclusion of water within the pores modifies the interaction of CO2 molecules with portlandite surfaces, as shown by the density profile plots (Figure 6a-c). For the (001) surface, the density profile of the water oxygen (Ow) and hydrogen (Hw) atoms show the first peak at 3.76 Å and 2.68 Å from the surface, respectively. The orientation of the dipole vector of H2O molecules in the first layer is ~170° with respect to the surface normal (Figure 6d). Given their location and orientation, it implies that the water molecules in the first layer interact directly with the (001) surface by donating a hydrogen bond to the oxygen atoms of the portlandite surface, in agreement with previous atomistic simulations of portlandite/water interface [62,93]. CO2 molecules occupy the central region of the pore at ~5 Å from the surface (Figure 6a), hence, they do not alter the local coordination environment of H2O molecules at the portlandite (001) surface. The small peaks displayed by the C atoms density profile indicate that the CO2 molecules have a well-defined ordering on the H2O/CO2 interface. There is only a negligible distortion of the surface hydroxyl groups and calcium atoms as demonstrated by the highly ordered sharp peaks of Oh, H, and Ca atomic density profiles.
For the (100) surface, the Ow and Hw atomic density profiles of the H2O molecules show sharp peaks located between the surface oxygen and hydrogen atoms (Figure 6b). These water molecules occupy the surface sites created by the missing hydroxyl groups due to the five-fold coordination of the calcium atoms on the (100) surface. As shown in Figure 6e, the corresponding H2O molecules lie flat, parallel to the surface, their dipole vectors forming ~90° with the surface normal. Since these water molecules occupy surface sites and thus physically bound to the surface, they are not expected to hinder the interaction of the CO2 molecules with the (100) surface. The second peak of the water oxygen atoms is located at ~4 Å, coinciding with that of the C and Oc atoms. Therefore, there is no formation of a well-defined water layer between the (100) surface and the CO2 molecules. Consequently, this could imply a more facile reactivity of CO2 with the (100) portlandite surface. The CO2 molecules occupy the central region of the pore and, similar to the (001) pore, their density profile shows peaks suggesting a definite ordering. The surface atom density profiles (Oh, H, and Ca) exhibit well-ordered sharp peaks ( Figure 6b) and it is concluded that the H2O/CO2 mixture within the (100) pore does not significantly affect the surface structure. In contrast, the density profile calculated for the (110) pore demonstrates a highly distorted surface, indicated by the surface hydroxyl group density profile (Figure 6c). Some of the hydroxyl groups show peaks in the middle of the pore. The migration of the surface hydrogen atoms up to the center of the pores is understood to follow a Grotthusstype mechanism, through proton transfers mediated by the H2O molecules. Water molecules form two distinct peaks at z = 0 Å and z = 2.05 Å with a shoulder at z = 1.32 Å. We interpret the water Ow atom peak at z = 0 Å and the shoulder at z = 1.32 Å to be formed by H2O molecules coordinating directly with the undercoordinated surface calcium atoms. The orientation of the H2O dipole vector with respect to the surface normal of the water molecules contributing to the first peak is ~50°, while for the shoulder and the second peak it is ~70° (Figure 6f). In the middle of the pore, the density profile of the H2O molecules displays a uniform distribution, although it is somewhat less than the bulk water density of ~1 g/cm 3 . The CO2 molecules form a peak at the central region of the pore and, similar to the (001) and (100) pores, no dissolution of the calcium atoms is observed.

Free-Energy Calculations
To quantify the reactivity of CO2 with the portlandite surfaces in the presence of water, metadynamics calculations were performed for all interfaces considered so far. Unlike the case of pure scCO2 discussed in section 3.1, we did not observe any barrierless reaction of CO2 with the portlandite surfaces since the H2O molecules hydrate the surfaces and hinder direct interaction of CO2 molecules with the surfaces. Figure 7a shows the evolution of the free energy as function of the carbon atom coordination number to the oxygen atoms. This collective variable takes into account the oxygen atoms of both H2O molecules and surface hydroxyl groups. The calculated free energies of reaction are 16.43 kcal/mol, 10.08 kcal/mol, and 18.13 kcal/mol for the (001), (100) and (110) surfaces, respectively. The carbonation reaction product in all cases is a bicarbonate ion, illustrated in the snapshots in Figure 7b. The values obtained by our approach are in good agreement with the reaction energy barriers predicted for the formation of bicarbonate complexes [80][81][82] and carbonic acid [82,[94][95][96][97]. The free-energy barrier of 18.13 kcal/mol for the reaction on the (110) surface is the highest for the surfaces investigated here. The carbonation reaction on the (110) surface is preceded by many proton transfers among the water molecules and the surface hydroxyl groups. This is demonstrated by the migration of hydroxyl groups up to the middle of the pores shown by the density profile plots reported in the section 3.2.1. It directly involves concerted interaction of two H2O molecules, a hydroxyl group and a CO2 molecule (illustrated as the reaction path 2 in Figure 8). The reaction barrier is in good agreement with the 17.4 kcal/mol [95] and 19.3 kcal/mol [82] values estimated for the formation of carbonic acid under similar conditions. For both the (001) and (100) surfaces, the carbonation reaction involved one H2O molecule, a surface OH group and a CO2 molecule (see the reaction path 1 in Figure 8). The H2O molecule donates a proton to the surface hydroxyl group and the resultant hydroxyl ion in turn interacts with a CO2 molecule to form a bicarbonate ion. The carbonation reaction barrier within the (100) surface is, however, noticeably lower and is within the range of values obtained for the reaction of CO2 with the (001) surface in the absence of water, as discussed in section 3.1. We relate the variation of the free-energy barriers in the formation of the bicarbonate complex to the interfacial structure of water [98,99] which is, in turn, dictated by the surface chemistry and the degree of nanoconfinement [87,88]. A threshold H2O film thickness is required for carbonate precipitation, and it has been recently shown [100] that carbonate precipitation becomes predominant at about 1.5 monolayers of water, below which the reaction products are limited to (bi)carbonate surface complexes. The arrangement of the water molecules in the vicinity of the three surfaces is distinct (see the density profiles and H2O orientation plots in Figure 6), and this explains the variation in the carbonation reaction mechanisms. Notably, within the (001) and (110) pores, H2O molecules form a distinct layer between the CO2 molecules and the surfaces in contrast to the (100) pore where most of the water molecules occupy surface sites created by the missing hydroxyl groups due to the five-fold coordination of the surface Ca atoms, with only a small percentage occupying the rest of the pore space (see the density profiles and H2O orientation plots, Figure 4b,e). Thus, there is no formation of a water layer between the CO2 and the (100) surface. Moreover, the adsorbed H2O molecules are strongly surfaceassociated and are not sufficiently free to reorient and to participate in the hydration of CO2 [98]. The low reaction barrier for the carbonation reaction on the (100) surface is hence associated with the chemisorbed H2O molecules and the absence of a full water layer between the CO2 and the surface.

Analysis of the Dynamical Properties
Analysis of the diffusional properties of H2O and CO2 molecules was performed by plotting their MSD as a function of time (Figure 9a-c). CO2 molecules diffuse faster than H2O molecules, noticeably within the xy plane. In addition to the small size of the pores, the presence of H2O molecules within the pores is shown to further attenuate the CO2 molecule diffusion in comparison to the pure scCO2 systems discussed in section 3.1.1. By contact with CO2-rich fluids [101], the CO2 diffusivity plays here a key role in controlling the degree of cement reactivity, and these results indicate that the extent of CO2 molecular transport within cement is limited for fully water-filled pores like the inter/intralayer and small gel pores of the C-S-H phase. The diffusion of H2O molecules is slow in general (Figure 9a-c), which we relate to their strong interaction with the portlandite surfaces. Due to their location within the pores, water molecules interact directly with the portlandite surfaces ( Figure 6 and Figure  S3 in the supporting information). The parallel diffusion of H2O molecules within the (001) pore is, however, relatively fast in comparison to the (100) and (110) pores, and this can be associated with the relatively weak water-surface interactions [62]. The negligible parallel diffusion of H2O molecules within the (100) and the (110) pores (Figure 9b,c) is attributed to the strong interaction of the water molecules with those surfaces since some of the H2O molecules occupy surface sites. The diffusion in the z direction is negligible in all cases for both H2O and CO2 molecules, ascribed to nanoconfinement in the z direction.

Conclusions
In this work, we investigated the reaction mechanisms of dry scCO2 and H2O/CO2 mixtures with portlandite surfaces, relevant for understanding the interaction of CO2 with cement in subsurface operations. Our simulations demonstrate that dry scCO2 undergoes a rapid barrierless carbonation reaction with the edge surfaces of portlandite crystals. However, the carbonation reaction soon ceases due to the deposition of (bi)carbonate surface complexes which form a carbonate layer. We therefore propose the reaction of scCO2 with portlandite to be primarily interface-controlled, i.e., it depends on the exposed surface and then secondly to be diffusion-controlled in that the rate-limiting step is the formation of the surface carbonate barrier which controls the diffusion of CO2 through the layer. The implications of the formation of the carbonate passivating layer are manifold: while it has positive effects by blocking the cracks which could act as CO2 leakage pathways and attenuate cement degradation, it would have a detrimental effect in mineral carbonation CCS operations where CO2 is permanently stored through reaction with host rocks. It is important to emphasize that this conclusion is limited to the relatively small pore sizes studied in this work and we might expect that increasing the pore sizes could favor diffusion of the CO2/(bi)carbonate complexes thus enhancing the CO2 reactivity.
The presence of water alters the interaction of CO2 with portlandite surfaces. Our results show that H2O molecules coordinate directly to the portlandite surfaces forming well-structured water layers which CO2 molecules cannot easily penetrate. It is deduced that the rate-limiting step in the carbonation reaction of portlandite and the H2O/CO2 mixture is the structure and the water content within the portlandite pores. As such, CO2 reactivity for pores with highly structured water-surface layers (with no bulk-like water) is expected to be limited due to the attenuated inward diffusion of the CO2 molecules.
Dissolution of calcium atoms can enhance the carbonation reaction by facilitating their interaction with the dissolved (bi)carbonate/carbonic acid to form calcite in the presence of water [31]. Dissolution-precipitation carbonation reaction mechanism is expected to be dominant in the presence of water. However, no dissolution of calcium atoms was observed in our simulations. This is attributed to the small size of the atomistic models and short simulation times permitted by the AIMD approach. To overcome this drawback, future work should involve the use of a reactive molecular dynamics approach using an empirical reactive force field, such as ReaxFF [102], which would allow atomistically model relatively large systems over longer time scale. This approach, however, is more strongly dependent on the choice and accuracy of the force field.
Supplementary Materials: The following are available online at www.mdpi.com/2075-163X/11/5/509/s1, Figure S1: (a) The variation of the total energy with the plane-wave basis cutoff energy. (b) Energy change relative to the previous step. The value of 400 Ry is selected since there is negligible change in energy beyond this value; Figure S2: Variation of the Carbon atom coordination number (CV) with time for the metadynamics simulations of the portlandite/CO2/H2O systems; Figure S3

Data Availability Statement:
The data presented in this study are available within the article.