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

Abstract

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₂ 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₂ with the basal and edge surfaces of a model hydrated cement phase—portlandite—in dry scCO₂ and water-rich conditions. Our simulations show that in dry scCO₂ conditions, the undercoordinated edge surfaces of portlandite experience a fast barrierless reaction with CO₂, 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₂ conditions is the formation of the surface carbonate barrier which controls the diffusion of CO₂ through the layer. The presence of water hinders direct interaction of CO₂ with portlandite as H₂O molecules form well-structured surface layers. In the water-rich environment, CO₂ undergoes a concerted reaction with H₂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.

1. 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 CO₂, it is important to understand geochemical interactions between CO₂ and cement at relevant subsurface operation conditions.

The carbonation reactions of the injected CO₂ with wellbore cement pose serious well integrity issues because of the acid attack on cement [16]. CO₂-saturated fluids are thermodynamically incompatible with cement, since cement is highly alkaline (pH > 12.5) and not in chemical equilibrium with the acidic CO₂-rich fluids (pH < 6) [17]. Depending on the subsurface conditions (e.g., temperature, pressure, salinity), studies report different degrees of cement degradation due to CO₂-rich fluids [11,13,18,19,20,21]. For instance, in dry supercritical CO₂ (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)₂) 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 CO₂ with the portlandite phase. Portlandite is a hydrous mineral which is very reactive with CO₂, making it a good model for understanding carbonation reaction mechanisms in cement and in other hydrous minerals with similar structure (e.g., Mg(OH)₂). 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 (CaSiO₃). In contrast, Vance et. al [30] suggests that so long as suitable reaction conditions are provisioned, the direct carbonation of Ca(OH)₂ 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 scCO₂ and H₂O/CO₂ 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 CO₂–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 CO₂ with the portlandite surfaces. We show that in dry scCO₂ 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 CO₂ with portlandite surface, with the free-energy barrier being governed by the thickness of the interfacial water layer.

2. Methods

2.1. Computational Details

To model the reaction mechanisms of CO₂ 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 CO₂ 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⁻⁶ 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 CO₂ molecules with the portlandite surfaces in pure supercritical state and in a CO₂/H₂O mixture, we chose the carbon atom (for all the CO₂ molecules) coordination number (CN) to oxygen atoms of both the surface hydroxyls (O_h) and of the neighboring H₂O molecules (O_w) 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, C_ij(r_ij):

$$C_{ij}\left(r_{ij}\right)=\frac{1-{\left(\frac{r_{ij}-d_0}{r_0}\right)}^n}{1-{\left(\frac{r_{ij}-d_0}{r_0}\right)}^m}$$

(1)

where r_ij is the instantaneous distance between atoms i and j, r₀ is the distance beyond which the bond is considered to be broken and d₀ is the central value of the switching function (i.e., C_ij(d₀) = 1). m and n are parameters of the switching function chosen in such a way that C_ij(r_ij) tends to zero beyond r₀. 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 DFT-based 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]:

$$F =-\frac{1}{\beta }ln\left({\displaystyle \int }ds e^{-\beta F\left(s\right)}\right)$$

(2)

where $\beta ={(k_{B}T)}^{-1}$, $k_B$ 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:

$$\gamma =\frac{1}{2A}\left(E_{\mathrm{slab}}-n{E}_{\mathrm{bulk}}\right)$$

(3)

where E_slab is the total energy of the slab and E_bulk 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].

2.2. Atomistic Models

Portlandite, Ca(OH)₂, is a layered mineral, with the layers consisting of distorted edge-sharing CaO₆ 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), $\left(1\overline{1}0\right)$, (110), (101) and $\left(1\overline{1}\overline{1}\right)$, 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² for the (001), (100), $\left(1\overline{1}0\right)$, (110), (101) and $\left(1\overline{1}\overline{1}\right)$ surfaces, respectively. We note that the (100) and $\left(1\overline{1}0\right)$ 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² value obtained in our calculations compares well with the 0.051 J/m² value obtained from other DFT simulations [61], 0.073 J/m² by molecular mechanics calculations [62] and 0.107 J/m² 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 four-fold 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 CO₂ with portlandite surfaces, we analyzed two scenarios: portlandite nanopores filled with pure scCO₂ and a H₂O/CO₂ mixture. The portlandite/scCO2 and portlandite/H₂O/CO₂ atomistic models were built by inserting the CO₂ and H₂O/CO₂ 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³ and 0.95 g/cm³ for the pure scCO₂ and H₂O/CO₂ mixture, respectively (

Table 1. The number of CO₂ and H₂O molecules within the portlandite pores cleaved along the (001), (100) and (110) planes for both the pure scCO₂ and H₂O/CO₂ mixture systems.

System	Number of CO2 and H2O Molecules within the (001), (100) and (110) Pores
	(001)		(100)		(110)
Slab size (Å3)	4411.49		4478.17		8188.60
	CO2	H2O	CO2	H2O	CO2	H2O
Pure scCO2	24	-	24	-	44	-
H2O/CO2 mixture	12	38	12	38	12	88

). 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 CO₂ 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⁻⁵ 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.

