Aqueous CO2 Solutions at Silica Surfaces and within Nanopore

Pressure-driven supercritical CO 2 transport through a silica nanochannel. Bing Liu , Xiaoqi Li , Chao Qi , Tingyi Mai , Kaiyun Zhan , Li Zhao , Yue S...
0 downloads 0 Views 2MB Size
Article pubs.acs.org/JPCC

Aqueous CO2 Solutions at Silica Surfaces and within Nanopore Environments. Insights from Isobaric−Isothermal Molecular Dynamics Ariel A. Chialvo,*,† Lukas Vlcek,† and David R. Cole‡ †

Chemical Sciences Division Geochemistry & Interfacial Sciences Group, Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831-6110, United States ‡ School of Earth Sciences and Department of Chemistry, The Ohio State University, Columbus, Ohio 43210, United States S Supporting Information *

ABSTRACT: We present a detailed molecular-based characterization via isobaric−isothermal molecular dynamics simulation of the microstructure and dynamics of water-rich aqueous CO2 solutions at silica surfaces and under extreme confinement between finite silica plates at state conditions relevant to geologic capture and sequestration of carbon dioxide. The study comprises three types of slit-pore plates to represent two extreme cases of surface polarity and a mismatched pair of plates to interrogate the fluid behavior at and confined between heterogeneous surfaces. We found layer formation of H2O and CO2 whose strength depends on the nature of the plate surface, i.e., stronger H2O layering at hydrophilic than at hydrophobic plates with simultaneous weaker water-mediated CO2/hydrophilic-surface interactions. We observed the opposite behavior with the hydrophobic plates in which the weaker water layering results from the CO2-mediated H2O/hydrophobic-surface interactions. Moreover, we illustrate how the interplay between these types of interactions and extreme fluid confinement, i.e., strong overlapping of interfacial structures, can induce a drying out of the pore environment whose immediate consequence is a significant CO2 concentration enhancement relative to that of the bulk environment. Finally, we assessed the effect of the nature of the plate surfaces on the translational diffusion coefficient of water. We found that this property changes monotonically at purely interfacial regions but nonmonotonically under confinement.

1. INTRODUCTION The capture and sequestration of CO2(CCS) in deep geological reservoirs have been considered as a potential approach to mitigate its release into the atmosphere and reduce its contribution to the greenhouse effect on climate change.1 The feasibility and safety of this process for long-term storage of CO2 depends on the low hydraulic permeability of the caprock and its ability to hold the brines plus CO2 in its porous structure, i.e., its interfacial and confinement properties. Interfacial and confined fluids exhibit microstructural, dynamical, and thermophysical behavior that differ dramatically from their bulk counterparts.2−10 This behavior results from the molecular asymmetry between fluid−fluid and solid−fluid interactions that manifests as inhomogeneous local fluid density distributions, where the overlapping of these interfacial structures can induce confined environments. The immediate consequence is the inherent inability of current macroscopic modeling approaches to capture the actual fluid−solid (aqueous CO2−caprock) and fluid−fluid interfacial mechanisms underlying the geological CO2 sequestration. The injection of CO2 into confined aqueous environments can perturb the equilibrium between ionic species, with the © 2012 American Chemical Society

potential for either mineral surface dissolution or carbonate precipitation. Secure storage of CO2 in the subsurface depends on retention involving hydrodynamic, solubility, and mineral trapping mechanisms. The engineering of more effective retention processes depends on our ability to manipulate the fluid−solid interfacial behavior, where the latter relies on our atomic-level interpretation/description of the interplay between the nature of the minerals, their porosity, fluid composition, and state conditions. In fact, we currently lack a fundamental understanding of the structural and dynamic behavior of aqueous CO2 in contact with caprock minerals, i.e., fluid environments at nonambient state conditions interacting with substrates characterized by a wide range of nanoporosity and hydrodynamic permeability. Examples of such rock formations include (a) natural sandstones involving a pore-size distribution in the range 10 Å−50 μm according to the small-angle neutron scattering analysis combined with backscatter SEM and NMR methods;11 (b) tight shale samples exhibiting unimodal 0.3−60 Received: January 6, 2012 Revised: June 1, 2012 Published: June 7, 2012 13904

dx.doi.org/10.1021/jp3001948 | J. Phys. Chem. C 2012, 116, 13904−13916

The Journal of Physical Chemistry C

Article

Figure 1. (a) View of the NPT tetragonal simulation box for the slit-pore configuration involving two silica plates immersed in a water-rich CO2 solution. Yellow rectangle highlights the region defining the μVT subsystem, (b) Schematic representation of the formation of a confined environment as a result of the overlapping of two approaching solid−fluid interfacial structures.

interactions. While explicit polarizable models can successfully handle the asymmetry of species polarities, the study of solubility phenomena in the H2O−CO2 system has typically been done in terms of nonpolarizable intermolecular potential models.16−19 Under these conditions, the accurate prediction of species solubility becomes more challenging, as we have recently demonstrated for the case of the SPC/E and EPM2 models,20 for which we developed a coupling parameter approach to adjust the unlike-pair interactions in the context of fluid−fluid phase equilibria at representative CCS state conditions. The systems under investigationH2O−CO2 in contact with mineral surfaces and under confinement within porous mineral substratesare obviously open, i.e., subject to constant temperature and pressure with a fluctuating number of fluid particles and resulting composition toward the equilibration of the species chemical potentials with the surroundings.21 The physicochemical behavior of these systems can be efficiently studied by a combination of NPzzT and μPzzT22 or by a μVT23 molecular simulation approach involving infinitely replicated solid surfaces. Alternatively, for the case of finite-size solid surfaces immersed in fluid systems (Figure 1), their behavior can be successfully studied by molecular dynamics in a globally isobaric−isothermal ensemble (NPT−MD) while locally as a grand canonical ensemble (μVT).24−26 Thus, the overarching goal of this study is to use reliable and precisely defined intermolecular models to provide a detailed molecular-based characterization of the microstructure and dynamics of aqueous CO2 solutions at silica surfaces and under

