Article pubs.acs.org/Langmuir
Molecular Dynamics Simulations of CO2/Water/Quartz Interfacial Properties: Impact of CO2 Dissolution in Water Gina Javanbakht, Mohammad Sedghi, William Welch, and Lamia Goual* Department of Chemical and Petroleum Engineering, University of Wyoming, 1000 E. University Avenue, Laramie, Wyoming 82071, United States ABSTRACT: The safe trapping of carbon dioxide (CO2) in deep saline aquifers is one of the major concerns of CO2 sequestration. The amount of capillary trapping is dominated by the capillary pressure of water and CO2 inside the reservoir, which in turn is controlled by the interfacial tension (IFT) and the contact angle (CA) of CO2/water/rock systems. The measurement of IFT and CA could be very challenging at reservoir conditions, especially in the presence of toxic cocontaminants. Thus, the ability to accurately predict these interfacial properties at reservoir conditions is very advantageous. Although the majority of existing molecular dynamics (MD) studies of CO2/water/mineral systems were able to capture the trends in IFT and CA variations with pressure and temperature, their predictions often deviated from experimental data, possibly due to erroneous models and/or overlooked chemical reactions. The objective of this study was to improve the MD predictions of IFT and CA of CO2/ water/quartz systems at various pressure and temperature conditions by (i) considering the chemical reactions between CO2 and water and (ii) using a new molecular model for α-quartz surface. The results showed that the presence of carbonic acid at the CO2/water interface improved the predictions of IFT, especially at low temperature and high pressure where more CO2 dissolution occurs. On the other hand, the effect on CA was minor. The slight decrease in CA observed across the pressure range investigated could be attributed to an increase in the total number of H-bonds between fluid molecules and quartz surface. that among molecular models for water, TIP4P/200514 could produce the closest match to experimental data. Other researchers have reached the same conclusion by using a dummy particle to carry the negative charge of oxygen atom in various water models.14,15 Older water models (e.g., TIP3P,16 SPC,17 TIP4P,14 and SPC/E18) were not as accurate as TIP4P/ 2005 in predicting the behavior of water molecules during MD simulations. Several models were also suggested for CO2 such as MSM, EPM2, TraPPE, Errington, etc.5 In order to find the best model for CO2/water systems at high pressure (P) and temperature (T), Liu et al. examined several existing models for water (SPC, TIP4P, TIP4P2005) and CO2 (EPM2, TraPPE) and found that TraPPE/TIP4P2005 and EPM2/SPC combinations provided the most accurate predictions of CO2 solubility and IFT variations with P and T.6,16 Kvamme et al. used the SPC/E model for water and a simple EMP model for CO2 to study the IFT between CO2 and water at 278−335 K temperatures and 0.1−20 MPa pressures. Their results were in good agreement with experimental data and revealed a decrease in IFT with increasing pressure.7 More recently, Li et al. studied the IFT of CO2/brine systems at higher P and T using TIP4P and SPC/E models for water and EPM2 model for CO2.8 They showed that the SPC/E water model gives a better match between experimental and simulation data. However, the simulation errors were higher at elevated pressures. In addition to IFT, the CA between two fluid phases and a solid phase controls the capillary pressure of the system. Using
1. INTRODUCTION The successful sequestration of carbon dioxide (CO2) in deep geological formations requires a minimum leakage of CO2 through the cap rock. The amount of leakage could be estimated from the threshold capillary pressure, which in turn is affected by the interfacial tension (IFT) between CO2 and the present fluids, as well as the contact angle (CA) between fluids and the solid surface, according to the following equation1 Pc = [2γα , β cos(θα , β , δ)]/r
(1)
where γα,β is the interfacial tension between phases α and β and θα,β,δ is the contact angle between phases α and β and solid surface δ. Although this relationship applies to perfect capillary tubes, it is more complicated in pore spaces as pore geometry affects the capillary pressure.2,3 In addition to CO2 sequestration, the prediction of IFT in CO2 enhanced oil recovery (EOR) processes is important since IFT (or γ) controls the spreading coefficient of one fluid over the others4 according to S = γgas,brine − γoil,gas − γoil,brine
(2)
The measurement of IFT and CA at reservoir conditions is usually very challenging, especially in the presence of cocontaminants such as SO2 and NOx. Therefore, the ability to accurately predict these properties at reservoir conditions can be very useful. Several research groups have used molecular dynamics (MD) modeling to estimate the IFT and CA of CO2/ water/quartz systems under various conditions.5−12 In 2007, Vega and Miguel calculated the surface tension of water in order to compare different water models.13 Their study showed © XXXX American Chemical Society
Received: February 3, 2015 Revised: April 14, 2015
A
DOI: 10.1021/acs.langmuir.5b00445 Langmuir XXXX, XXX, XXX−XXX
Article
Langmuir the same models, Liu et al. investigated the effect of pressure on the wetting behavior of CO2/water/quartz systems by considering both β-quartz and quartz with Silanol groups.9 They found that the CA increased as pressure increased for both surfaces. However, unlike hydrophilic quartz, water droplets on hydrophobic quartz detached from the surface by increasing CO2 pressure. More recently, Iglauer et al. studied the IFT and CA of CO2/water/α-quartz systems using the EMP2/TIP4P2005 model combination up to 20 MPa pressure.10 The partial charges of quartz were 2.4 for Si and −1.2 for oxygen atoms. Their results showed an increase in CA with increasing pressure. IFT, on the other hand, decreased with increasing pressure and increased with increasing temperature. The model used in that study was accurate in predicting IFT versus pressure below the saturation pressure CO2 at 300 K (i.e., subcritical CO2) but overpredicted IFT by 20% at higher pressures (i.e., supercritical CO2). This study was continued by McCaughan et al., who examined the effect of divalent salts (such as CaCl2) on the wettability of CO2/H2S/ water/α-quartz systems during subsurface gas injection.11 They reported that salt addition did not affect the CA. On the other hand, H2S interacted more strongly with quartz than CO2, resulting in a higher CA under the same P and T conditions. Although previous MD studies could capture the trends in IFT and CA variations with P and T, the simulation results often deviated from experimental data, especially at high P and T conditions. A plausible explanation for this discrepancy is that, by default, MD simulations do not explicitly take into account the chemical reactions that may occur at fluid/fluid and fluid/solid interfaces. Moreover, the molecular models for the solid surface used in previous CA predictions might have been inaccurate. The goal of this study was to improve the predictions of IFT and CA of CO2/water/α-quartz systems under a wide range of T and P (from subcritical to supercritical CO2) by considering chemical reactions between CO2 and water and using a novel molecular model for the α-quartz surface.
HCO−3 ⇌ H+ + CO32 −
(6)
When combined with eq 4 the expression of K1 becomes a H+a HCO−3 K1 = K CO2PCO2 Since aH = a +
− HCO3
(8) +
due to charge balance, the activity of H reduces
to
a H+ =
K1K CO2PCO2
(9)
By increasing CO2 pressure from 3.4 to 12 MPa, the pH (i.e., −log aH+) decreases from 3.15 to 2.80 at 308 K and from 3.26 to 3.00 at 333 K, due to the presence of carbonic acid. Under these low pH conditions the dissociation of H2CO3 is negligible. For example at a pH of 2.88, the concentration of H3O+ is 10−2.88 M and is incredibly small in comparison to that of water (55 M). The stoichiometric ratio of hydronium/water is 2.4 × 10−5:1. In our simulations, which contain 4178 water molecules, this would be a small fraction of a hydronium ion. However, there is a considerable amount of solvated CO2/ carbonic acid at the interface, which can have a sizable impact on the IFT between CO2 and water. In our simulations, we used H2CO3 to simulate the solvation of CO2, as suggested by ref 20. At each P and T condition, the concentration of carbonic acid was calculated, and the corresponding number of H2CO3 molecules was added to the simulation box. 2.2. Intermolecular Forces and MD Simulations. The intermolecular potential forces used in this work are divided into van der Waals interactions and electrostatic interactions. The van der Waals interactions were calculated from Lennard-Jones potential
⎡⎛ σij ⎞12 ⎛ σij ⎞6 ⎤ Uij = 4εij⎢⎜ ⎟ − ⎜ ⎟ ⎥ ⎝r ⎠⎦ ⎣⎝ r ⎠
2.1. Chemical Reactions. Reactions involving carbonate systems control the pH of most natural waters. When carbon dioxide comes in contact with water, it starts to dissolve in the water until an equilibrium state is reached. The concentration of dissolved CO2 depends on the P and T conditions. The dissolution of CO2 in water and subsequent reaction is presented by the following20
(10)
where εij is the interaction strength between atoms i and j at a distance r, σij is the effective diameter of atoms i and j where σij = ((σii + σjj))/2 and εij = ((εiiεij)1/2), σii is the diameter of atom i, and εij is the interaction strength between atoms i and j. Electrostatic interactions were calculated from Coulomb equation
(3)
The right-hand side of the reaction above includes both carbonic acid and solvated CO2 in water. In fact the amount of solvated CO2 is much greater than carbonic acid. In other words, both carbonic acid and solvated CO2 are represented by H2CO3.20 The activity of H2CO3 can be obtained from the equilibrium constant (KCO2) and the pressure of carbon dioxide in the CO2-rich phase (PCO2) a H2CO3 = K CO2PCO2
(5)
The equilibrium constants K1 and K2 of reactions 5 and 6 are tabulated in ref 20 at various temperatures. At 308 K, K1 = 10−6.33 and K2 = 10−10.29 whereas at 333 K, K1 = 10−6.29 and K2 = 10−10.14. The constant K2 is very small compared to KCO2 and K1 and can therefore be neglected. Based on reaction 5, K1 can be written as a H+a HCO−3 K1 = a H2CO3 (7)
2. METHODS
CO2 (g) + H 2O ⇌ H 2CO3(aq)
H 2CO3 ⇌ H+ + HCO−3
Fij⃗ =
Ke*Q i*Q jriĵ rij2
(11)
where Ke is Coulomb’s constant, Qα is the absolute value of the partial charge of particle α, and rij is the separation distance between particles i and j. The long-range interactions of electrostatic forces were computed using particle-mesh Ewald summation (PME). The cutoff radios for LJ interactions were set to 1.5 nm, and the periodic boundary condition was applied in all directions. The interfacial tension of the CO2/water system was calculated based on the mechanical route from the diagonal components of the pressure tensor normal to the surface according to
(4)
where KCO2 is equal to 10−1.52 at 308 K and 10−1.78 at 333 K20 and PCO2 varies based on the applied pressure. H2CO3 is neutral and its activity coefficient is equal to unity, meaning that the activity of the carbonic acid in water is equal. Carbonic acid in aqueous media starts to dissociate into bicarbonate and hydrogen as the pH increases. Bicarbonate itself can further dissociate to carbonate and proton at high pH conditions according to20
γ=
⎤ 1 ⎡ 1 LZ ⎢PZZ − (PXX + PYY )⎥ ⎦ 2 ⎣ 2
(12)
where γ is the surface tension, LZ is the system length along the z-axis, and Pxx, Pyy, and PZZ are the components of the pressure tensor.21 B
DOI: 10.1021/acs.langmuir.5b00445 Langmuir XXXX, XXX, XXX−XXX
Article
Langmuir
Figure 1. α-Quartz surface: (a) top view and (b) side view. The solid surface consists of a fully coordinated α-quartz with Si− O−Si bridges on the surface, as seen in Figure 1. The surface density of Silanol groups is 4.9 nm−2, which is in line with experimental data reported for various silica surfaces.22,23 Half of the nonbridge oxygen atoms at the surface were protonated. Half-protonated quartz gave the best contact angle compared to fully protonated and nonprotonated surfaces. The unsaturated oxygen atoms as well as the Silanol groups on the surface were considered flexible with the bonded interactions presented in Table 1. The rest of atoms in the quartz were kept frozen
simulation for 500 ps with Berendsen thermostat and barostat.29 The production simulations were performed using the Nose−Hoover temperature bath and the Parrinello−Rahman pressure coupling30 for 10 ns. For the simulations containing carbonic acid, the concentration of H2CO3 was calculated at each pressure and temperature according to the procedure described in Section 2.1, and a corresponding number of carbonic acid molecules was added to the simulation box. The number of H2CO3 molecules varied between 70 and 260. For contact angle calculations, the simulation box had dimensions of 29 × 29 × 80 nm3 and contained 21 000 SiO4 units in the solid surface. We placed 3548 water molecules on the solid surface and used 29 000 and 58 000 CO2 molecules for low- and high-pressure simulations, respectively. All simulations were performed in NVT ensembles using the Nose−Hoover thermostat.31 The pressure of the CO2 phase was controlled by exerting an external pressure on a piston composed of virtual hard sphere particles in a face-centered cubic (fcc) arrangement. The hard sphere wall was modeled using repulsive interactions (i.e., LJ particles with low attraction potential). To verify the pressure of each simulation, the average CO2 density was calculated and compared with the experimental values, as shown in Table 2. The densities obtained from MD simulations closely match
Table 1. LJ Parameters and Partial Charges of Components Considered in MD Simulations Intermolecular Interactions atoms
ε (kJ/mol)
Si (quartz) 7.700 × 10−6 O (Si−O−Si) 0.6500 O (Si−OH) 0.6500 O (Si−O) 0.6500 H (H2CO3) 0.00 C (H2CO3) 0.43932 O (O−H in H2CO3) 0.71128 O (H2CO3) 0.87864 Intramolecular Interactions in bonds
Kb (kJ/Å2)
σ (Å) 3.3020 3.1655 3.1655 3.1655 0.00 3.7500 3.000 2.960 Silica
q (e) 2.10 −1.050 −0.950 −0.525 0.429 0.88 −0.57 −0.602
Table 2. Measured and Simulated Density of CO2 at Different P and T Conditions density (kg/m3)
r0 (Å)
Si−O O−H angle
2510.40 1627.50 Kθ (kJ/rad2)
1.630 0.950 θ0 (deg)
Si−O−Si
397.488
109.5
pressure (MPa) 3.4 6.2 9.0 11.7
during the simulations. Since the frozen atoms are separately thermostated and are not bonded to the surface atoms, we expect that freezing and restraining their positions will have a negligible impact on the simulation results. ClayFF24 parameters, as reported in ref 25, were used to represent nonbonded interactions for quartz (see Table 1). LJ parameters for carbonic acid atoms were chosen from optimized potentials for liquid simulations-all atoms (OPLS-aa) force field.26,27 To determine the partial charges of atoms in carbonic acid molecules, quantum mechanical (QM) calculations were performed to optimize the structure of H2CO3 molecule using the Gaussian software.28 For QM computations, MP2 level of theory combined with 6-31G(d) basis set was employed, and the charges were calculated by ESP method. The LJ potentials and partial charges of the atoms in the particles used in MD simulations are provided in Table 1. The size of the simulation box for IFT calculation was 5 × 5 × 10 nm3. The box contained 4178 molecules of water and 700 molecules of CO2. Initially, equilibrium conditions were reached by running NPT
3.4 6.2 9.0 11.7
simulation 308 K 79.34 175.24 645.8 761.9 333 K 69.39 135.8 230.02 378.92
exptl19 72.3 175.2 610.4 739.2 63 131.2 233.8 398
the experimental values. At one temperature and pressure condition, a simulation was run for 5 ns, and the contact angle was measured every 500 ps. The contact angle values did not change significantly after 4 ns, and hence all other simulations were performed for 5 ns to ensure the system has reached equilibrium. The contact angles were calculated by the following equation using ImageJ software, where a is the droplet height and b is the droplet base radius,32 as shown in Figure 2.
θ = arcsin C
2ab a 2 + b2
(13) DOI: 10.1021/acs.langmuir.5b00445 Langmuir XXXX, XXX, XXX−XXX
Article
Langmuir
Figure 2. Parameters used in eq 13. Periodic boundary conditions were applied in the X and Y directions. All the simulation were performed using GROMACS 4.6.333 and VMD software34 was used to visualize the simulations results.
3. RESULTS AND DISCUSSION We first ran MD simulations to calculate the surface tension of water and select the best models for water and CO2. We then applied those models to calculate the IFT between CO2 and water and the CA of CO2/water/quartz system (with and without reaction 3). The simulated IFT and CA were compared to experimental data reported by Saraji et al.,19 Sarmadivaleh et al.,35 and Farokhpoor et al.36 These reference data were selected due to their high precision and accuracy.37 3.1. Surface Tension of Water. We tested three different water models (SPC/E, TIP4P, and TIP4P2005) to calculate the surface tension of water. We placed 4050 water molecules in a 4 × 4 × 8 nm3 simulation box and conducted NVT simulations for 1 ns at a temperature of 300 K. Among all models, TIP4P/ 2005 presented the closest surface tensions to experimental data38 (see Table 3). Therefore, we continued our study in the Table 3. Comparison between Measured and Simulated Surface Tension of Water at 300 K Using Different Water Models experiment38
water model surface tension (mN/m)
TIP4P2005
TIP4P
SPC/E
68.3
50
63.8
71.8
next sections using TIP4P/2005 as the water model. We also selected the TraPPE model for CO2 based on its accuracy when combined with TIP4P/2005 water model.6 3.2. Interfacial Tension between CO2 and Water. 3.2.1. Without Reaction. MD simulations were first performed to estimate the IFT between CO2 and water at different pressures (3.4−11.7 MPa) and temperatures (308, 333 K), without considering any reactions between CO2 molecules and water (see dashed lines in Figure 3). The results revealed that IFT decreases with increasing pressure, in good agreement with experimental data measured by axisymmetric drop shape analysis (ADSA) and other methods.19,35 However, like previous MD studies,7,12 the IFT between supercritical CO2 and water was overestimated. The estimated error in the simulations was between 0.25 to 1.1 mN/m. The variation of IFT values with temperature is within the error estimate and is negligible. 3.2.2. With Reaction. Using the same procedure described earlier, NPT simulations of CO2/water IFT were performed by taking into consideration the dissolution of CO2 in water according to reaction 3. The concentration of H2CO3 in water was calculated at each P and T, and the corresponding number of molecules were added to the simulation box (see Section 2.1). A snapshot of CO2, water, and H2CO3 molecules is presented in Figure 4 a where H2CO3 molecules are depicted as red particles in the aqueous phase. This figure shows that some
Figure 3. Effect of pressure on IFT of the CO2/water system at (a) 308 K and (b) 333 K. MD simulations performed with and without reaction between CO2 and water according to reaction 3.
acid molecules accumulate at the CO2/water interface, affecting the IFT of the system. The density profiles of water, CO2, and carbonic acid are shown in Figure 5 at both temperatures and 11.7 MPa. In this figure, the density of H2CO3 at the 0.5 nm thick interfaces is 24% higher than the bulk value at 308 K and 12% higher at 333 K. This implies a lower accumulation of acid molecules at the interface as temperature increases. In the bulk phase, the density of CO2 in water at 308 K is 33 kg/m3, which corresponds to 43 molecules of CO2 in water. If the number of H2CO3 molecules is reduced by 43 (i.e., from 260 to 217), the IFT variation remains within the simulation error. This is also the case at 333 K. Therefore, it is safe to neglect the number of D
DOI: 10.1021/acs.langmuir.5b00445 Langmuir XXXX, XXX, XXX−XXX
Article
Langmuir
Figure 4. MD snapshots of CO2/water (11.7 MPa, 333 K) and CO2/water/quartz (11.7 MPa, 308 K) systems in the presence of carbonic acid. Colors in the fluid phase refer to CO2 (green), H2O (blue), and H2CO3 (red).
vary with temperature as the variations are within the simulation error. 3.3.2. With Reaction. MD simulations of water/CO2/quartz system were performed in the NVT ensemble by taking into consideration reaction 3 between CO2 and water. Quartz, on the other hand, is a relatively stable mineral that reacts slowly and to a small extent with water and CO2. Hence we did not consider any chemical reactions between quartz and the fluids in our system. A snapshot of the simulation is illustrated in Figure 4b where a water droplet surrounded by CO2 molecules is clearly seen on the quartz surface. The simulation results depicted as solid lines in Figure 6 show that the contact angle of water/CO2/quartz slightly decreased by 1−4 deg in the presence of H2CO3, across the pressure range investigated. A possible explanation of this behavior is the increased H-bonding between fluid (i.e., H2O and H2CO3) molecules and quartz surface (see Table 4), which increases the hydrophilicity of the surface. The extent of CA reduction could be loosely correlated to the number of H-bonds formed between fluid molecules and quartz surface. Taking into consideration reaction 3 in our simulations has improved the predictions at low temperature. Indeed, our results at 308 K are now comparable to equilibrium values measured in a previous study.36 They also lay between advancing and receding angles measured at the same temperature.19 At 333 K, the simulated CA decreased by 1−2 deg and became closer to measured receding angles. It is possible that the degree of protonation of the quartz surface decreases as temperature increases, which may explain the underestimated CA values. Overall, the effect of carbonic acid
CO2 molecules in water and represent both carbonic acid and solvated CO2 by H2CO3. Results of the MD simulations are presented in Figure 3 as solid lines. The addition of carbonic acid to the system had more effect on IFT at high pressures. This is due to the fact that the concentration of H2CO3 in water increased (see Table 4), due to a higher dissolution of CO2 in water. As a result, the IFT values decreased and became in better agreement with experimental data.19,35 Under the same pressure, the impact of reaction 3 on IFT is smaller at 333 K than 308 K. The reason is twofold: (i) there is a lower number of H2CO3 molecules in water at 333 K since CO2 dissolution decreases as temperature increases, and (ii) the accumulation of acid molecules at the interface decreases at 333 K due to the fact that the entropic contribution to the system’s free energy becomes stronger, which favors more homogeneous mixtures (see Figure 5). 3.3. Contact Angle of CO2/Water/Quartz System. 3.3.1. Without Reaction. The equilibrium contact angles were calculated at different P and T using ImageJ software and eq 13. Five snapshots of MD simulations were analyzed at each temperature and pressure. The average values are presented in Figure 6 with a maximum error of 1 deg. At 308 K, the predicted contact angles are similar to experimental advancing angles and follow the same trend.19 The increase in CA with pressure could be interpreted as a wettability alteration of the quartz surface to less water-wet due to an increase in CO2 density and the associated higher CO2−quartz interactions.10 At 333 K, the simulated angles are closer to experimental receding ones.19 The contact angles did not E
DOI: 10.1021/acs.langmuir.5b00445 Langmuir XXXX, XXX, XXX−XXX
Article
Langmuir
Figure 5. Density profiles of water, CO2, and carbonic acid at 11.7 MPa and different temperatures: (a) 308 K, (b) 333 K.
Figure 6. Effect of pressure on contact angle of CO2/water/quartz system at (a) 308 K and (b) 333 K. MD simulations performed with and without reaction between CO2 and water according to reaction 3.
Table 4. Effect of H2CO3 Concentration in Water on Number of H-Bonds between Fluid Molecules and Quartz and CA Decrease (Δθ) at Different P and T Conditions pressure (MPa)
[H2CO3] (mol/nm3)
3.4 6.2 9.0 11.7
1.04 1.87 2.7 3.54
3.4 6.2 9.0 11.7
0.57 1.029 1.49 1.95
avg no. H-bonds between H2CO3 and quartz 308 K 32 74 84 73 333 K 32 60 87 68
avg no. H-bonds between H2CO3 + water and quartz
Δθ (deg)
2253 2735 2498 2194
1.9 3.1 3.7 2.0
2821 2809 2906 2415
2.4 1.2 2.8 1.2
on IFT and CA appears to be more pronounced at low temperature and elevated pressure due to the higher dissolution of CO2 in water.
4. CONCLUSION The interfacial tension and contact angle of water/CO2/quartz systems were investigated using MD simulations with TIP4P2005 and TraPPE models for water and CO 2 , respectively, and a new molecular model for α-quartz surface. Two sets of simulations were performed at different temperature and pressure conditions: (i) one without any carbonic acid (i.e., no reaction between CO2 and water), which was similar to previous MD studies, and (ii) another with carbonic acid by considering chemical reactions between dissolved CO2 F
DOI: 10.1021/acs.langmuir.5b00445 Langmuir XXXX, XXX, XXX−XXX
Article
Langmuir
(13) Vega, C.; de Miguel, E. Surface Tension of the Most Popular Models of Water by Using the Test-Area Simulation Method. J. Chem. Phys. 2007, 126, 154707. (14) Abascal, J. L. F.; Vega, C. A General Purpose Model for the Condensed Phases of Water: TIP4P/2005. J. Chem. Phys. 2005, 123, 234505. (15) Liu, Y.; Lafitte, T.; Panagiotopoulos, A. Z.; Debenedetti, P. G. Simulations of Vapor−Liquid Phase Equilibrium and Interfacial Tension in the CO2−H2O−NaCl System. AIChE J. 2013, 59, 3514−3522. (16) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926−935. (17) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. In Intermolecular Forces; Pullman, B., Eds.; Reidel: Dordrecht, 1981; p 331. (18) Jorgensen, W. L.; Jenson, C. Temperature Dependence of TIP3P, SPC, and TIP4P Water from NPT Monte Carlo Simulations: Seeking Temperatures of Maximum Density. J. Comput. Chem. 1998, 19, 1179−1186. (19) Saraji, S.; Goual, L.; Piri, M.; Plancher, H. Wettability of Supercritical Carbon Dioxide/Water/Quartz Systems: Simultaneous Measurement of Contact Angle and Interfacial Tension at Reservoir Conditions. Langmuir 2013, 29, 6856−6866. (20) Drever, J. I. The geochemistry of natural waters: Surface and groundwater environments, 3rd ed.; Prentice Hall: Upper Saddle River, NJ, 1997. (21) Alejandre, J.; Tildesley, D. J.; Chapela, G. A. Molecular Dynamics Simulation of the Orthobaric Densities and Surface Tension of Water. J. Chem. Phys. 1995, 102, 4574−4583. (22) Mkhonto, D.; Leeuw, N. H. de. The Effect of Surface Silanol Groups on the Deposition of Apatite onto Silica Surfaces: A Computer Simulation Study. J. Mater. Sci: Mater. Med. 2008, 19, 203−216. (23) Zhuravlev, L. T. Concentration of Hydroxyl Groups on the Surface of Amorphous Silicas. Langmuir 1987, 3 (3), 316−318. (24) Cygan, R. T.; Liang, J.-J.; Kalinichev, A. G. Molecular Models of Hydroxide, Oxyhydroxide, and Clay Phases and the Development of a General Force Field. J. Phys. Chem. B 2004, 108, 1255−1266. (25) Skelton, A. A.; Fenter, P.; Kubicki, J. D.; Wesolowski, D. J.; Cummings, P. T. Simulations of the Quartz(101̅1)/Water Interface: A Comparison of Classical Force Fields, Ab Initio Molecular Dynamics, and X-ray Reflectivity Experiments. J. Phys. Chem. C 2011, 115, 2076− 2088. (26) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. J. Am. Chem. Soc. 1996, 118, 11225−11236. (27) Jorgensen, W. L.; Laird, E. R.; Nguyen, T. B.; Tirado-Rive, J. Monte Carlo Simulations of Pure Liquid Substituted Benzenes with OPLS Potential Functions. J. Comput. Chem. 1993, 14, 206−215. (28) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, Ö .; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09; Gaussian, Inc.: Wallingford, CT, 2009. (29) Berendsen, H. J. C.; Postma, J. P. M.; Gunsteren, W. F. v.; DiNola, A.; Haak, J. R. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81, 3684−3690.
and water. While the presence of carbonic acid in water had a minor effect on quartz wettability, it improved the IFT predictions at low temperature and elevated pressure due to the higher dissolution of CO2 in water. The good agreement between modeled and measured properties also validates the models adopted in this work. Future work will be directed toward understanding the impact of temperature on the extent of surface protonation. The proposed methodology could then be used as a platform to accurately predict the impact of brine chemistry and CO2 cocontaminants on the interfacial properties of CO2/water/quartz systems in order to determine the risk of gas leakage though the cap rock under a wide range of pressures and temperatures.
■
AUTHOR INFORMATION
Corresponding Author
*E-mail: [email protected]. Phone: 307-766-3278. Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS The authors thank the School of Energy Resources at the University of Wyoming for financial support and Dr. Mohammad Piri for providing access to his computer cluster. We greatly appreciate the insightful discussions with Dr. Carrick Eggleston on geochemical reactions.
■
REFERENCES
(1) Chiquet, P.; Broseta, D.; Thibeau, S. Wettability Alteration of Caprock Minerals by Carbon Dioxide. Geofluids 2007, 7, 112−122. (2) Iglauer, S.; Pentland, C. H.; Busch, A. CO2 Wettability of Seal and Reservoir Rocks and the Implications for Carbon GeoSequestration. Water Resour. Res. 2015, 51 (1), 729−774. (3) Purcell, W. R. Interpretation of Capillary Pressure Data. J. Pet. Technol. 1950, 2 (08), 11−12. (4) Harkins, W. D.; Feldman, A. Films. The Spreading of Liquids and the Spreading Coefficient. J. Am. Chem. Soc. 1922, 44, 2665−2685. (5) Zhang, Z.; Duan, Z. An Optimized Molecular Potential for Carbon Dioxide. J. Chem. Phys. 2005, 122, 214507. (6) Liu, Y.; Panagiotopoulos, A. Z.; Debenedetti, P. G. Monte Carlo Simulations of High-Pressure Phase Equilibria of CO2−H2O Mixtures. J. Phys. Chem. B 2011, 115, 6629−6635. (7) Kvamme, B.; Kuznetsova, T.; Hebach, A.; Oberhof, A.; Lunde, E. Measurements and Modelling of Interfacial Tension for Water + Carbon Dioxide Systems at Elevated Pressures. Comput. Mater. Sci. 2007, 38, 506−513. (8) Li, X.; Ross, D. A.; Trusler, J. P. M.; Maitland, G. C.; Boek, E. S. Molecular Dynamics Simulations of CO2 and Brine Interfacial Tension at High Temperatures and Pressures. J. Phys. Chem. B 2013, 117, 5647−5652. (9) Liu, S.; Yang, X.; Qin, Y. Molecular Dynamics Simulation of Wetting Behavior at CO2/Water/Solid Interfaces. Chin. Sci. Bull. 2010, 55, 2252−2257. (10) Iglauer, S.; Mathew, M. S.; Bresme, F. Molecular Dynamics Computations of Brine−CO2 Interfacial Tensions and Brine−CO2− Quartz Contact Angles and Their Effects on Structural and Residual Trapping Mechanisms in Carbon Geo-Sequestration. J. Colloid Interface Sci. 2012, 386, 405−414. (11) McCaughan, J.; Iglauer, S.; Bresme, F. Molecular Dynamics Simulation of Water/CO2−Quartz Interfacial Properties: Application to Subsurface Gas Injection. Energy Procedia 2013, 37, 5387−5402. (12) Nielsen, L. C.; Bourg, I. C.; Sposito, G. Predicting CO2−Water Interfacial Tension under Pressure and Temperature Conditions of Geologic CO2 Storage. Geochim. Cosmochim. Acta 2012, 81, 28−38. G
DOI: 10.1021/acs.langmuir.5b00445 Langmuir XXXX, XXX, XXX−XXX
Article
Langmuir (30) Parrinello, M.; Rahman, A. Strain Fluctuations and Elastic Constants. J. Chem. Phys. 1982, 76, 2662−2666. (31) Nosé, S. A Molecular Dynamics Method for Simulations in the Canonical Ensemble. Mol. Phys. 1984, 52, 255−268. (32) Schneider, C. A.; Rasband, W. S.; Eliceiri, K. W. NIH Image to ImageJ: 25 Years of Image Analysis. Nat. Methods 2012, 9, 671−675. (33) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435−447. (34) Humphrey, W.; Dalke, A.; Schulten, K. VMDVisual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33−38. (35) Sarmadivaleh, M.; Al-Yaseri, A. Z.; Iglauer, S. Influence of Temperature and Pressure on Quartz−Water−CO2 Contact Angle and CO2−Water Interfacial Tension. J. Colloid Interface Sci. 2015, 441, 59−64. (36) Farokhpoora, R.; Bjørkvik, B.; Lindebergb, E.; Torsætera, O. Wettability Behaviour of CO2 at Storage Conditions. Int. J. Greenhouse Gas Control 2013, 12, 18−25. (37) Iglauer, S.; Salamah, A.; Sarmadivaleh, M.; Liu, K.; Phan, C. Contamination of Silica Surfaces: Impact on Water−CO2−Quartz and Glass Contact Angle Measurements. Int. J. Greenhouse Gas Control 2014, 22, 325−328. (38) Vargaftik, N. B.; Volkov, B. N.; Voljak, L. D. International Tables of the Surface Tension of Water. J. Phys. Chem. Ref. Data 1983, 12, 817−820.
H
DOI: 10.1021/acs.langmuir.5b00445 Langmuir XXXX, XXX, XXX−XXX
Molecular Dynamics Simulations of CO2/Water/Quartz Interfacial Properties: Impact of CO2 Dissolution in Water Gina Javanbakht, Mohammad Sedghi, William Welch, and Lamia Goual* Department of Chemical and Petroleum Engineering, University of Wyoming, 1000 E. University Avenue, Laramie, Wyoming 82071, United States ABSTRACT: The safe trapping of carbon dioxide (CO2) in deep saline aquifers is one of the major concerns of CO2 sequestration. The amount of capillary trapping is dominated by the capillary pressure of water and CO2 inside the reservoir, which in turn is controlled by the interfacial tension (IFT) and the contact angle (CA) of CO2/water/rock systems. The measurement of IFT and CA could be very challenging at reservoir conditions, especially in the presence of toxic cocontaminants. Thus, the ability to accurately predict these interfacial properties at reservoir conditions is very advantageous. Although the majority of existing molecular dynamics (MD) studies of CO2/water/mineral systems were able to capture the trends in IFT and CA variations with pressure and temperature, their predictions often deviated from experimental data, possibly due to erroneous models and/or overlooked chemical reactions. The objective of this study was to improve the MD predictions of IFT and CA of CO2/ water/quartz systems at various pressure and temperature conditions by (i) considering the chemical reactions between CO2 and water and (ii) using a new molecular model for α-quartz surface. The results showed that the presence of carbonic acid at the CO2/water interface improved the predictions of IFT, especially at low temperature and high pressure where more CO2 dissolution occurs. On the other hand, the effect on CA was minor. The slight decrease in CA observed across the pressure range investigated could be attributed to an increase in the total number of H-bonds between fluid molecules and quartz surface. that among molecular models for water, TIP4P/200514 could produce the closest match to experimental data. Other researchers have reached the same conclusion by using a dummy particle to carry the negative charge of oxygen atom in various water models.14,15 Older water models (e.g., TIP3P,16 SPC,17 TIP4P,14 and SPC/E18) were not as accurate as TIP4P/ 2005 in predicting the behavior of water molecules during MD simulations. Several models were also suggested for CO2 such as MSM, EPM2, TraPPE, Errington, etc.5 In order to find the best model for CO2/water systems at high pressure (P) and temperature (T), Liu et al. examined several existing models for water (SPC, TIP4P, TIP4P2005) and CO2 (EPM2, TraPPE) and found that TraPPE/TIP4P2005 and EPM2/SPC combinations provided the most accurate predictions of CO2 solubility and IFT variations with P and T.6,16 Kvamme et al. used the SPC/E model for water and a simple EMP model for CO2 to study the IFT between CO2 and water at 278−335 K temperatures and 0.1−20 MPa pressures. Their results were in good agreement with experimental data and revealed a decrease in IFT with increasing pressure.7 More recently, Li et al. studied the IFT of CO2/brine systems at higher P and T using TIP4P and SPC/E models for water and EPM2 model for CO2.8 They showed that the SPC/E water model gives a better match between experimental and simulation data. However, the simulation errors were higher at elevated pressures. In addition to IFT, the CA between two fluid phases and a solid phase controls the capillary pressure of the system. Using
1. INTRODUCTION The successful sequestration of carbon dioxide (CO2) in deep geological formations requires a minimum leakage of CO2 through the cap rock. The amount of leakage could be estimated from the threshold capillary pressure, which in turn is affected by the interfacial tension (IFT) between CO2 and the present fluids, as well as the contact angle (CA) between fluids and the solid surface, according to the following equation1 Pc = [2γα , β cos(θα , β , δ)]/r
(1)
where γα,β is the interfacial tension between phases α and β and θα,β,δ is the contact angle between phases α and β and solid surface δ. Although this relationship applies to perfect capillary tubes, it is more complicated in pore spaces as pore geometry affects the capillary pressure.2,3 In addition to CO2 sequestration, the prediction of IFT in CO2 enhanced oil recovery (EOR) processes is important since IFT (or γ) controls the spreading coefficient of one fluid over the others4 according to S = γgas,brine − γoil,gas − γoil,brine
(2)
The measurement of IFT and CA at reservoir conditions is usually very challenging, especially in the presence of cocontaminants such as SO2 and NOx. Therefore, the ability to accurately predict these properties at reservoir conditions can be very useful. Several research groups have used molecular dynamics (MD) modeling to estimate the IFT and CA of CO2/ water/quartz systems under various conditions.5−12 In 2007, Vega and Miguel calculated the surface tension of water in order to compare different water models.13 Their study showed © XXXX American Chemical Society
Received: February 3, 2015 Revised: April 14, 2015
A
DOI: 10.1021/acs.langmuir.5b00445 Langmuir XXXX, XXX, XXX−XXX
Article
Langmuir the same models, Liu et al. investigated the effect of pressure on the wetting behavior of CO2/water/quartz systems by considering both β-quartz and quartz with Silanol groups.9 They found that the CA increased as pressure increased for both surfaces. However, unlike hydrophilic quartz, water droplets on hydrophobic quartz detached from the surface by increasing CO2 pressure. More recently, Iglauer et al. studied the IFT and CA of CO2/water/α-quartz systems using the EMP2/TIP4P2005 model combination up to 20 MPa pressure.10 The partial charges of quartz were 2.4 for Si and −1.2 for oxygen atoms. Their results showed an increase in CA with increasing pressure. IFT, on the other hand, decreased with increasing pressure and increased with increasing temperature. The model used in that study was accurate in predicting IFT versus pressure below the saturation pressure CO2 at 300 K (i.e., subcritical CO2) but overpredicted IFT by 20% at higher pressures (i.e., supercritical CO2). This study was continued by McCaughan et al., who examined the effect of divalent salts (such as CaCl2) on the wettability of CO2/H2S/ water/α-quartz systems during subsurface gas injection.11 They reported that salt addition did not affect the CA. On the other hand, H2S interacted more strongly with quartz than CO2, resulting in a higher CA under the same P and T conditions. Although previous MD studies could capture the trends in IFT and CA variations with P and T, the simulation results often deviated from experimental data, especially at high P and T conditions. A plausible explanation for this discrepancy is that, by default, MD simulations do not explicitly take into account the chemical reactions that may occur at fluid/fluid and fluid/solid interfaces. Moreover, the molecular models for the solid surface used in previous CA predictions might have been inaccurate. The goal of this study was to improve the predictions of IFT and CA of CO2/water/α-quartz systems under a wide range of T and P (from subcritical to supercritical CO2) by considering chemical reactions between CO2 and water and using a novel molecular model for the α-quartz surface.
HCO−3 ⇌ H+ + CO32 −
(6)
When combined with eq 4 the expression of K1 becomes a H+a HCO−3 K1 = K CO2PCO2 Since aH = a +
− HCO3
(8) +
due to charge balance, the activity of H reduces
to
a H+ =
K1K CO2PCO2
(9)
By increasing CO2 pressure from 3.4 to 12 MPa, the pH (i.e., −log aH+) decreases from 3.15 to 2.80 at 308 K and from 3.26 to 3.00 at 333 K, due to the presence of carbonic acid. Under these low pH conditions the dissociation of H2CO3 is negligible. For example at a pH of 2.88, the concentration of H3O+ is 10−2.88 M and is incredibly small in comparison to that of water (55 M). The stoichiometric ratio of hydronium/water is 2.4 × 10−5:1. In our simulations, which contain 4178 water molecules, this would be a small fraction of a hydronium ion. However, there is a considerable amount of solvated CO2/ carbonic acid at the interface, which can have a sizable impact on the IFT between CO2 and water. In our simulations, we used H2CO3 to simulate the solvation of CO2, as suggested by ref 20. At each P and T condition, the concentration of carbonic acid was calculated, and the corresponding number of H2CO3 molecules was added to the simulation box. 2.2. Intermolecular Forces and MD Simulations. The intermolecular potential forces used in this work are divided into van der Waals interactions and electrostatic interactions. The van der Waals interactions were calculated from Lennard-Jones potential
⎡⎛ σij ⎞12 ⎛ σij ⎞6 ⎤ Uij = 4εij⎢⎜ ⎟ − ⎜ ⎟ ⎥ ⎝r ⎠⎦ ⎣⎝ r ⎠
2.1. Chemical Reactions. Reactions involving carbonate systems control the pH of most natural waters. When carbon dioxide comes in contact with water, it starts to dissolve in the water until an equilibrium state is reached. The concentration of dissolved CO2 depends on the P and T conditions. The dissolution of CO2 in water and subsequent reaction is presented by the following20
(10)
where εij is the interaction strength between atoms i and j at a distance r, σij is the effective diameter of atoms i and j where σij = ((σii + σjj))/2 and εij = ((εiiεij)1/2), σii is the diameter of atom i, and εij is the interaction strength between atoms i and j. Electrostatic interactions were calculated from Coulomb equation
(3)
The right-hand side of the reaction above includes both carbonic acid and solvated CO2 in water. In fact the amount of solvated CO2 is much greater than carbonic acid. In other words, both carbonic acid and solvated CO2 are represented by H2CO3.20 The activity of H2CO3 can be obtained from the equilibrium constant (KCO2) and the pressure of carbon dioxide in the CO2-rich phase (PCO2) a H2CO3 = K CO2PCO2
(5)
The equilibrium constants K1 and K2 of reactions 5 and 6 are tabulated in ref 20 at various temperatures. At 308 K, K1 = 10−6.33 and K2 = 10−10.29 whereas at 333 K, K1 = 10−6.29 and K2 = 10−10.14. The constant K2 is very small compared to KCO2 and K1 and can therefore be neglected. Based on reaction 5, K1 can be written as a H+a HCO−3 K1 = a H2CO3 (7)
2. METHODS
CO2 (g) + H 2O ⇌ H 2CO3(aq)
H 2CO3 ⇌ H+ + HCO−3
Fij⃗ =
Ke*Q i*Q jriĵ rij2
(11)
where Ke is Coulomb’s constant, Qα is the absolute value of the partial charge of particle α, and rij is the separation distance between particles i and j. The long-range interactions of electrostatic forces were computed using particle-mesh Ewald summation (PME). The cutoff radios for LJ interactions were set to 1.5 nm, and the periodic boundary condition was applied in all directions. The interfacial tension of the CO2/water system was calculated based on the mechanical route from the diagonal components of the pressure tensor normal to the surface according to
(4)
where KCO2 is equal to 10−1.52 at 308 K and 10−1.78 at 333 K20 and PCO2 varies based on the applied pressure. H2CO3 is neutral and its activity coefficient is equal to unity, meaning that the activity of the carbonic acid in water is equal. Carbonic acid in aqueous media starts to dissociate into bicarbonate and hydrogen as the pH increases. Bicarbonate itself can further dissociate to carbonate and proton at high pH conditions according to20
γ=
⎤ 1 ⎡ 1 LZ ⎢PZZ − (PXX + PYY )⎥ ⎦ 2 ⎣ 2
(12)
where γ is the surface tension, LZ is the system length along the z-axis, and Pxx, Pyy, and PZZ are the components of the pressure tensor.21 B
DOI: 10.1021/acs.langmuir.5b00445 Langmuir XXXX, XXX, XXX−XXX
Article
Langmuir
Figure 1. α-Quartz surface: (a) top view and (b) side view. The solid surface consists of a fully coordinated α-quartz with Si− O−Si bridges on the surface, as seen in Figure 1. The surface density of Silanol groups is 4.9 nm−2, which is in line with experimental data reported for various silica surfaces.22,23 Half of the nonbridge oxygen atoms at the surface were protonated. Half-protonated quartz gave the best contact angle compared to fully protonated and nonprotonated surfaces. The unsaturated oxygen atoms as well as the Silanol groups on the surface were considered flexible with the bonded interactions presented in Table 1. The rest of atoms in the quartz were kept frozen
simulation for 500 ps with Berendsen thermostat and barostat.29 The production simulations were performed using the Nose−Hoover temperature bath and the Parrinello−Rahman pressure coupling30 for 10 ns. For the simulations containing carbonic acid, the concentration of H2CO3 was calculated at each pressure and temperature according to the procedure described in Section 2.1, and a corresponding number of carbonic acid molecules was added to the simulation box. The number of H2CO3 molecules varied between 70 and 260. For contact angle calculations, the simulation box had dimensions of 29 × 29 × 80 nm3 and contained 21 000 SiO4 units in the solid surface. We placed 3548 water molecules on the solid surface and used 29 000 and 58 000 CO2 molecules for low- and high-pressure simulations, respectively. All simulations were performed in NVT ensembles using the Nose−Hoover thermostat.31 The pressure of the CO2 phase was controlled by exerting an external pressure on a piston composed of virtual hard sphere particles in a face-centered cubic (fcc) arrangement. The hard sphere wall was modeled using repulsive interactions (i.e., LJ particles with low attraction potential). To verify the pressure of each simulation, the average CO2 density was calculated and compared with the experimental values, as shown in Table 2. The densities obtained from MD simulations closely match
Table 1. LJ Parameters and Partial Charges of Components Considered in MD Simulations Intermolecular Interactions atoms
ε (kJ/mol)
Si (quartz) 7.700 × 10−6 O (Si−O−Si) 0.6500 O (Si−OH) 0.6500 O (Si−O) 0.6500 H (H2CO3) 0.00 C (H2CO3) 0.43932 O (O−H in H2CO3) 0.71128 O (H2CO3) 0.87864 Intramolecular Interactions in bonds
Kb (kJ/Å2)
σ (Å) 3.3020 3.1655 3.1655 3.1655 0.00 3.7500 3.000 2.960 Silica
q (e) 2.10 −1.050 −0.950 −0.525 0.429 0.88 −0.57 −0.602
Table 2. Measured and Simulated Density of CO2 at Different P and T Conditions density (kg/m3)
r0 (Å)
Si−O O−H angle
2510.40 1627.50 Kθ (kJ/rad2)
1.630 0.950 θ0 (deg)
Si−O−Si
397.488
109.5
pressure (MPa) 3.4 6.2 9.0 11.7
during the simulations. Since the frozen atoms are separately thermostated and are not bonded to the surface atoms, we expect that freezing and restraining their positions will have a negligible impact on the simulation results. ClayFF24 parameters, as reported in ref 25, were used to represent nonbonded interactions for quartz (see Table 1). LJ parameters for carbonic acid atoms were chosen from optimized potentials for liquid simulations-all atoms (OPLS-aa) force field.26,27 To determine the partial charges of atoms in carbonic acid molecules, quantum mechanical (QM) calculations were performed to optimize the structure of H2CO3 molecule using the Gaussian software.28 For QM computations, MP2 level of theory combined with 6-31G(d) basis set was employed, and the charges were calculated by ESP method. The LJ potentials and partial charges of the atoms in the particles used in MD simulations are provided in Table 1. The size of the simulation box for IFT calculation was 5 × 5 × 10 nm3. The box contained 4178 molecules of water and 700 molecules of CO2. Initially, equilibrium conditions were reached by running NPT
3.4 6.2 9.0 11.7
simulation 308 K 79.34 175.24 645.8 761.9 333 K 69.39 135.8 230.02 378.92
exptl19 72.3 175.2 610.4 739.2 63 131.2 233.8 398
the experimental values. At one temperature and pressure condition, a simulation was run for 5 ns, and the contact angle was measured every 500 ps. The contact angle values did not change significantly after 4 ns, and hence all other simulations were performed for 5 ns to ensure the system has reached equilibrium. The contact angles were calculated by the following equation using ImageJ software, where a is the droplet height and b is the droplet base radius,32 as shown in Figure 2.
θ = arcsin C
2ab a 2 + b2
(13) DOI: 10.1021/acs.langmuir.5b00445 Langmuir XXXX, XXX, XXX−XXX
Article
Langmuir
Figure 2. Parameters used in eq 13. Periodic boundary conditions were applied in the X and Y directions. All the simulation were performed using GROMACS 4.6.333 and VMD software34 was used to visualize the simulations results.
3. RESULTS AND DISCUSSION We first ran MD simulations to calculate the surface tension of water and select the best models for water and CO2. We then applied those models to calculate the IFT between CO2 and water and the CA of CO2/water/quartz system (with and without reaction 3). The simulated IFT and CA were compared to experimental data reported by Saraji et al.,19 Sarmadivaleh et al.,35 and Farokhpoor et al.36 These reference data were selected due to their high precision and accuracy.37 3.1. Surface Tension of Water. We tested three different water models (SPC/E, TIP4P, and TIP4P2005) to calculate the surface tension of water. We placed 4050 water molecules in a 4 × 4 × 8 nm3 simulation box and conducted NVT simulations for 1 ns at a temperature of 300 K. Among all models, TIP4P/ 2005 presented the closest surface tensions to experimental data38 (see Table 3). Therefore, we continued our study in the Table 3. Comparison between Measured and Simulated Surface Tension of Water at 300 K Using Different Water Models experiment38
water model surface tension (mN/m)
TIP4P2005
TIP4P
SPC/E
68.3
50
63.8
71.8
next sections using TIP4P/2005 as the water model. We also selected the TraPPE model for CO2 based on its accuracy when combined with TIP4P/2005 water model.6 3.2. Interfacial Tension between CO2 and Water. 3.2.1. Without Reaction. MD simulations were first performed to estimate the IFT between CO2 and water at different pressures (3.4−11.7 MPa) and temperatures (308, 333 K), without considering any reactions between CO2 molecules and water (see dashed lines in Figure 3). The results revealed that IFT decreases with increasing pressure, in good agreement with experimental data measured by axisymmetric drop shape analysis (ADSA) and other methods.19,35 However, like previous MD studies,7,12 the IFT between supercritical CO2 and water was overestimated. The estimated error in the simulations was between 0.25 to 1.1 mN/m. The variation of IFT values with temperature is within the error estimate and is negligible. 3.2.2. With Reaction. Using the same procedure described earlier, NPT simulations of CO2/water IFT were performed by taking into consideration the dissolution of CO2 in water according to reaction 3. The concentration of H2CO3 in water was calculated at each P and T, and the corresponding number of molecules were added to the simulation box (see Section 2.1). A snapshot of CO2, water, and H2CO3 molecules is presented in Figure 4 a where H2CO3 molecules are depicted as red particles in the aqueous phase. This figure shows that some
Figure 3. Effect of pressure on IFT of the CO2/water system at (a) 308 K and (b) 333 K. MD simulations performed with and without reaction between CO2 and water according to reaction 3.
acid molecules accumulate at the CO2/water interface, affecting the IFT of the system. The density profiles of water, CO2, and carbonic acid are shown in Figure 5 at both temperatures and 11.7 MPa. In this figure, the density of H2CO3 at the 0.5 nm thick interfaces is 24% higher than the bulk value at 308 K and 12% higher at 333 K. This implies a lower accumulation of acid molecules at the interface as temperature increases. In the bulk phase, the density of CO2 in water at 308 K is 33 kg/m3, which corresponds to 43 molecules of CO2 in water. If the number of H2CO3 molecules is reduced by 43 (i.e., from 260 to 217), the IFT variation remains within the simulation error. This is also the case at 333 K. Therefore, it is safe to neglect the number of D
DOI: 10.1021/acs.langmuir.5b00445 Langmuir XXXX, XXX, XXX−XXX
Article
Langmuir
Figure 4. MD snapshots of CO2/water (11.7 MPa, 333 K) and CO2/water/quartz (11.7 MPa, 308 K) systems in the presence of carbonic acid. Colors in the fluid phase refer to CO2 (green), H2O (blue), and H2CO3 (red).
vary with temperature as the variations are within the simulation error. 3.3.2. With Reaction. MD simulations of water/CO2/quartz system were performed in the NVT ensemble by taking into consideration reaction 3 between CO2 and water. Quartz, on the other hand, is a relatively stable mineral that reacts slowly and to a small extent with water and CO2. Hence we did not consider any chemical reactions between quartz and the fluids in our system. A snapshot of the simulation is illustrated in Figure 4b where a water droplet surrounded by CO2 molecules is clearly seen on the quartz surface. The simulation results depicted as solid lines in Figure 6 show that the contact angle of water/CO2/quartz slightly decreased by 1−4 deg in the presence of H2CO3, across the pressure range investigated. A possible explanation of this behavior is the increased H-bonding between fluid (i.e., H2O and H2CO3) molecules and quartz surface (see Table 4), which increases the hydrophilicity of the surface. The extent of CA reduction could be loosely correlated to the number of H-bonds formed between fluid molecules and quartz surface. Taking into consideration reaction 3 in our simulations has improved the predictions at low temperature. Indeed, our results at 308 K are now comparable to equilibrium values measured in a previous study.36 They also lay between advancing and receding angles measured at the same temperature.19 At 333 K, the simulated CA decreased by 1−2 deg and became closer to measured receding angles. It is possible that the degree of protonation of the quartz surface decreases as temperature increases, which may explain the underestimated CA values. Overall, the effect of carbonic acid
CO2 molecules in water and represent both carbonic acid and solvated CO2 by H2CO3. Results of the MD simulations are presented in Figure 3 as solid lines. The addition of carbonic acid to the system had more effect on IFT at high pressures. This is due to the fact that the concentration of H2CO3 in water increased (see Table 4), due to a higher dissolution of CO2 in water. As a result, the IFT values decreased and became in better agreement with experimental data.19,35 Under the same pressure, the impact of reaction 3 on IFT is smaller at 333 K than 308 K. The reason is twofold: (i) there is a lower number of H2CO3 molecules in water at 333 K since CO2 dissolution decreases as temperature increases, and (ii) the accumulation of acid molecules at the interface decreases at 333 K due to the fact that the entropic contribution to the system’s free energy becomes stronger, which favors more homogeneous mixtures (see Figure 5). 3.3. Contact Angle of CO2/Water/Quartz System. 3.3.1. Without Reaction. The equilibrium contact angles were calculated at different P and T using ImageJ software and eq 13. Five snapshots of MD simulations were analyzed at each temperature and pressure. The average values are presented in Figure 6 with a maximum error of 1 deg. At 308 K, the predicted contact angles are similar to experimental advancing angles and follow the same trend.19 The increase in CA with pressure could be interpreted as a wettability alteration of the quartz surface to less water-wet due to an increase in CO2 density and the associated higher CO2−quartz interactions.10 At 333 K, the simulated angles are closer to experimental receding ones.19 The contact angles did not E
DOI: 10.1021/acs.langmuir.5b00445 Langmuir XXXX, XXX, XXX−XXX
Article
Langmuir
Figure 5. Density profiles of water, CO2, and carbonic acid at 11.7 MPa and different temperatures: (a) 308 K, (b) 333 K.
Figure 6. Effect of pressure on contact angle of CO2/water/quartz system at (a) 308 K and (b) 333 K. MD simulations performed with and without reaction between CO2 and water according to reaction 3.
Table 4. Effect of H2CO3 Concentration in Water on Number of H-Bonds between Fluid Molecules and Quartz and CA Decrease (Δθ) at Different P and T Conditions pressure (MPa)
[H2CO3] (mol/nm3)
3.4 6.2 9.0 11.7
1.04 1.87 2.7 3.54
3.4 6.2 9.0 11.7
0.57 1.029 1.49 1.95
avg no. H-bonds between H2CO3 and quartz 308 K 32 74 84 73 333 K 32 60 87 68
avg no. H-bonds between H2CO3 + water and quartz
Δθ (deg)
2253 2735 2498 2194
1.9 3.1 3.7 2.0
2821 2809 2906 2415
2.4 1.2 2.8 1.2
on IFT and CA appears to be more pronounced at low temperature and elevated pressure due to the higher dissolution of CO2 in water.
4. CONCLUSION The interfacial tension and contact angle of water/CO2/quartz systems were investigated using MD simulations with TIP4P2005 and TraPPE models for water and CO 2 , respectively, and a new molecular model for α-quartz surface. Two sets of simulations were performed at different temperature and pressure conditions: (i) one without any carbonic acid (i.e., no reaction between CO2 and water), which was similar to previous MD studies, and (ii) another with carbonic acid by considering chemical reactions between dissolved CO2 F
DOI: 10.1021/acs.langmuir.5b00445 Langmuir XXXX, XXX, XXX−XXX
Article
Langmuir
(13) Vega, C.; de Miguel, E. Surface Tension of the Most Popular Models of Water by Using the Test-Area Simulation Method. J. Chem. Phys. 2007, 126, 154707. (14) Abascal, J. L. F.; Vega, C. A General Purpose Model for the Condensed Phases of Water: TIP4P/2005. J. Chem. Phys. 2005, 123, 234505. (15) Liu, Y.; Lafitte, T.; Panagiotopoulos, A. Z.; Debenedetti, P. G. Simulations of Vapor−Liquid Phase Equilibrium and Interfacial Tension in the CO2−H2O−NaCl System. AIChE J. 2013, 59, 3514−3522. (16) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926−935. (17) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. In Intermolecular Forces; Pullman, B., Eds.; Reidel: Dordrecht, 1981; p 331. (18) Jorgensen, W. L.; Jenson, C. Temperature Dependence of TIP3P, SPC, and TIP4P Water from NPT Monte Carlo Simulations: Seeking Temperatures of Maximum Density. J. Comput. Chem. 1998, 19, 1179−1186. (19) Saraji, S.; Goual, L.; Piri, M.; Plancher, H. Wettability of Supercritical Carbon Dioxide/Water/Quartz Systems: Simultaneous Measurement of Contact Angle and Interfacial Tension at Reservoir Conditions. Langmuir 2013, 29, 6856−6866. (20) Drever, J. I. The geochemistry of natural waters: Surface and groundwater environments, 3rd ed.; Prentice Hall: Upper Saddle River, NJ, 1997. (21) Alejandre, J.; Tildesley, D. J.; Chapela, G. A. Molecular Dynamics Simulation of the Orthobaric Densities and Surface Tension of Water. J. Chem. Phys. 1995, 102, 4574−4583. (22) Mkhonto, D.; Leeuw, N. H. de. The Effect of Surface Silanol Groups on the Deposition of Apatite onto Silica Surfaces: A Computer Simulation Study. J. Mater. Sci: Mater. Med. 2008, 19, 203−216. (23) Zhuravlev, L. T. Concentration of Hydroxyl Groups on the Surface of Amorphous Silicas. Langmuir 1987, 3 (3), 316−318. (24) Cygan, R. T.; Liang, J.-J.; Kalinichev, A. G. Molecular Models of Hydroxide, Oxyhydroxide, and Clay Phases and the Development of a General Force Field. J. Phys. Chem. B 2004, 108, 1255−1266. (25) Skelton, A. A.; Fenter, P.; Kubicki, J. D.; Wesolowski, D. J.; Cummings, P. T. Simulations of the Quartz(101̅1)/Water Interface: A Comparison of Classical Force Fields, Ab Initio Molecular Dynamics, and X-ray Reflectivity Experiments. J. Phys. Chem. C 2011, 115, 2076− 2088. (26) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. J. Am. Chem. Soc. 1996, 118, 11225−11236. (27) Jorgensen, W. L.; Laird, E. R.; Nguyen, T. B.; Tirado-Rive, J. Monte Carlo Simulations of Pure Liquid Substituted Benzenes with OPLS Potential Functions. J. Comput. Chem. 1993, 14, 206−215. (28) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, Ö .; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09; Gaussian, Inc.: Wallingford, CT, 2009. (29) Berendsen, H. J. C.; Postma, J. P. M.; Gunsteren, W. F. v.; DiNola, A.; Haak, J. R. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81, 3684−3690.
and water. While the presence of carbonic acid in water had a minor effect on quartz wettability, it improved the IFT predictions at low temperature and elevated pressure due to the higher dissolution of CO2 in water. The good agreement between modeled and measured properties also validates the models adopted in this work. Future work will be directed toward understanding the impact of temperature on the extent of surface protonation. The proposed methodology could then be used as a platform to accurately predict the impact of brine chemistry and CO2 cocontaminants on the interfacial properties of CO2/water/quartz systems in order to determine the risk of gas leakage though the cap rock under a wide range of pressures and temperatures.
■
AUTHOR INFORMATION
Corresponding Author
*E-mail: [email protected]. Phone: 307-766-3278. Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS The authors thank the School of Energy Resources at the University of Wyoming for financial support and Dr. Mohammad Piri for providing access to his computer cluster. We greatly appreciate the insightful discussions with Dr. Carrick Eggleston on geochemical reactions.
■
REFERENCES
(1) Chiquet, P.; Broseta, D.; Thibeau, S. Wettability Alteration of Caprock Minerals by Carbon Dioxide. Geofluids 2007, 7, 112−122. (2) Iglauer, S.; Pentland, C. H.; Busch, A. CO2 Wettability of Seal and Reservoir Rocks and the Implications for Carbon GeoSequestration. Water Resour. Res. 2015, 51 (1), 729−774. (3) Purcell, W. R. Interpretation of Capillary Pressure Data. J. Pet. Technol. 1950, 2 (08), 11−12. (4) Harkins, W. D.; Feldman, A. Films. The Spreading of Liquids and the Spreading Coefficient. J. Am. Chem. Soc. 1922, 44, 2665−2685. (5) Zhang, Z.; Duan, Z. An Optimized Molecular Potential for Carbon Dioxide. J. Chem. Phys. 2005, 122, 214507. (6) Liu, Y.; Panagiotopoulos, A. Z.; Debenedetti, P. G. Monte Carlo Simulations of High-Pressure Phase Equilibria of CO2−H2O Mixtures. J. Phys. Chem. B 2011, 115, 6629−6635. (7) Kvamme, B.; Kuznetsova, T.; Hebach, A.; Oberhof, A.; Lunde, E. Measurements and Modelling of Interfacial Tension for Water + Carbon Dioxide Systems at Elevated Pressures. Comput. Mater. Sci. 2007, 38, 506−513. (8) Li, X.; Ross, D. A.; Trusler, J. P. M.; Maitland, G. C.; Boek, E. S. Molecular Dynamics Simulations of CO2 and Brine Interfacial Tension at High Temperatures and Pressures. J. Phys. Chem. B 2013, 117, 5647−5652. (9) Liu, S.; Yang, X.; Qin, Y. Molecular Dynamics Simulation of Wetting Behavior at CO2/Water/Solid Interfaces. Chin. Sci. Bull. 2010, 55, 2252−2257. (10) Iglauer, S.; Mathew, M. S.; Bresme, F. Molecular Dynamics Computations of Brine−CO2 Interfacial Tensions and Brine−CO2− Quartz Contact Angles and Their Effects on Structural and Residual Trapping Mechanisms in Carbon Geo-Sequestration. J. Colloid Interface Sci. 2012, 386, 405−414. (11) McCaughan, J.; Iglauer, S.; Bresme, F. Molecular Dynamics Simulation of Water/CO2−Quartz Interfacial Properties: Application to Subsurface Gas Injection. Energy Procedia 2013, 37, 5387−5402. (12) Nielsen, L. C.; Bourg, I. C.; Sposito, G. Predicting CO2−Water Interfacial Tension under Pressure and Temperature Conditions of Geologic CO2 Storage. Geochim. Cosmochim. Acta 2012, 81, 28−38. G
DOI: 10.1021/acs.langmuir.5b00445 Langmuir XXXX, XXX, XXX−XXX
Article
Langmuir (30) Parrinello, M.; Rahman, A. Strain Fluctuations and Elastic Constants. J. Chem. Phys. 1982, 76, 2662−2666. (31) Nosé, S. A Molecular Dynamics Method for Simulations in the Canonical Ensemble. Mol. Phys. 1984, 52, 255−268. (32) Schneider, C. A.; Rasband, W. S.; Eliceiri, K. W. NIH Image to ImageJ: 25 Years of Image Analysis. Nat. Methods 2012, 9, 671−675. (33) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435−447. (34) Humphrey, W.; Dalke, A.; Schulten, K. VMDVisual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33−38. (35) Sarmadivaleh, M.; Al-Yaseri, A. Z.; Iglauer, S. Influence of Temperature and Pressure on Quartz−Water−CO2 Contact Angle and CO2−Water Interfacial Tension. J. Colloid Interface Sci. 2015, 441, 59−64. (36) Farokhpoora, R.; Bjørkvik, B.; Lindebergb, E.; Torsætera, O. Wettability Behaviour of CO2 at Storage Conditions. Int. J. Greenhouse Gas Control 2013, 12, 18−25. (37) Iglauer, S.; Salamah, A.; Sarmadivaleh, M.; Liu, K.; Phan, C. Contamination of Silica Surfaces: Impact on Water−CO2−Quartz and Glass Contact Angle Measurements. Int. J. Greenhouse Gas Control 2014, 22, 325−328. (38) Vargaftik, N. B.; Volkov, B. N.; Voljak, L. D. International Tables of the Surface Tension of Water. J. Phys. Chem. Ref. Data 1983, 12, 817−820.
H
DOI: 10.1021/acs.langmuir.5b00445 Langmuir XXXX, XXX, XXX−XXX