3. Results and Discussion

3.1. The Reactivity of Pure ScCO₂ 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 scCO₂, 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 H₂O molecule. Some protons hop to the adjacent layers in the bulk region, facilitating formation of H₂O 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]:

$$\mathrm{Ca}{\left(\mathrm{OH}\right)}_2+{\mathrm{CO}}_2⟶{\mathrm{CaCO}}_3+{\mathrm{H}}_2\mathrm{O}$$

(4)

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 CO₂. Consequently, the carbonation reaction ceases, in accordance with previous studies [37,73,74,77,78]. There is no penetration of CO₂ 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 scCO₂ 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 CO₂ 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 CO₂ with the (001) surface to form a bicarbonate ion. The free-energy reaction barrier for the reaction

$${\mathrm{CO}}_2+{\mathrm{OH}}^{-}⟶{\mathrm{HCO}}_3^{-}$$

(5)

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.

3.1.1. Structural and Dynamical Analysis of ScCO₂ 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 CO₂ 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 CO₂ 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⁻¹, carbonates exhibit out-of-plane and in-plane bending absorptions (${\nu }_{{\mathrm{CO}}_3}$) [83,84]. These bending modes are present in all the three systems at 802 cm⁻¹, 818 cm⁻¹ and 823 cm⁻¹ for the (001), (100) and (110) surfaces, respectively. The weaker bands referred to as the totally symmetric stretching frequency, (${\nu }_{\mathrm{C}–\mathrm{O}}$, tot sym) are IR active only when the D_3h 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 D_3h to C₁ asymmetries in the local hydrogen bonding environment, and this is demonstrated by the occurrence of the ${\nu }_{\mathrm{C}–\mathrm{O}}$ tot sym bands at 1047 cm⁻¹ and 1052 cm⁻¹ for the (100) and (110) surfaces, respectively.

The C–O–H bending vibration (${\delta }_{\mathrm{C}–\mathrm{O}–\mathrm{H}}$), symmetric C-O stretch (${\nu }_{\mathrm{C}–\mathrm{O},\mathrm{sym}}$) and asymmetric C–O stretch (${\nu }_{\mathrm{C}-\mathrm{O}, \mathrm{asym}}$) bands experimentally assigned to the formation of a bicarbonate ion on hydroxylated surfaces at 1220, 1410, and 1622 cm⁻¹ [85,86], respectively, are exhibited at 1218, 1413, and 1640 cm⁻¹ for the (001) surfaces. The calculated frequencies agree with the experimental values within 20 cm⁻¹. For the (100) and (110) surfaces, the ${\delta }_{\mathrm{C}–\mathrm{O}–\mathrm{H}}$ vibrational mode is shifted to ~1195 cm⁻¹—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 ${\nu }_{\mathrm{C}–\mathrm{O},\mathrm{sym}}$ and ${\nu }_{\mathrm{C}-\mathrm{O}, \mathrm{asym}}$ vibrational bands associated with the carbonate anion is an indication that the D_3h symmetry of the anion is broken by coordination to cations [84]. For the (100) and (110) surfaces, the ${\nu }_{\mathrm{C}–\mathrm{O},\mathrm{sym}}$ band is situated at ~1332 and 1342 cm⁻¹ and the ${\nu }_{\mathrm{C}–\mathrm{O}, \mathrm{asym}}$ band at ~1480 and 1500 cm⁻¹. 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 CO₂ carbon (C), CO₂ oxygen (O_c), hydroxyl oxygen (O_h), 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 scCO₂ and their density profile corresponds to their crystallographic positions. Similarly, the CO₂ 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 CO₂ 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 O_c atomic density peaks are coinciding, which suggests that the CO₂ molecules are lying almost parallel to the surface. The small peaks of the C and O_c 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 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 ions while the peaks beyond ~2 Å from the surface to the unreacted CO₂ 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 CO₂ molecules and ${\mathrm{CO}}_3^{2-}$ 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 CO₂ 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, CO₂ 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 CO₂ 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 CO₂-rich fluids within cement is that CO₂ will diffuse faster within the pores of this nature, without necessarily readily reacting. The MSD of the ${\mathrm{CO}}_3^{2-}$ 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 ${\mathrm{CO}}_3^{2-}$ ions is very slow, and it can be concluded that these ions are static, at least within the time scales of our simulations.

3.2. The Reactivity of H₂O/CO₂ Mixture with Portlandite Surfaces

3.2.1. Density Profiles and H₂O Molecular Orientation