nm pore-size distributions based on XRD, SEM, and petrographical thin-section analysis;12 (c) natural silica minerals such as chalcedony, chert, quartz, and volcanic glass that exhibit 1.5− 3.5 Å pore-size distributions according to positron annihilation lifetime spectroscopy (PALS);13 and (d) sedimentary (equigranular quartzite) sands whose NMR and PALS analyses indicated bimodal 4−200 Å pore-size distributions.14 A typical shale caprock can be described as quartz−clay rock,15 where for modeling purposes the interrogation of the behavior of aqueous CO2 in contact with fine-grained quartz component can be analyzed in terms of a pure silica substrate proxy. Fractures and cracks in the caprock allow the migration of the aqueous CO2 into and through these caprock nanoenvironments, facilitate the chemical interactions (dissolution/precipitation) between the fluid and the solid substrates, and can potentially affect the substrate permeability. Therefore, it is essential to gain microscopic insights into the interfacial aqueous CO2−caprock phenomena, such as the wetting and adsorption for variable pore sizes and geometries as well as for pore wall chemistry (hydrophobicity/hydrophilicity) at the actual CCS state conditions. A key parameter underlying the CO 2 capture and sequestration is its aqueous solubility, not just at bulk conditions, but more importantly, within nanoconfined environments. From a simulation viewpoint, the proper description of H2O−CO2 solubility requires intermolecular potential models able to depict the solvation behavior of polar species in nonpolar environments (and vice versa) through an appropriate description of the corresponding unlike-pair 13905

dx.doi.org/10.1021/jp3001948 | J. Phys. Chem. C 2012, 116, 13904−13916

The Journal of Physical Chemistry C

Article

extreme confinement between finite silica plates. Specifically, we address some fundamental issues, including (a) how the degree of surface hydrophobicity (i.e., surface polarity) affects the interfacial structure, (b) how the overlapping of approaching interfacial structures affects the confined fluid composition (i.e., relative solubility), (c) how the hydrogenbond strength and tetrahedral nature of the hydrogen-bonded network change under mineral confinement, and (d) how the fluid mobility is affected by all of the above. To confront these issues, we perform extensive molecular dynamics simulations of atomistically explicit descriptions for the aqueous solutions and the silica plates according to the models and methodology presented in Section 2. In Section 3, we present and discuss the simulated microstructural and dynamical results in terms of the full characterization of the inhomogeneous environments. Finally, in Section 4 we discuss the significance of our findings and summarize their implications on the modeling of nanopore aqueous environments.

Table 1. Lennard-Jones Potential Parameters and Partial Charges for the H2O and CO2 Models and Their CrossInteractions iiinteraction

εii/k (K)

σii (Å)

qi (e)

Ow−Ow Oc−Oc C−C Ow−Oc

78.23 80.507 28.129

3.166 3.033 2.757

−0.8476 −0.3256 0.6512

ηij

ξij

0.973

1.0

Ow−C

0.973

1.0

Ow−Oc Ow−C

1.0171 0.9594

1.1351 1.4134

ref 28 29 29 case a in this work case a in this work case b 20 case b 20

invoked the second set of combining rules, i.e., case b in Table 1, to be used throughout the current investigation. To analyze the effect of the surface chemistry on the interfacial/confinement behavior, we consider two types of silica surfaces: namely, a hydroxylated (hydrophilic) and a nonhydroxylated (hydrophobic) surface.24 The hydroxylated surface is built by attaching a hydrogen atom to each surface oxygen atom to form polar Si−O−H groups carrying partial electrostatic charges,36 as described in Table 2, where the −H

2. MODELS AND SIMULATION APPROACH To study the behavior of water-rich aqueous CO2 solutions at the silica-fluid interfaces (free-standing plate configuration) and under confinement between silica plates (finite-size, slit-pore configuration), we performed isobaric−isothermal molecular dynamics simulations of aqueous CO2 solutions at the representative subsurface sequestration state conditions of P = 20 MPa and T = 318 K as well as a global mole fraction xCO2 ≃ 0.023 corresponding to the water-rich phase in liquid−liquid equilibria (LLE) at the above state conditions.27 For that purpose, we have chosen simple intermolecular potential models such as the rigid SPC/E for H2O28 and the EPM2 for CO229 because these models provide reasonably good descriptions of the behavior of pure water30−33 and pure carbon dioxide,34 respectively. For the H2O−CO2 cross-interactions, we have analyzed two sets of “adjusted” combining rules; namely: (a) Berthelot’s combining rule to describe the Lennard-Jones energy parameters in conjunction with an adjusted Lorentz combining rule to match the simulated equilibrium composition of the water-rich phase to the experimental value in the corresponding LLE and (b) the recently proposed (ξ,η)-perturbed and polarization-corrected Lorentz−Berthelot combining rules that are able to describe simultaneously the species solubility in the two fluid phases at equilibrium.20 Note also that there is no need to include reactive interactions in the description of the fluid phase equilibria of H2O−CO2 at these state conditions.35 Initially, we chose to adjust the Lennard-Jones size parameters for the H2O−CO2 interactions as a deviation, ηij, from the corresponding values given by the Lorentz rule; i.e., ηij = 2σij/(σii + σjj) while keeping the Berthelot rule ξij = εij/ (εiiεjj)0.5 = 1.0 for the two unlike pair interactions C−Ow and Oc−Ow, where the subindices c and w denote CO2 and H2O molecules, respectively. The resulting adjusted parameters from the isobaric−isothermal Gibbs ensemble Monte Carlo simulations at P = 20 MPa and T = 318 K are ηOwC = ηOwCc ≃ 0.973, as illustrated in Figure SI-1 of the Supporting Information (see case a in Table 1). Although this set of combining rules predicts the correct CO2 solubility in the H2O-rich phase at the state condition of interest, not achievable with the conventionally used Lorentz−Berthelot rules, it does not allow the simultaneous accurate prediction of the H2O solubility in the corresponding CO2-rich phase at LLE.20 For that purpose, we

Table 2. Lennard-Jones Potential Parameters and Partial Charges for the Silica Sites36,a ii-interaction

εii/k (K)

σii (Å)

Si−Si Os−Os(bulk) Os−Os(surface) Hs−Hs(surface)

64.20 78.04 78.04

3.795 3.154 3.154

qi (e) 0.31b −0.71b 0.40b

a Unlike-pair interactions parameters are based on Lorentz−Berthelot combining rules. bElectrostatic charges on Si−O−H units for hydrophilic surfaces.

site is constrained to move in a circle centered around the oxygen site with a fixed bond length, S OH = 1 Å, and a fixed bond angle, ∠Si−O−H = 109.47°, using a shake−rattle routine37 (see Figure 2d). The nonhydroxylated (hydrophobic) surface counterpart is produced by detaching the −H from the silanol group and reducing to zero all their electrostatic charges so that the surrounding fluid and the silica surface will interact solely via nonelectrostatics interactions. On the basis of these two types of plate surfaces, we built three types of slit-pore confinements: the hydrophobic−hydrophobic, the hydrophilic−hydrophilic, and the mismatched (or Janus) hydrophobic−hydrophilic configurations (see Figure 3) The isobaric−isothermal molecular dynamics simulations of the aqueous solutions involved N = 2048 molecules and were carried out according to our own implementation of a Nosé− Poincare algorithm38,39 for the integration of the Newton− Euler equations of motion within a tetragonal, Lz = 2Lx = 2Ly, simulation box subject to 3D periodic boundary conditions. These systems comprise the aqueous solutions where we place two silica plates separated by a distance 6 ≤ h (Å) ≤ 14 and approximately equidistant from the center of the tetragonal box. The interplate distance, h, is defined as that between the planes containing the hydrogen sites of the surface silanol groups, even for the case of the hydrophobic surfaces, where these planes will 13906

dx.doi.org/10.1021/jp3001948 | J. Phys. Chem. C 2012, 116, 13904−13916

The Journal of Physical Chemistry C

Article

Figure 2. Views of the (a) front and (b, c) sides of the hydrophilic silica plate, as well as (d) a close-up of the silanol group.

still be located at S OH sin(19.47°) above the plane containing the corresponding silanol’s oxygen sites. Each silica plate, built with CrystalMaker software,40 comprises four rigid SiO2 layers representing the (111) face of cristobalite with a total of 314 and 278 atoms, whose approximate dimensions are 20.3 Å × 14.7 Å × 8.9 Å and 20.3 Å × 14.7 Å × 8.3 Å for the hydroxylated and nonhydroxylated plates, respectively (Figure 3). These plates are always surrounded by the aqueous solutions because the isobaric− isothermal system volume oscillates around the average value so that we are able to analyze simultaneously the aqueous CO2 bulk solutions in equilibrium with the corresponding interfacial and confined environments; i.e., these interfacial and confined regions behave effectively as grand canonical systems. The initial simulations at the chosen global composition were started from random fluid configurations generated by the Packmol utility41 and equilibrated by at least 2.0 ns, followed by production runs of about 24−30 ns, using an integration time step of 3.0 fs. All interactions were truncated at a cutoff radius rc ≈ 4.5σOO, where the configurational energy and virial were corrected by the corresponding long-range contributions, and the electrostatic interactions were accounted for by an Ewald summation whose convergence parameters were chosen to obtain an error smaller than 5 × 10−5εOO for both the real and reciprocal spaces.42 For the analysis of the fluid local environments at the silica− aqueous interfaces and the characterization of the confinement behavior, we first determine the axial (perpendicular to the silica plates(s)) profiles for the species densities ρi(z) and the

relative orientation of the water and carbon dioxide molecules ρθ(z) as follows, Ρ(z) ≡

∑ Ρi Βc (zi , zi − 0.5Δ, zi + 0.5Δ) (LxLyΔ) i

(1)

where Pi denotes either δ(z−zi) or the angle θi, Lα represents the lateral dimension of the plate along the α-axis, Δ ∼ 0.025σOO is the bin-size, Bc(x, a, b) ≡ [Θ(x−a) − Θ(x−b)] is the “boxcar” function,43 δ(···) is the Dirac delta function, Θ(···) is the Heaviside function, and ⟨···⟩ indicates a time average over the simulation trajectory. For the determination of the relative molecular orientation and because we are dealing with linear (carbon dioxide) and nonlinear (water) molecules, these angles are given by θi = cos−1(μi·ẑ/|μi|) for water and θi = cos−1(δi·ẑ/ |δi|) for carbon dioxide, where μi is the dipole moment of the water molecule, ẑ is the outward unit vector normal to the plate surface, and δi is the vector direction along the carbon dioxide O−C−O bond. Then, we characterize the local water packing behavior in terms of three relevant structural descriptorsthe first coordination number, the hydrogen bonding strength, and the tetrahedral order parameterto aid in the interpretation of the resulting inhomogeneous fluid environment. For that purpose, we invoked the inhomogeneous version of the tetrahedral order parameter introduced by Chau and Hardwick,44 written in the normalized form suggested by Errington and Debenedetti,45 i.e., 13907

dx.doi.org/10.1021/jp3001948 | J. Phys. Chem. C 2012, 116, 13904−13916

The Journal of Physical Chemistry C

Article

Figure 3. Schematic of the three silica-plate surface configurations.

qT (z) ≡ 1 −

neighbors within a specified shell radius, rc, given by the location of the first valley of the corresponding radial distribution functions, gOwOw(r) and gOwC(r), respectively.20 To study the effect of mineral confinement on the fluid composition, we calculate the average species densities within the interplate region as a function of the interplate distance h, ⟨ρi(h)⟩ = ⟨ni(h)⟩ × mwi/Vslit(h), where ni(h) is the number of molecules of the species i within the confined space of volume Vslit(h), mwi is the corresponding molecular weight, and ⟨···⟩ denotes an ensemble average. Note that this definition of average species densities involves a certain degree of arbitrariness resulting from the choice of the corresponding confined (accessible) volume. Throughout this study, we will use the volume of the slit pore defined in terms of the plates' dimensions (see above) and the interplate distance, h; i.e., Vslit(h) Å3 = 20.3 Å × 14.7 Å × h Å. To complete the characterization of the fluid behavior in the inhomogeneous environment, we calculated the self-diffusion coefficient of water (for all environments) and carbon dioxide (for the case of hydrophobic extreme confined environment) according to the local in-plane mean-squared displacement approach proposed by Liu et al.48 and based on the solution of the anisotropic Smoluchowski equation as follows,

∑ Ρi Βc (zi , zi − 0.5Δ, zi + 0.5Δ) i

N (z ) (2)

for i such that (0 ≤ xi ≤ Lx, 0 ≤ yi ≤ Ly), where Pi = 3/8∑j =3 1 ∑4k=j+1(cos ϕijk − cos(109.47°))2, cos ϕijk = rij·rik/(rijrik), N(z) = Δ∑i δ(z−zi) Bc(zi, zi − 0.5Δ, zi + 0.5Δ), with a bin-size Δ = 0.25σOO, the double summation comprises the four nearest neighbors (j, k) of the water molecule located at the axial distance zi, and the angle brackets denote time average over the simulation trajectory. Therefore, this parameter spans from qT(z) = 0, i.e., uncorrelated 4-molecule configurations, to qT(z) = 1, representing perfectly tetrahedral configurations of water oxygens. For the characterization of the axial profile of the hydrogen bonding strength, i.e., the average number of hydrogen-bonded molecules to any water molecule located at an axial position z, we use the rβ-geometric definition46 based on the interpretation by Wernet et al.47 of X-ray absorption spectroscopy; i.e., given by a pair of water molecules whose oxygen sites are separated by a distance rOO(β) = 3.3 − 0.00044β2 (Å), where the angle β (degrees) ≡ ∠Hα − Oα···Oγ is defined by the intramolecular Oα−Hα bond vector on molecule α and the intermolecular Oα−Oγ vector formed with molecule γ. During the sorting of pair distances during the calculation of qT(z), we also determine the corresponding axial profiles of the water coordination numbers around species, i.e., nOwOw(z, rc) and nOwC(z, rc) on the basis of the number of water nearest