The inclusion of water within the pores modifies the interaction of CO₂ 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 (O_w) and hydrogen (H_w) atoms show the first peak at 3.76 Å and 2.68 Å from the surface, respectively. The orientation of the dipole vector of H₂O 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]. CO₂ 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 H₂O molecules at the portlandite (001) surface. The small peaks displayed by the C atoms density profile indicate that the CO₂ molecules have a well-defined ordering on the H₂O/CO₂ interface. There is only a negligible distortion of the surface hydroxyl groups and calcium atoms as demonstrated by the highly ordered sharp peaks of O_h, H, and Ca atomic density profiles.

For the (100) surface, the O_w and H_w atomic density profiles of the H₂O 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 H₂O 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 CO₂ molecules with the (100) surface. The second peak of the water oxygen atoms is located at ~4 Å, coinciding with that of the C and O_c atoms. Therefore, there is no formation of a well-defined water layer between the (100) surface and the CO₂ molecules. Consequently, this could imply a more facile reactivity of CO₂ with the (100) portlandite surface. The CO₂ 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 (O_h, H, and Ca) exhibit well-ordered sharp peaks (Figure 6b) and it is concluded that the H₂O/CO₂ 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 Grotthuss-type mechanism, through proton transfers mediated by the H₂O 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 O_w atom peak at z = 0 Å and the shoulder at z = 1.32 Å to be formed by H₂O molecules coordinating directly with the undercoordinated surface calcium atoms. The orientation of the H₂O 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 H₂O molecules displays a uniform distribution, although it is somewhat less than the bulk water density of ~1 g/cm³. The CO₂ 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.

3.2.2. Free-Energy Calculations

To quantify the reactivity of CO₂ with the portlandite surfaces in the presence of water, metadynamics calculations were performed for all interfaces considered so far. Unlike the case of pure scCO₂ discussed in Section 3.1, we did not observe any barrierless reaction of CO₂ with the portlandite surfaces since the H₂O molecules hydrate the surfaces and hinder direct interaction of CO₂ 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 H₂O 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 H₂O molecules, a hydroxyl group and a CO₂ 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 H₂O molecule, a surface OH group and a CO₂ molecule (see the reaction path 1 in Figure 8). The H₂O molecule donates a proton to the surface hydroxyl group and the resultant hydroxyl ion in turn interacts with a CO₂ 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 CO₂ 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 H₂O 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 H₂O orientation plots in Figure 6), and this explains the variation in the carbonation reaction mechanisms. Notably, within the (001) and (110) pores, H₂O molecules form a distinct layer between the CO₂ 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 H₂O orientation plots, Figure 4b,e). Thus, there is no formation of a water layer between the CO₂ and the (100) surface. Moreover, the adsorbed H₂O molecules are strongly surface-associated and are not sufficiently free to reorient and to participate in the hydration of CO₂ [98]. The low reaction barrier for the carbonation reaction on the (100) surface is hence associated with the chemisorbed H₂O molecules and the absence of a full water layer between the CO₂ and the surface.

3.2.3. Analysis of the Dynamical Properties

Analysis of the diffusional properties of H₂O and CO₂ molecules was performed by plotting their MSD as a function of time (Figure 9a–c). CO₂ molecules diffuse faster than H₂O molecules, noticeably within the xy plane. In addition to the small size of the pores, the presence of H₂O molecules within the pores is shown to further attenuate the CO₂ molecule diffusion in comparison to the pure scCO₂ systems discussed in Section 3.1.1. By contact with CO₂-rich fluids [101], the CO₂ diffusivity plays here a key role in controlling the degree of cement reactivity, and these results indicate that the extent of CO₂ 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 H₂O 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 H₂O 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 H₂O 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 H₂O molecules occupy surface sites. The diffusion in the z direction is negligible in all cases for both H₂O and CO₂ molecules, ascribed to nanoconfinement in the z direction.

4. Conclusions

In this work, we investigated the reaction mechanisms of dry scCO₂ and H₂O/CO₂ mixtures with portlandite surfaces, relevant for understanding the interaction of CO₂ with cement in subsurface operations. Our simulations demonstrate that dry scCO₂ 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 scCO₂ 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 CO₂ 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 CO₂ leakage pathways and attenuate cement degradation, it would have a detrimental effect in mineral carbonation CCS operations where CO₂ 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 CO₂/(bi)carbonate complexes thus enhancing the CO₂ reactivity.

The presence of water alters the interaction of CO₂ with portlandite surfaces. Our results show that H₂O molecules coordinate directly to the portlandite surfaces forming well-structured water layers which CO₂ molecules cannot easily penetrate. It is deduced that the rate-limiting step in the carbonation reaction of portlandite and the H₂O/CO₂ mixture is the structure and the water content within the portlandite pores. As such, CO₂ 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 CO₂ 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.