Τ 2α

⟨ΔS (τ )⟩{z} =

∑ N(t )−1 ∑ Ρi t=1

i

Βc (zi , zi − 0.5Δ, zi + 0.5Δ)

(3)

where Pi = ∑i∈{(t+τ)}[(xi(t + τ) − xi(t)) + (yi(t + τ) − yi(t))2]{z}, the subscript {z} indicates the location of the slab of thickness Δ = σOO perpendicular to the plane where the 2

13908

dx.doi.org/10.1021/jp3001948 | J. Phys. Chem. C 2012, 116, 13904−13916

The Journal of Physical Chemistry C

Article

diffusion takes place, {(t + τ)} represents the set of α-type molecules N(t + τ) that stay within the slab centered at z during the time interval t to t + τ, T is the total number of time steps involved in the averaging process, N(t) denotes the number of α-type molecules residing in the slab at time t, and ⟨···⟩ is the corresponding ensemble average. Because not all N(t) molecules will stay the entire time τ within the slab centered at z, the fraction of the original number N(t) that actually would contribute to the local diffusion process is given by the survival probability P(τ) = ∑t =T 1N(t + τ)/N(t); consequently, the “adapted” Einstein relation for the determination of axial dependence of the in-plane diffusion coefficient becomes D α (z) = lim ⟨ΔS α 2(τ )⟩{z} /4τ Ρ(τ ) τ →∞

Table 3. Global, Bulk, and Confined Fluid Densities of the Aqueous Systems at T = 318 K and P = 20 MPa system philic−philic

phobic− phobic

(4)

philic−phobic

For the application of eq 4, we first analyzed the axial behavior of the local density, ρi(z), and according to the location of the relevant peaks in both the interfacial and confined regions, we determined the corresponding in-plane species translational diffusivity. To enhance the sampling in the averaging of the in-plane mean squared displacement without extending the simulation length, we sample simultaneously several time-origins, t, and take advantage of the obvious axial symmetry of the plate configurations, except for the mismatched case when that symmetry is broken.

h (Å)

⟨ρ⟩ (g/cm3)a

⟨ρslit(h)⟩ (g/cm3)b

6 8 10 12 14 6

0.906 0.905 0.903 0.903 0.903 0.848

± ± ± ± ± ±

0.001 0.001 0.001 0.001 0.001 0.001

0.871 0.825 0.807 0.797 0.794 0.206

1.051 1.052 1.047 1.048 1.049 1.038

± ± ± ± ± ±

0.003 0.003 0.003 0.003 0.003 0.003

8 10 12 14 6 8 10 12 14

0.847 0.849 0.850 0.850 0.879 0.879 0.877 0.877 0.877

± ± ± ± ± ± ± ± ±

0.002 0.002 0.002 0.001 0.001 0.001 0.001 0.002 0.001

0.282 0.530 0.564 0.575 0.626 0.650 0.673 0.695 0.711

1.039 1.032 1.031 1.032 1.042 1.045 1.042 1.041 1.041

± ± ± ± ± ± ± ± ±

0.003 0.004 0.004 0.004 0.004 0.003 0.004 0.004 0.005

⟨ρbulk⟩ (g/cm3)c

NPT average global fluid density. bAverage confined-fluid density. Average axial fluid density away from the interfaces.

a c

3. DISCUSSION OF SIMULATION RESULTS Here, we address the characterization of the silica−aqueous CO2 interfacial environments as a function of the interplate separation, h, in terms of the contrasting microstructural and dynamical behavior between the interfacial and confined regions. Although confinement (i.e., in a slit pore in our case) involves at least two interfacial regions, we make the distinction between a (purely) interfacial and a confined region to highlight that the latter comprises actually the overlapping of two interfacial structures corresponding to the plates conforming the slit-pore configuration (see how the approaching of the left (L) and right (R) interfaces results in their overlapping with the formation of the confined fluid environment in Figure 1b). As we will illustrate below, the degree of interfacial overlapping (∼h−1) affects the inhomogeneous density distribution and the species solvation, as well as the hydrogen bonding, and depends strongly on the nature of the plate surface. 3.1. Interfacial and Confined Fluid Structure. In Table 3, we summarize the relevant fluid density information regarding the systems under investigation. In Figures 4 and 5, we display the comparison of the axial distribution functions for the H2O sites as a function of the interplate separation h for the hydrophilic and hydrophobic silica plates at the same global state conditions of P = 20 MPa and T = 318 K. The CO2 counterparts are shown in Figures 6a, b as axial local density distributions rather than distribution functions since the extremely low global solute concentration makes it difficult to determine the proper values for the normalizing bulk CO2 densities. These figures highlight the contrasting interfacial behavior between the external (purely interfacial) and the internal (confined) regions. In particular, and regardless of the nature of the plate surfaces, we observe the layer formation of H2O and CO2 within the first 10−15 Å from the plane containing the

Figure 4. Interplate dependence of the water sites axial distribution functions for interfacial and confined water-rich CO2 solutions involving hydrophilic silica plates at T = 318 K and P = 20 MPa .

silanol’s hydrogens on the silica surface, beyond which the axial structure is completely lost. The strength of the observed layering depends on the nature of the plate surface; i.e., H2O exhibits stronger layering when interfacing with hydrophilic than with hydrophobic plates, clearly depicted as higher peaks in the corresponding axial distributions (Figure 4). This behavior results from the favorable dipolar water-surface interactions when the aqueous CO2 systems interface with hydroxylated surfaces, leading to a simultaneous repulsion of the nonpolar CO2 away from the surface (Figure 6b). 13909

dx.doi.org/10.1021/jp3001948 | J. Phys. Chem. C 2012, 116, 13904−13916

The Journal of Physical Chemistry C

Article

i.e., it destroys the nearest layer observed in the hydrophilic case and weakens the next-nearest layers with increasing peak overlapping to the point that a drying out of the confined environment occurs as h → 6 Å (Figure 5). On the other hand, in comparison to H2O and as a consequence of its nonpolar nature, CO2 shows the opposite interfacial trend in that it interacts (a) rather weakly with the hydrophilic plates as a fully hydrated species, i.e., through water mediating layers (Figure 6b), but (b) significantly with the hydrophobic plates by pushing the water molecules away from the plate surfaces (Figure 6a). In fact, the preferential interaction of CO2 with the hydrophobic plate surface is clearly manifested in Figure 6a as single strong adsorption peaks, in contrast to the multiple peaks observed in the hydrophilic plates, i.e., Figure 6b. This preferential CO2 adsorption at the external and internal interfaces translates into local densities over 1 order of magnitude larger than those of the corresponding bulk counterparts. The implications of this behavior for the CCS process will be discussed below. An additional feature in these figures is the evolution of pronounced layering structure in the interplate region, i.e., the manifestation of confinement effect on the inhomogeneous density distribution of species, while the interfacial structure for the external regions remains mostly unaffected. In fact, further analysis of the data in Figures 4−6 indicates that the location (for both types of plates) and strength (especially for the hydrophilic plates) of the species distribution peaks are independent of the interplate separation, i.e., a behavior that suggests negligible plate-mediated fluid−fluid interactions. This can be unambiguously illustrated in Figures 7, 8, and 9, where

Figure 5. Interplate dependence of the water sites axial distribution functions for water-rich CO2 solutions interfacing with, and under confinement between, hydrophobic silica plates at T = 318 K and P = 20 MPa .

Figure 6. Axial distribution of the carbon dioxide (carbon site) axial density for water-rich CO2 solutions interfacing with, and under confinement between, (a) hydrophobic and (b) hydrophilic silica plates at T = 318 K and P = 20 MPa .

Figure 7. Shifted axial distribution of the water (oxygen site) axial distribution functions for the external interfacial region of water-rich CO2 solutions in contact with (a) hydrophobic and (b) hydrophilic silica plates at T = 318 K and P = 20 MPa .

Moreover, Figures 4−6 demonstrate how differently these solutions behave under increasing confinement between either hydrophobic or hydrophilic surface plates. On one hand, water layering under hydrophilic confinement exhibits four welldefined peaks (with a narrow quasi-bulk region) for h = 14 Å and decreases to three peaks after the strong and gradual additive overlapping of a pair of them as h → 6 Å (Figure 4). In contrast, the hydrophobic confinement attenuates the water layering due to obvious weaker substrate−water interactions;

we replot the distributions of Figures 4−6 by shifting the axial coordinates as either z → z + (0.5h + S z) for z < 0 or z → z + (0.5h + S z) for z > 0, when the slit pore is axially centered at z = 0 and where S z represents the thickness of the silica plate. Although a quick comparison with the available simulation results from the available literature (for a recent review see ref 13910

dx.doi.org/10.1021/jp3001948 | J. Phys. Chem. C 2012, 116, 13904−13916

The Journal of Physical Chemistry C

Article

Figures 7a and 9a indicate that the main H2O layers for the external interfaces occur at ∼2.5 Å as well as ∼5.5 Å for water and at ∼1.8 Å for carbon dioxide, from the corresponding (hydrophobic) nonhydroxylated surfaces. Moreover, on the basis of the peak location in the water-site distributions (Figures 4, 5), and the corresponding axial distribution of θi = cos−1(μi·ẑ/|μi|) for the interfacial water molecules (Figure 10a), the water molecules in the first

Figure 8. Shifted axial distribution of the water (hydrogen site) axial distribution functions for the external interfacial region of water-rich CO2 solutions in contact with (a) hydrophobic and (b) hydrophilic silica plates at T = 318 K and P = 20 MPa .

Figure 10. Axial distribution of the relative orientation of the (a) water molecules and (b) carbon dioxide molecules for water-rich CO2 solutions interfacing with, and under confinement between, hydrophobic and hydrophilic silica plates at T = 318 K and P = 20 MPa .

interfacial layer interacting with the hydroxylated surface exhibit a significant but distinctive change in the dipolar orientation, depending on whether they are in the external or confined region. In fact, we note that water exhibits a quasi-random axial 2O dipolar distribution away from the interfaces, i.e., θHrand ≈ 90°, but as it approaches the interfacial region, the water dipole starts gradually pointing away from the surface. As water reaches the first adsorption layer and due to its strong dipolar interaction with the surface, the dipolar orientation reverses its trend, i.e., pointing toward the hydroxylated surface. However, we highlight that for the case of the internal interfaces, there is a second reversal of the orientational trend that may be ascribed to the overlapping of interfacial structures in the confined region. This behavior is clearly reflected by the appearance of a shoulder in the water-oxygen peak at the external first water layer, a feature that is magnified in the confined region (Figure 4). In contrast, the corresponding dipole orientational distribution for the nonhydroxylated surface is almost featureless, showing a small orientational order at the first water layer, dipoles slightly tilted toward the hydrophobic surface, a distinctive behavior in comparison with the longrange decay toward orientational randomness in the case of the hydrophilic plate surfaces. In contrast to the water, carbon dioxide shows an oscillatory orientational preference around randomness for a linear −1 1/2 2 molecule, i.e., θCO ≈ 54° (Figure 10b), rand = cos (1/3) where the peaks follow closely the evolution of the corresponding local (C site) density distributions (Figures 6

Figure 9. Shifted axial distribution of the carbon dioxide (C site) axial distribution functions for the external interfacial region of water-rich CO2 solutions in contact with (a) hydrophobic and (b) hydrophilic silica plates at T = 318 K and P = 20 MPa .

49) suggests that the highly dilute CO2 has a rather week effect on the water layering, this observation is based on simulation results involving pure water under slightly larger confinement and rather different state conditions (e.g., lower temperature and higher pressure24) that portray the layering of pure water only at the internal interfaces. For example, according to Figures 7b and 9b, the water molecules lie preferentially in the external interfacial layers centered at ∼0.5, ∼3.0, and ∼5.5 Å, whereas the carbon dioxide molecules are located at ∼3.37 and ∼6.45 Å from the plane containing the silanol’s hydrogen sites in the external interfacial structures of (hydrophilic) hydroxylated surfaces. Similarly, 13911

dx.doi.org/10.1021/jp3001948 | J. Phys. Chem. C 2012, 116, 13904−13916

The Journal of Physical Chemistry C

Article

and 9). As carbon dioxide approaches the first adsorption layer, 2 the angle θ deviates from θCO rand, depending on the nature of the CO2 2 plate surface, i.e., θ > θrand and θ < θCO rand, for the hydrophobic and hydrophilic surface plates, respectively. Regardless of the nature of the plate surfaces, we observe a similar lack of h-effect on the axial profiles in the external interfaces for the average tetrahedral order parameter qT(z); number of hydrogen bonds per molecule, nHB(z); first coordination numbers, nOwOw(z); and nOwC(z), as illustrated in Figures 11 and 12 for the hydrophilic and hydrophobic plate

Figure 12. Axial profiles of (a) qT(z), (b) nHB(z), (c) nOwOw(z), and (d) nOwC(z) for the water-rich CO2 solutions interfacing with, and under confinement between, hydrophilic silica plates at T = 318 K and P = 20 MPa .

Figures7−9) In fact, the distortion of the confined hydrogenbonded network is clearly illustrated in Figures 11 and 12 for the hydrophobic and hydrophilic plates, respectively, where we observe the contrasting behavior of all three structural parameters (qT(z), nHB(z), and nOwα(z)) with decreasing interplate separation; i.e., (a) the gradual reduction of their magnitude for the hydrophilic plate surface in comparison with (b) their profound perturbation for the hydrophobic case. In particular, note that for the hydrophobic slit pore with h = 6 Å, the observed drying out of the confined environment translates into qT(z) ≈ 0, with a significant reduction (from the bulk value) in the hydrogen bonding strength ΔnHB(z) ≈ −1.0, followed by the corresponding weakening of the H2O coordination around itself and around the CO2, i.e., ΔnOwOw(z) ≈ −2.0 and ΔnOwC(z) ≈ −13.0. In the Supporting Information document, we provide a comparison between and short discussion of the confined behavior of aqueous CO2 and simple aqueous electrolyte solutions. The contrasting confinement effect described above translates into a rather remarkable fluid behavior with significant implications in the analysis of caprock reliability in CCS processes. In fact, the nature of the confining mineral surface can have a profound effect on the equilibrium composition of the pore fluid, as we show in Figure 13a−b, where we plot the pore-averaged species densities and the corresponding CO2 mole fraction for the aqueous solutions confined between hydrophilic and hydrophobic silica plates, respectively, as a function of the interplate distance, h. In particular, note that H2O and CO2 pore-averaged densities as well as the corresponding pore composition (mole fractions in this case) are weakly dependent on h when involving hydrophilic plates. However, for hydrophobic plates, we find a completely different behavior in which the species pore-average densities exhibit a

Figure 11. Axial profiles of (a) qT(z), (b) nHB(z), (c) nOwOw(z), and (d) nOwC(z) for the water-rich CO2 solutions interfacing with, and under confinement between, hydrophobic silica plates at T = 318 K and P = 20 MPa .

surfaces, respectively. The observed similarity should not be unexpected if we consider that these three quantities are all linked to the ability of water to preserve its hydrogen-bonded network. Moreover, all these axial profiles exhibit two common features, namely: (a) a steep decrease at the first adsorbed water layer centered at either ∼0.5 or ∼2.5 Å for the hydrophilic or hydrophobic surface, respectively (see Figures 7a, b and 8a, b), caused by the spatial constraint imposed by the plate surface, resulting in a disruption of the hydrogen bonding whose effect depends on the magnitude of the surface dipole moment (i.e., hydrophobic versus hydrophilic surfaces); and (b) an asymptotic approach to the expected properties of bulk fluid environments for z > 15 Å away from the plane containing the silanol’s hydrogen sites. In contrast to the fluid behavior in the external interfacial regions, the structure of the confined environments undergos significant changes caused by the interpenetration (overlapping) of the pair of solid−fluid interfacial structures with decreasing interplate distance h (see Figures 4−6), i.e., this “confinement” effect becomes stronger with decreasing h, regardless of the nature of the plate surface, yet with a negligible impact on the fluid−solid exclude volume (e.g., see 13912

dx.doi.org/10.1021/jp3001948 | J. Phys. Chem. C 2012, 116, 13904−13916

The Journal of Physical Chemistry C

Article

Figure 14. Axial behavior of the in-plane translational diffusion coefficient of water in water-rich CO2 solutions interfacing with, and under confinement between, hydrophilic silica plates at T = 318 K and P = 20 MPa .

Figure 13. Interplate dependence of the pore-averaged species densities and corresponding CO2 mole fraction for the aqueous solutions confined between (a) hydrophobic and (b) hydrophilic silica plates at T = 318 K and P = 20 MPa .

strong h dependency, but in opposite directions, leading to the occurrence of a drying out of the pore environment and consequent enhancement of the pore CO2 concentration relative to that of the bulk with increasing confinement. 3.2. Interfacial and Confined Fluid Mobility. In Figures 14 and 15, we display the calculated in-plane translational diffusion coefficient of water, DH∥ 2O(z), at specific axial positions with respect to the mineral surface, i.e., the location of relevant water density peaks, for the three interplate distances and the two types of plate surfaces. We display only the left half of the profiles for the sake of clarity, since the systems are axially symmetric along z = 0. Regardless of the type of plate surface, we find that DH∥ 2O(z) decreases monotonically with a decreasing axial distance from the plate surface in the external interfacial regions, a behavior that is consistent with previous simulation studies involving pure water interfacing atomistically detailed silica surfaces.36,50 Note, however, that within the confined region, DH∥ 2O(z) exhibits a nonmonotonic behavior, especially for the nanopores comprising hydrophobic plates, in which the nanoconfined environment undergoes large density and compositional changes (see Figures 5 and 13b). In fact, for h = 6 Å, the average fluid density between hydrophobic plates is ρslit ≃ 0.2 g/cc (see Table 3) with a mole fraction xCO2 ≃ 0.88, and the main CO2 peak is characterized by a local (slab) density ρslab ≈ 0.7 g/cc and diffusivity DH∥ 2O(z = 0) ≃ 2.8 × 10−9 m2/s, in comparison with the experimental value DCO2(318 K) ≃ 3.2 × 10−9 m2/s.51

Figure 15. Axial behavior of the in-plane translational diffusion coefficient of water in water-rich CO2 solutions interfacing with, and under confinement between, hydrophobic silica plates at T = 318 K and P = 20 MPa .

hydroxylated or nonhydroxylated surfaces as representations of the two extreme of silica surface polarity,50 i.e., hydrophilic and hydrophobic, respectively. Actual caprock’s nanopores, however, consist of a variety of surface types, including partially (patchy or patterned) and fully mismatched surfaces that will

4. CONCLUDING REMARKS In the previous sections, we have analyzed the behavior of aqueous CO2 interfacing with, and under confinement of, axially symmetric silica plates described in terms of either 13913

dx.doi.org/10.1021/jp3001948 | J. Phys. Chem. C 2012, 116, 13904−13916

The Journal of Physical Chemistry C

Article

contribute differently to the interfacial phenomena under consideration. Therefore, considering that the quality of the plate surfaces in conjunction with severe fluid confinement can have significant effects on the equilibrium and dynamical properties of the confined environments (as illustrated in Figures 13b and 15), it appears appropriate to check how surface mismatch (as depicted by type c in Figure 3) might affect those properties. In Figures 16, 17, and 18, we display the profiles for water sites axial distribution functions, carbon dioxide (C-site) axial

Figure 18. Axial profiles of (a) qT(z), (b) nHB (z), (c) nOwOw(z), and (d) nOwC(z) for the water-rich CO2 solutions interfacing with, and under confinement between, mismatched silica plates at T = 318 K and P = 20 MPa .

(hydrophobic) slit-pore counterparts. In contrast, the corresponding internal (confined) environments exhibit axial profiles that transition from those characterizing the hydrophilic plates (left) to those of the hydrophobic plates (right). This evolution becomes more evident when the mismatched profiles are plotted simultaneously with those of the hydrophilic and the hydrophobic slit pores, as done in Figures 19 and SI-2ab (see Supporting Information), where we can clearly observe that while for z < 0, the “mismatched” profiles approach those of the hydrophilic pores, for z > 0, they approach their hydrophobic counterparts. Moreover, if we compare Figure 20 with Figure 13a and b, we find that the pore-averaged fluid properties confined within the mismatched slit pore become “average” behavior between those observed for the hydrophilic and hydrophobic slit pores. Note that these hydrophobic and mismatched surfaces are typically found during the formation of new fractures and cracks in the caprock substrates, as well as during the drying out of wet surfaces in contact with supercritical CO2 environments.52 The implications of this behavior are manifold, including (a) the slit-pore fluid environment (i.e., density and composition) cannot be represented by that of a bulk counterpart at the prevailing state conditions when either interpreting or modeling the process at a macroscopic level, (b) the chemical processes occurring in these nanopore aqueous environments depend on the nature of the confining mineral surfaces, (c) the fluid environment in contact with the mineral surface can be considerably different from either that a few molecular diameters away or that of the corresponding bulk, and (d) the competitive coadsorption of CO2 impurities at the mineral surfaces can radically modify the interfacial equilibrium

Figure 16. Interplate dependence of the water sites axial distribution functions for interfacial and confined water-rich CO2 solutions involving mismatched silica plates at T = 318 K and P = 20 MPa .

Figure 17. Axial distribution of the carbon dioxide (carbon site) axial density for water-rich CO2 solutions interfacing with, and under confinement between, mismatched silica plates at T = 318 K and P = 20 MPa .

density distributions, and the corresponding microstructural descriptors (qT(z), nHB(z), and nOwα(z)) resulting from the simulation of the aqueous CO2 for a mismatched slit-pore configuration comprising a hydrophilic plate (left) and a hydrophobic plate (right), respectively. A careful analysis of these results indicates that the left (right) external interfacial structures behave as those of the corresponding hydrophilic 13914

dx.doi.org/10.1021/jp3001948 | J. Phys. Chem. C 2012, 116, 13904−13916

The Journal of Physical Chemistry C



Article

AUTHOR INFORMATION

Corresponding Author

*Fax: 865-574-4961. E-mail: [email protected]). Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS Support for this work comes from the US Department of Energy through the Lawrence Berkeley National Laboratory “Center for Nanoscale Control of Geologic CO2” (FWP ERKCC67) under contract DE-AC05-00OR22725 to Oak Ridge National Laboratory, managed and operated by UTBattelle, LLC.



Figure 19. Comparison of the behavior of qT(z), (b) nHB(z), (c) nOwOw(z), and (d) nOwC(z) between the mismatched, the hydrophobic and the hydrophilic slit-pore configurations for h = 10 Å.

Figure 20. Interplate dependence of the pore-averaged species densities and corresponding CO2 mole fraction for the aqueous solutions confined between mismatched silica plates.

and dynamical behavior and, consequently, their chemical stability.53



REFERENCES

(1) Carbon Dioxide Capture and Storage: Special Report of the Intergovernmental Panel on Climate Change; Metz, B., Davidson, O., de Coninck, H., Loos, M., Meyer, L., Eds.; Cambridge University Press: New York, 2005. (2) Scodinu, A.; Fourkas, J. T. J. Phys. Chem. B 2002, 106, 10292− 10295. (3) Sofos, F. D.; Karakasidis, T. E.; Liakopoulos, A. Phys. Rev. E 2009, 79. (4) Alba-Simionesco, C.; Coasne, B.; Dosseh, G.; Dudziak, G.; Gubbins, K. E.; Radhakrishnan, R.; Sliwinska-Bartkowiak, M. J. Phys.: Condens. Matter 2006, 18, R15−R68. (5) Verdaguer, A.; Sacha, G. M.; Bluhm, H.; Salmeron, M. Chem. Rev. 2006, 106, 1478−1510. (6) Mancinelli, R.; Bruni, F.; Ricci, M. A. J. Phys. Chem. Lett. 2010, 1, 1277−1282. (7) Mamontov, E.; Wesolowski, D. J.; Vlcek, L.; Cummings, P. T.; Rosenqvist, J.; Wang, W.; Cole, D. R. J. Phys. Chem. C 2008, 112, 12334−12341. (8) Malani, A.; Ayappa, K. G.; Murad, S. J. Phys. Chem. B 2009, 113, 13825−13839. (9) Castrillon, S. R. V.; Giovambattista, N.; Aksay, I. A.; Debenedetti, P. G. J. Phys. Chem. B 2009, 113, 7973−7976. (10) Argyris, D.; Cole, D. R.; Striolo, A. Langmuir 2009, 25, 8025− 8035. (11) Radlinski, A. P.; Ioannidis, M. A.; Hinde, A. L.; Hainbuchner, M.; Baron, M.; Rauch, H.; Kline, S. R. J. Colloid Interface Sci. 2004, 274, 607−612. (12) Katsube, T. J.; Williamson, M. A. Clay Miner. 1994, 29, 451− 461. (13) Yoshizawa, K.; Sato, K.; Murakami, H.; Shikazono, N.; Fujimoto, K.; Nakata, M. Open-nano pores in natural minerals studied by positron lifetime spectroscopy. In Positron and Positronium Chemistry; Wang, S. J., Chen, Z. Q., Wang, B., Jean, Y. C., Eds., 2009; Vol. 607; pp 189−191. (14) Chesta, M. A.; Ramia, M. E.; Jeandrevin, S.; Martin, C. A. Appl. Phys. A 2009, 97, 301−307. (15) Cole, D. R.; Chialvo, A. A.; Rother, G.; Vlcek, L.; Cummings, P. T. Philos. Mag. 2010, 90, 2339−2363. (16) Lisal, M.; Smith, W. R.; Aim, K. Fluid Phase Equilib. 2004, 226, 161−172. (17) da Rocha, S. R. P.; Johnston, K. P.; Westacott, R. E.; Rossky, P. J. J. Phys. Chem. B 2001, 105, 12092−12104. (18) Vorholz, J.; Harismiadis, V. I.; Rumpf, B.; Panagiotopoulos, A. Z.; Maurer, G. Fluid Phase Equilib. 2000, 170, 203−234. (19) Duan, Z. H.; Zhang, Z. G. Geochim. Cosmochim. Acta 2006, 70, 2311−2324. (20) Vlcek, L.; Chialvo, A. A.; Cole, D. R. J. Phys. Chem. B 2011, 115, 8775−8784. (21) O’ Connell, J. P.; Haile, J. M. Thermodynamics: Fundamentals for Applications; Cambridge University Press: New York, 2005. (22) Odriozola, G.; Aguilar, J. F.; Lopez-Lemus, J. J. Chem. Phys. 2004, 121, 4266−4275.

ASSOCIATED CONTENT

S Supporting Information *

A set of seven additional figures to illustrate (a) the effect of deviations of the Lorentz combining rule on the equilibrium CO2 mole fraction in the water-rich phase from Gibbs ensemble Monte Carlo simulations, (b) the water-sites axial distributions for slit pores involving mismatched silica plates, and (c) the surface density distributions of water oxygen as well as carbon dioxide oxygen sites interfacing either the hydrophilic or the hydrophobic silica plate. In addition, this document includes a discussion regarding the comparison between the confined behavior of aqueous CO2 and simple aqueous electrolyte solutions. This material is available free of charge via the Internet at http://pubs.acs.org. 13915

dx.doi.org/10.1021/jp3001948 | J. Phys. Chem. C 2012, 116, 13904−13916

The Journal of Physical Chemistry C

Article

(23) Botan, A.; Rotenberg, B.; Marry, V.; Turq, P.; Noetinger, B. J. Phys. Chem. C 2010, 114, 14962−14969. (24) Giovambattista, N.; Rossky, P. J.; Debenedetti, P. G. Phys. Rev. E 2006, 73. (25) Rodriguez, J.; Elola, M. D.; Laria, D. J. Phys. Chem. B 2009, 113, 12744−12749. (26) Chialvo, A. A.; Cummings, P. T. J. Phys. Chem. A 2011, 115, 5918−5927. (27) Spycher, N.; Pruess, K.; Ennis-King, J. Geochim. Cosmochim. Acta 2003, 67, 3015−3031. (28) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269−6271. (29) Harris, J. G.; Yung, K. H. J. Phys. Chem. 1995, 99, 12021−12024. (30) Errington, J. R.; Panagiotopoulos, A. Z. J. Phys. Chem. B 1998, 102, 7470−7475. (31) Guillot, B. J. Mol. Liq. 2002, 101, 219−260. (32) Gonzalez, M. A.; Abascal, J. L. F. J. Chem. Phys. 2010, 132, 096101. (33) Hayward, T. M.; Svishchev, I. M. Fluid Phase Equilib. 2001, 182, 65−73. (34) Nieto-Draghi, C.; de Bruin, T.; Perez-Pellitero, J.; Avalos, J. B.; Mackie, A. D. J. Chem. Phys. 2007, 126, 064509. (35) The reactions between CO2 and H2O play a crucial role in many natural and technological processes, but they are not significantly involved in all thermodynamic aspects of H2O−CO2 mixtures. For example, the bulk fluid phase equilibria or relative composition of H2O and CO2 at solid−fluid interfaces is dominated by the overwhelmingly larger concentrations of molecular aqueous carbon dioxide than that of carbonic acid (H2CO3). The constant of equilibrium for the formation of carbonic acid at 315.8 K is ∼1.7 × 10−3 [Wang, X.; Conway, W.; Burns, R.; McCann, N.; Maeder, M. J. Phys. Chem. A 2010, 114, 1734− 1740 ], and therefore, there would be only about 1 H2CO3 molecule per 588 CO2 molecules in aqueous solution. Although the equilibrium constant varies with temperature and pressure, the percentage of H2CO3 molecules in solution is always insignificant in the context of our simulations(there are only 48 CO2 molecules dissolved in our aqueous environment). We also note that the rate of H2O + CO2 → H2CO3 reaction is extremely low (4.2 × 10−3 L/mol/s, correspondingto CO2 molecule half-life of about 3 s; i.e., the probabilityof a CO2 molecule converted to H2CO3 during our simulation run is about 10−8), and therefore, we cannot expect to observe a statistically significant number of CO2 molecules converted to H2CO3 during our simulation run. From this perspective, our system can be safely and successfully described as nonreactive; consequently, our treatment is appropriate. The study of reactive interactions between fluid molecules and surfaces deserves proper attention,but is beyond the scope of our present investigation. (36) Lee, S. H.; Rossky, P. J. J. Chem. Phys. 1994, 100, 3334−3345. (37) Leimkuhler, B. J.; Reich, S. Simulating Hamiltonian Dynamics; Cambridge University Press: Cambridge, 2005. (38) Nose, S. J. Phys. Soc. Jpn. 2001, 70, 75−77. (39) Okumura, H.; Itoh, S. G.; Okamoto, Y. J. Chem. Phys. 2007, 126, 084103. (40) CrystalMaker, 2009; http://www.crystalmaker.com/index.html (accessed December 2011). (41) Martinez, J. M.; Martinez, L. J. Comput. Chem. 2003, 24, 819− 825. (42) Fincham, D. Mol. Simul. 1994, 13, 1−9. (43) von Seggern, D. Standard Curves and Surfaces; CRC Press: Boca Raton, 1993. (44) Chau, P. L.; Hardwick, A. J. Mol. Phys. 1998, 93, 511−518. (45) Errington, J. R.; Debenedetti, P. G. Nature 2001, 409, 318−321. (46) Kumar, R.; Schmidt, J. R.; Skinner, J. L. J. Chem. Phys. 2007, 126, 204107. (47) Wernet, P.; Nordlund, D.; Bergmann, U.; Cavalleri, M.; Odelius, M.; Ogasawara, H.; Naslund, L. A.; Hirsch, T. K.; Ojamae, L.; Glatzel, P.; Pettersson, L. G. M.; Nilsson, A. Science 2004, 304, 995−999. (48) Liu, P.; Harder, E.; Berne, B. J. J. Phys. Chem. B 2004, 108, 6595−6602.

(49) Striolo, A. Adsorp. Sci. Technol. 2011, 29, 211−258. (50) Castrillon, S. R. V.; Giovambattista, N.; Aksay, I. A.; Debenedetti, P. G. J. Phys. Chem. B 2009, 113, 1438−1446. (51) Tamimi, A.; Rinker, E. B.; Sandall, O. C. J. Chem. Eng. Data 1994, 39, 330−332. (52) Gaus, I. Int. J. Greenhouse Gas Control 2010, 4, 73−89. (53) Chialvo, A. A.; Vlcek, L.; Cole, D. R. CO2 Confined Environments between Mismatched Silica Surfaces and Containing SO2/H2S Contaminants; Lawrence Berkeley National Laboratory: Berkeley, CA, 2011.

13916

dx.doi.org/10.1021/jp3001948 | J. Phys. Chem. C 2012, 116, 13904−13916