Article pubs.acs.org/JPCB

Carbon Dioxide Hydrate Phase Equilibrium and Cage Occupancy Calculations Using Ab Initio Intermolecular Potentials Srinath C. Velaga†,‡ and Brian J. Anderson*,†,‡ †

National Energy Technological Laboratory, Morgantown, West Virginia 26506, United States Department of Chemical Engineering, West Virginia University, Morgantown, West Virginia 26506, United States

‡

ABSTRACT: Gas hydrate deposits are receiving increased attention as potential locations for CO2 sequestration, with CO2 replacing the methane that is recovered as an energy source. In this scenario, it is very important to correctly characterize the cage occupancies of CO2 to correctly assess the sequestration potential as well as the methane recoverability. In order to predict accurate cage occupancies, the guest− host interaction potential must be represented properly. Earlier, these potential parameters were obtained by fitting to experimental equilibrium data and these fitted parameters do not match with those obtained by second virial coefficient or gas viscosity data. Ab initio quantum mechanical calculations provide an independent means to directly obtain accurate intermolecular potentials. A potential energy surface (PES) between H2O and CO2 was computed at the MP2/aug-cc-pVTZ level and corrected for basis set superposition error (BSSE), an error caused due to the lower basis set, by using the half counterpoise method. Intermolecular potentials were obtained by fitting Exponential-6 and Lennard-Jones 6-12 models to the ab initio PES, correcting for many-body interactions. We denoted this model as the “VAS” model. Reference parameters for structure I carbon dioxide hydrate were calculated using the VAS model (site−site ab initio intermolecular potentials) as Δμ0w = 1206 ± 2 J/mol and ΔH0w = 1260 ± 12 J/mol. With these reference parameters and the VAS model, pure CO2 hydrate equilibrium pressure was predicted with an average absolute deviation of less than 3.2% from the experimental data. Predictions of the small cage occupancy ranged from 32 to 51%, and the large cage is more than 98% occupied. The intermolecular potentials were also tested by calculating the pure CO2 density and diffusion of CO2 in water using molecular dynamics simulations.

■

INTRODUCTION Clathrate (gas) hydrates are non-stoichiometric crystalline inclusion compounds consisting of water host cages and small guest molecules. The water molecule network is built by hydrogen bonds, and polyhedral cavities are formed. These cavities are stabilized by the encapsulated guest molecules. Large deposits of CH4 hydrates are found in the permafrost and on the continental margins; these could be potentially a source of energy if recovered in an economical way.1 Hydrates have also been considered as a possible solution to the CO2 problem. The idea of sequestrating the carbon dioxide in the oceanic sediments as CO2 hydrates to hold the increase in greenhouse gas in the atmosphere has been proposed.2 It is also suggested that CO2 can be used for the recovery of the CH4 gas from the naturally occurring CH4 hydrate by replacing CH4 hydrate with CO2.3 This replacement process (swapping process) allows both CH4 recovery and CO2 sequestration, as the formation conditions for the CO2 hydrate are known to be more stable than that of CH4 hydrate.4 Therefore, the thermodynamics of gas hydrates are important to predict the swapping process between two gaseous guests. Also, the cage occupancies must be correctly characterized to assess the sequestration potential and also the methane recoverability with the swapping process. Natural gas hydrates are known to occur in three structural forms: structure I (sI), structure II (sII), and structure H (sH). Pure carbon dioxide and methane each form structure I hydrate. © 2013 American Chemical Society

A basic thermodynamic model for the phase equilibrium of the gas hydrates was developed by van der Waals and Platteeuw5 (vdWP) from statistical thermodynamics along with the Lennard-Jones and Devonshire (LJD) spherical cell approximation.6,7 This model was modified by Parrish and Praushnitz8 using the Kihara intermolecular potential to describe the guest and surrounding water molecule interactions, and they also extended the model for mixed gas hydrates. A significant amount of studies have been done over the years to improve the original vdWP model. Studies by Holder and co-workers9−11 and by Tester and co-workers12−14 showed that the addition of terms to account for the contribution of second and subsequent water shells to the potential energy of the guest−host interactions was large enough to significantly influence the phase equilibrium predictions. The intermolecular potential parameters between the guest and host water molecules, often the Kihara potentials, are commonly fitted to experimental phase equilibrium data for pure or mixed hydrate systems when there are data.1 These intermolecular potentials are then used to predict thermodynamic equilibrium conditions. These fitted parameters work very well in the experimental range but fail when extended Received: October 17, 2013 Revised: December 9, 2013 Published: December 13, 2013 577

dx.doi.org/10.1021/jp410306v | J. Phys. Chem. B 2014, 118, 577−589

The Journal of Physical Chemistry B

Article

outside the range of fit, including when applied to mixed hydrate thermodynamics.14−16 In the widely used form of the vdWP model, the guest−guest interactions were neglected but it was shown that these interactions are significant enough to the contribution of the Langmuir constant and can influence the phase equilibrium calculations.17,18 Also, the Langmuir constant was evaluated using the spherical smoothed cell approximation, which introduces inaccuracies in evaluating the configurational integral.13,19 Sparks et al.13 showed that the geometrically explicit water positions of the cage structure can be taken into account by computing the full configurational integral through the use of the Monte Carlo technique or by the multi-interval Gauss−Legendre quadrature formula. The composition of the hydrate at a given temperature and pressure must be known in order to fully understand the thermodynamics of the gas hydrate formation and decomposition. The cage occupancy, i.e., the fraction of each cage occupied by the guest molecule, reported for CO2 hydrate using experimental techniques is nearly 100% for the large cage, but there is a large discrepancy in measurements and predictions of the small cage occupancy.20,21 In all of the experimental measurements20−22 except those by Ripmeester and Ratcliff,23 the CO2 hydrate samples prepared for determining the cage occupancies and hydration numbers were well above the equilibrium pressures, and these higher pressures during the synthesis are likely to produce higher occupancies. Ripmeester and Ractliff,23 with hydrate prepared near equilibrium conditions, report a lower limit to the hydration number of 7.0 for CO2 hydrate and a measurement of small cage to large cage ratio (θS/θL) equal to 0.32. The cage occupancies reported by the published theoretical models1,24,25 overestimated the cage occupancy for CO2 hydrate; mainly the small cage occupancies compare to the data of Ripmeester and Ractliff. Physically based intermolecular potentials are needed to more accurately predict phase equilibria and cage occupancies and for use in molecular simulation studies of hydrates. Cao et al.,26 Klauda and Sandler,24 and Anderson et al.27 computed guest− host intermolecular potentials from ab initio quantum mechanical calculations. With these ab initio intermolecular potentials, Anderson et al.27 computed the Langmuir constant and then calculated accurate phase equilibrium and cage occupancies for methane hydrate. Ab initio quantum mechanical calculations provide an independent means to directly obtain accurate guest− host intermolecular potentials. In this work, we have determined the effective interaction energies between the CO2 guest molecule and the water host molecules by developing guest−host pair potentials using an ab initio potential energy surface. These ab initio intermolecular potentials were then used to calculate the Langmuir constants by including the contributions of interactions between the CO2 guest and the host molecules from the first water shell to the fourth water shell and also the contributions from guest−guest interactions. A full configurational integral was evaluated in the Langmuir constant calculation. Using these Langmuir constants, the phase equilibrium and cage occupancy of the CO2 hydrate were predicted. To test the developed intermolecular potentials, molecular dynamics simulations were also performed using the developed potential model to calculate the densities of pure CO2 liquid and vapor, and the diffusion coefficient of CO2 in water.

μwH (T , P) = μwL, α (T , P)

(1)

where μHw (T, P) is the chemical potential of water (w) in the hydrate phase (H) and μL,α w (T, P) is the chemical potential of water in the water rich (L) or ice phase (α) at a given temperature, T, and pressure, P. The water rich liquid or ice phase is dependent on whether the temperature is above the freezing point or not. Using μβw, the chemical potential of a hypothetical empty hydrate lattice with no cages occupied by guest molecules, the condition for equilibrium can be written as in eq 2. Δμwβ − H = Δμwβ − L, α

(2)

where Δμβ−H ≡ μββ − μHw and Δμβ−L,α ≡ μβw − μL,α w w w . The statistical thermodynamics model proposed by van der Waals and Platteeuw5 for the difference in chemical potential between the empty lattice (β) and fully filled hydrate lattice is 2

Δμwβ− H = −RT ∑ νj ln(1 − j=1

∑ θji) (3)

i

where vj is the number of j-type cavities per water molecule and θji is the fractional occupancy of j-type cavities with i-type guest molecules

θji =

Cjifi 1 + ∑i Cjifi

(4)

where f i is the fugacity of gas component i, calculated using the Peng−Robinson equation of state.28 Cji is the temperaturedependent Langmuir constant species i in cavity j. The Langmuir constant for a nonspherical and nonlinear guest species i is defined as Cji =

Zji kT

=

1 8π 2kT

⎛

⎞

ϕ, α , β , γ) ⎟ ∫V exp⎜⎝− Φ(r , θ , kT ⎠

r 2 sin θ dr dθ dϕ dα sin β dβ dγ

(5)

where Zji is the configurational integral. The total interaction potential, Φ, is pairwise additive and includes the guest−host as well as the guest−guest interactions. Φ is a function of both the angular position of the center of mass of the molecule and the Euler angles (α, β, γ) for the orientation. Carbon dioxide is a linear molecule, so there are only five degrees of freedom. The Langmuir constant for a linear molecule of species i in cavity j is given as Cji =

Zji kT

=

1 4πkT

⎛

⎞

, ϕ, β , γ) ⎟ ∫V exp⎜⎝− Φ(r , θkT ⎠

r 2 sin θ dr dθ dϕ sin β dβ dγ

(6)

In the Langmuir constant calculation, the contributions of the interactions between guest and host molecules from first water shell to fourth water shell were included. The coordinates of the host water molecules were obtained from the crystallographic data.29,30 A full five-dimensional integral, the guest molecule orientation in the cage, was evaluated in calculating the Langmuir constant. The guest−guest interaction energy contribution to the Langmuir constant was ignored in previous thermodynamics models. However, it was shown that the guest−guest interactions can greatly affect the total interaction energy.18,31 The Langmuir constant is the combination of interaction energies from

■

THERMODYNAMIC MODEL The criterion for the phase equilibrium is the equality of chemical potentials (μ) of each component in the coexisting phases. At equilibrium, 578

dx.doi.org/10.1021/jp410306v | J. Phys. Chem. B 2014, 118, 577−589

The Journal of Physical Chemistry B

Article

guest−host (gh) interactions and guest−guest (gg) interactions. The Langmuir constant is then given as31 ⎛ Φgg ⎞ 1 ⎟ Cji = exp⎜ − ⎝ kT ⎠ 4πkT

⎛

∫V exp⎜⎝− Φ

gh

ΔH wβ − L, α = ΔH w0(T0) +

(r , θ , ϕ , β , γ ) ⎞ ⎟ kT ⎠

ΔCp = ΔC p0(T0) + b(T − T0)

2

r sin θ dr dθ dϕ sin β dβ dγ

(7)

gg Φgg = θjilike Φgg like + θji unlike Φ unlike

where the guest−guest interaction energy (Φgg) in cavity j is split into contributions from the like cages and unlike cages. Φgg like is the interaction energy of a CO2 molecule in a cavity of the same type. Φgg unlike is the interaction energy of a CO2 molecule in unlike cavities. The guest−guest interaction energy was calculated using the CO2−CO2 intermolecular potentials (VAS model) given in Table 3. The positions of guest CO2 molecules were kept at the center of the cages, and the orientation of the guest molecule positions in the cages were obtained from a Molecular Dynamics (2 × 2 × 2 unit cell of sI CO2 hydrate) energy minimization step using the steepest decent algorithm by the program GROMACS.32 An iterative procedure was used to determine the cage occupancy; first, the cage occupancy is calculated by assuming no guest−guest interactions and then iterating over eqs 4, 7, and 8 until there is a negligible change in the cage occupancy. The guest−host multidimensional configurational integral in eq 7 was calculated using 10-point multi-interval Gauss−Legendre quardrature formulas.13,33 The chemical potential difference, Δμβ−L,α , between the hypow thetical empty hydrate lattice and water in the hydrate phase is given by Holder et al.34 as

kT

kT +

∫0

P

∫T

T

0

ΔH wβ − L, α kT 2

ΔVwβ − L, α dP − ln a w kT

Table 1. Reference Properties between the Empty Hydrate Lattice and Fluid Phase (Liquid Water or Ice) for Structure I structure I 3 ΔVβ−α w (m /mol) 3 ΔVL−α (m /mol) w L−α ΔHw (J/mol) ΔCβ−α (J/mol·K) p ΔCβ−L (J/mol·K) p β−α 2

b (J/mol·K ) bβ−L (J/mol·K2) Δμ0w (J/mol) ΔH0w (J/mol)

dT

(9)

3.0 × 10 −1.598 × 10−6 6009.5 0.565 −37.32 0.002 0.179 1206 1260

reference 14 14 15 34 34 34 34 this work this work

■

CALCULATION OF AB INITIO INTERMOLECULAR POTENTIALS An ab initio quantum mechanical calculation provides an independent method to directly obtain intermolecular potentials which can be used in gas hydrate modeling. The exact value of the system energy and other properties can be obtained by solving the time-independent Schrödinger equation described below (13) H Ψ = EΨ where H is the Hamiltonian operator for the system of nuclei and electrons, E is the energy of the system, and Ψ is the electron

f wL fw

−6

have determined the two critical thermodynamic reference parameters described in eqs 9 and 11, Δμ0w and ΔH0w. Several methods, both experimental and analytical, have been adopted in the past to determine the reference parameters. The reference parameters, Δμ0w and ΔH0w, given by earlier researchers for structure I are given in Table 1 of Anderson et al.27 Holder et al.35 suggested that the reference chemical potential difference, Δμ0w, varies with the size of the guest molecule instead of using a single value for all the guest molecules, as there is a distortion in the lattice when the size of the guest molecule is increased.36 Pradhan found that the reference chemical potential difference value increases with the increase in size of the guest molecule by fitting the experimental data while slightly adjusting the Kihara parameters for some guest molecules.37 The molecular diameter of CO2 molecule is 5.12 Å and that for CH4 is 4.36 Å. The large size of the carbon dioxide molecule compared to methane might cause lattice distortion.

where Δμβ−L,α (T0, 0) is the reference chemical potential difw ference at the reference temperature, T0, and zero pressure. The reference temperature, T0, is the ice point temperature. In the case of methane hydrate, the ice point temperature T0 is 273.15 K, and in the case of carbon dioxide hydrate, T0 is 271.75 K. The depression in the ice point temperature for CO2 hydrate is due to the solubility of carbon dioxide in water. The last term in eq 9 is the activity of water, aw, defined as aw =

(12)

where is the reference heat capacity difference at the reference temperature T0. The constant, b, represents the dependence of heat capacity on the temperature. Two different expressions must be used for the water in liquid phase and in solid phase. The volume difference, ΔVβ−L,α , is assumed to be w constant. Reference Parameters. The thermodynamic reference properties used in this work are given in Table 1. Many investigators

(8)

−

(11)

ΔC0p(T0)

Φgg is the total interaction energy between the carbon dioxide molecule in cavity j and with the rest of the carbon dioxide molecules. The guest−guest interaction energy depends on the extent of surrounding cage occupancy

Δμwβ − L, α (T0 , 0)

ΔCp dT

where ΔH0w(T0) is the reference enthalpy difference between the empty hydrate lattice and the pure water phase at reference temperature T0. The heat capacity difference between the empty hydrate lattice and the pure water phase, ΔCp, is also temperature dependent, and it is approximated by the following expression

∫

=

T

0

r 2 sin θ dr dθ dϕ sin β dβ dγ ⎛ Φgh(r , θ , ϕ , β , γ ) ⎞ 1 ⎟ exp⎜ − = C jigg 4πkT V kT ⎠ ⎝

Δμwβ − L, α (T , P)

∫T

(10)

f Lw

where is the fugacity of water in the water rich aqueous phase and f w is the water fugacity at the reference state, the pure water phase. The temperature-dependent enthalpy difference is given by eq 11. 579

dx.doi.org/10.1021/jp410306v | J. Phys. Chem. B 2014, 118, 577−589

The Journal of Physical Chemistry B

Article

source of error in the quantum calculations is due to the electron correlation used and the error in using the incomplete basis set. The effect of the electron correlation energy was corrected by using the second-order Møller−Plesset (MP2) perturbation method.39 The CO2−H2O geometry was optimized for the CO2−H2O complex using the MP2 level of theory and the 6-31G basis set. The optimized structure obtained was a T-shaped structure with minimum energy distance between the carbon of carbon dioxide and oxygen of water molecule, r(C···O) = 2.788 Å, which is in agreement with the experimentally determined r(C···O) = 2.836 Å by Peterson and Klemperer.40 The calculated r value is somewhat less, and this might be because of the smaller basis set used for the geometry optimization. Optimum Method for PES Calculation. In order to determine the basis set that balances accuracy and computational cost, binding energy calculations were performed on the geometry optimized (minimum energy) H2O−CO2 complex using the MP2 level of electron correlation at increasing basis set. The calculated uncorrected binding energies and BSSE corrected binding energies by the full counterpoise method are plotted in

wave function. However, for any system except the smallest system, exact solutions to the Schrödinger equation are not computationally practical. Therefore, a great number of approximate methods strive to achieve the best trade-off between accuracy and computational cost. The simplest type of the ab initio electronic structure calculation is the Hartree−Fock (HF) scheme, in which the instantaneous Coulombic electron−electron repulsion is not specifically taken into account; only its average effect is included in the calculations. The energy obtained with this inaccurate approximation is always equal to or greater than the exact energy and tends to a limiting value called the Hartree−Fock limit as the basis set size increases. A basis set is a mathematical representation of the molecular orbital within a molecule. The basis set can be interpreted as restricting each electron to a particular region of space through the use of probability functions. The use of larger basis sets includes more probability density functions and thus imposes fewer constraints on electrons, allowing more flexibility to occupy orbitals and more accurately approximate exact molecular orbitals. The basis set for computing the potential energy hypersurface was carefully selected, considering accuracy and the computational cost. The interaction energy(ΔE) is the difference in energies between the dimer (H2O−CO2) and the monomers (CO2, H2O). ΔE BE = E H2O − CO2 − E H2O − ECO2

(14)

The calculation of interaction energies is subject to basis set superposition error (BSSE) when a finite basis set is used.38 When calculating the interaction energy, the atoms of interacting molecules approach one another and their basis functions overlap, resulting in more basis functions for the complex than the individual molecules. The difference in energy obtained from using different basis sets in each individual calculation results in BSSE. The counterpoise (CP) approach calculates the BSSE by reperforming all the calculations using the mixed basis sets, through introducing “ghost” atoms. A ghost atom has no protons or electrons but has the same number of molecular orbitals as the original atom, so it has the basis functions which can be used by another atom to increase the available basis functions. The full counterpoise corrected binding energy between the CO2 and H2O molecules can be expressed as ΔE BE,CP = E H2O − CO2 − E H2O,CP − ECO2,CP

Figure 1. Effect of increasing basis set size on the BSSE.

Figure 1. From Figure 1, it can be observed that the binding energy calculated with full counterpoise and without counterpoise (uncorrected) exhibit opposite trends as the basis set is increased from cc-pVDZ to aug-cc-pV5Z. The binding energy computed decreases as the number of basis functions increases for the uncorrected energy, making the binding weaker, while, for the BSSE corrected with the full counterpoise method, the binding energy increases with the number of basis functions, making the binding stronger. As the basis set size is increased from cc-pVDZ to aug-cc-pV5Z, the deviations of binding energies calculated with the full counterpoise method and without the counterpoise method decrease and will finally merge if a still higher basis set is used. We can see from Figure 1 that the average value (half counterpoise method) of the corrected and uncorrected binding energies gives the exact binding energy of the dimer. The binding energies obtained using the half counterpoise method using eq 16 are shown in Figure 2 along with the computational time required to calculate the binding energies. From Figure 2, one can determine the optimum basis set, with trade-offs between the calculation time and accuracy, to be the aug-cc-pVTZ basis set. The binding energies obtained using the half counterpoise method by aug-cc-pVTZ fall within the energy calculated by aug-cc-pV5Z. This aug-cc-pVTZ basis set therefore was used for the calculation of the potential energy surface. Calculation of Potential Energy Surface. The orientations between the water molecule and the carbon dioxide

(15)

EH2O−CO2 is the total energy of the CO2 and H2O dimer calculated at any particular configuration. EH2O,CP is the counterpoise corrected energy for a water molecule with CO2 as a ghost atom, and ECO2,CP is the counterpoise corrected energy for a carbon dioxide molecule with H2O as a ghost atom. The half counterpoise method is the average of uncorrected energy (ΔEBE) and corrected energy by the full counterpoise method (ΔEBE,CP) and is given in eq 16. Gaussian 03 was used for all quantum mechanical calculations. ΔE BE,halfCP =

1 (ΔE BE + ΔE BE,CP) 2

(16)

Determination of Pair Potentials. Guest−host pair potentials for CO2 can be developed by fitting atomic site−site potentials calculated using analytical potential models to the ab initio potential energy surface (PES). The optimum method and basis set must be determined for the accurate and efficient calculation of the CO2 and H2O interaction energies, which can then be used for the calculation of the potential energy surface. The general 580

dx.doi.org/10.1021/jp410306v | J. Phys. Chem. B 2014, 118, 577−589

The Journal of Physical Chemistry B

Article

Figure 2. Calculation time and binding energy (BSSE corrected by the half counterpoise method) at each basis set for the CO2−H2O complex.

molecules were varied, similar to that of Cao et al.14 for CH4− H2O, based on the actual clathrate structure. In this, the position of the water molecule was fixed and the orientation of the carbon dioxide molecule was varied in five dimensions for the interaction energy calculations. By inspecting the ball-and-stick model of the structure I clathrate hydrate, the relative orientations between guest and host water molecule fall into two types characterizing the plane containing the water molecule, as shown in Figure 3.

Figure 4. Five-dimensional orientation of carbon dioxide and water complex.

2.4−6.0 Å with 10 equally separated points sampled in the radial direction. (2) The ranges of the polar angle ξ and azimuthal angle φ were determined by moving a guest molecule some minimum distance inside a cage not too close to the cage wall. The ξ and φ angles were selected to range from −40 to 40° with five equally spaced angular points sampled for the carbon dioxide−water space sampling. (3) The rigid body of the carbon dioxide guest molecule was rotated in its own internal coordinates described by Euler angles β and γ. The angle β was spaced between 0 and 360° with four points sampled at 0, 90, 180, and 270°. The γ was spaced between 0 and 120° with three points sampled at 0, 40, and 80°. For each radial position r, there are two different water planes and 5 × 5 × 4 × 3 = 300 angular orientations of the carbon dioxide molecule. Anderson et al.27 developed a site−site method, which can be used to obtain the potentials from the 300 × 2 = 600 interaction energies at each radial position. These potential functions can be incorporated into the configurational integral to evaluate the Langmuir constant. Prior to the work of Anderson et al.,27 Cao et al.14 performed the angle average using the Boltzmann-weighted average over the five angular orientations (ξ, φ, α, β, γ) at each radial point r for methane. The work of Cao et al. was corrected by Anderson et al. by employing the site−site method instead of averaging all five orientational and rotational degrees of freedom. This resulted in better prediction of phase equilibria and cage occupancies. Therefore, a site−site model developed by Anderson et al.27 was used for the CO2−H2O interactions to account for all degrees of freedom. The radial position, r, is equally spaced in 10 different points comprised of 10 × 600 = 6000 configurations between the CO2 and H2O.

Figure 3. Planar orientation of a water molecule: (a) water plane parallel to the page, plane-1; (b) water plane perpendicular to the page, plane-2.

The different orientations are generated by fixing the plane of the water molecule in one orientation and moving the carbon dioxide guest molecule in a five-dimensional grid to different positions inside the cage. The center of mass of the carbon dioxide molecule was moved in a polar coordinate system (r, ξ, φ) with the oxygen of the water molecule as the origin. This coordinate system defines the position of carbon dioxide carbon with the water oxygen as the origin. Furthermore, the rigid body of the carbon dioxide guest molecule was rotated in its own internal coordinates, Euler angles α, β, γ. The five-dimensional grid between the CO2 and water molecules is illustrated in Figure 4. The Euler angle, β, is the degrees of angle rotated in counter clockwise direction along the ỳ-axis, and γ is the degrees of angle rotated in counter clockwise direction along the z̀-axis. CO2 is symmetric along the x̀-axis so there will not be any change in energy by varying the Euler angle, α. Thus, only Euler angles β and γ were considered in the generation of a grid. The limits of the r, ξ, and φ dimensions and the Euler angles are determined in the following manner: (1) The r distance between the center of the carbon of the CO2 molecule and the oxygen of the water molecule cannot be greater than the maximum diameter of the cage, because the guest molecule is entrapped in a cage. The interaction energy will become extremely repulsive when the separation distance is very small. Considering the above limitation, the separation distance, r, was set at

■

RESULTS AND DISCUSSION Potential Fit to Intermolecular Energies. The interaction energies calculated at 6000 different CO2−H2O configurations were fitted to site−site potentials based on the interactions between the center of the guest (carbon of CO2) molecule and the oxygen on the water (O2C−OH2), the oxygen of carbon dioxide with oxygen on the water (OCO−OH2), and the carbon 581

dx.doi.org/10.1021/jp410306v | J. Phys. Chem. B 2014, 118, 577−589

The Journal of Physical Chemistry B

Article

where 1/4πε0 is the proportionality constant with ε0 as the electric constant. The proportionality constant value is 8.98755178 × 109 N-m2/C2. qi and qj are the charges on the ith and jth atoms. A nonlinear least-squares fitting method was used to fit the ab initio potential data to the potential models. The O2C−OH2 interaction was modeled using the exponential-6 potential form, the OCO−OH2 interaction was modeled using the LennardJones 6-12 potential form, and the OCO−M, O2C−HOH, O2C−M, and OCO−HOH interactions were modeled using the Coulombic charge−charge form. The geometry assumed by the TIP4P/ice model42 for H2O was used in fitting the potential. Studies by Vega42−44 have shown that the TIP4P/ice model gives the most accurate predictions for the melting point of ice and methane hydrate dissociation pressure compared to other water models. The TIP4P/ice model as an additional interacting site “M” was located on the bisector of the H−O−H angle 0.1577 Å away from the oxygen atom toward the hydrogen atom. The TIP4P/ice42 water model is given in Figure 14, where q1 = 0.5897 is the charge on the hydrogen atoms and the charge on M is −1.1794. The 6000 ab initio energies were fitted to site−site potentials minimizing the Boltzmann factor-weighted objective function χ given in eq 20, thereby optimizing the parameters of L-J 6-12 and exponential-6. The parameter k is the Boltzmann constant, and the temperature T is 273.15 K. The temperature, T, is taken as 273.15 K because, in nature, hydrates generally occur in and around the temperature, T, of 273.15 K. Also, there is no effect of variation of temperature in the minimization of the objective function. The charges were held constant and are given in Table 2.

on the CO2 with the hydrogen of the water (O2C−HOH). The atoms marked with bold indicate the location of the interaction site for use in the site−site potential. The O2C−OH2 potential captures the guest position effects, the O2C−HOH potential captures the ξ and φ orientations in addition to the radial dimension, r, and the OCO−OH2 potential captures the Euler orientation of the carbon dioxide molecule with respect to the water molecule. In order to implement the ab initio interaction potentials in the configurational integral for the thermodynamics calculations, the calculated potential must be accurately represented in a mathematical form. An intermolecular potential model was used to represent the entire guest−host interaction potential accurately. The Lennard-Jones 6-1241 and exponential641 potential models along with the Coulombic charge−charge formula were used to fit the ab initio potential. The potential forms are as follows: Lennard-Jones 6-12: ⎡⎛ ⎞12 ⎛ ⎞6 ⎤ σ σ ΦLJ(rij) = 4εij⎢⎢⎜⎜ ⎟⎟ − ⎜⎜ ⎟⎟ ⎥⎥ r r ⎝ ij ⎠ ⎦ ⎣⎝ ij ⎠

(17)

Exponential-6: 6⎤ ⎡ ⎡ ⎛ rij ⎞⎤ ⎛ rm ⎞ ⎥ 6 ⎢ ⎜ ⎟ ⎢ ⎥ Φexp ‐6(rij) = exp γ ⎜1 − ⎟ − ⎜ ⎟ ⎥ ⎢⎣ ⎝ 1 − (6/γ ) ⎢ γ rm ⎠⎥⎦ ⎝ rij ⎠ ⎣ ⎦

εij

(18) 41

The Coulombic charge−charge formula is given as 1 qiqj Φcharge(rij) = 4πε0 rij

(19)

# data points

χ=

Table 2. CO2−H2O Potential Parameters by the VAS Model (Atoms Marked in Bold Indicate Interaction Sites) Exp-6 O2C−OH2 OCO−OH2 CO2 OCO H2O HOH

j

(20)

L-J 6-12

ε/k (K)

rm (Å)

γ

81.58

3.775

10.698

ε/k (K)

σ (Å)

77.95

3.0808

∑

2 ⎡ ⎛ ⎞ ⎛ −ΔEcp, j ⎞ ⎤ ⎢exp⎜ −ΔEcp, j ⎟ ⎥ − exp⎜ ⎟ ⎢ ⎝ kT ⎠ ⎝ kT ⎠QM ⎥⎦ ⎣ pred

The results of the prediction of the site−site binding energies and the ab initio binding energies are shown by a diagonal plot in Figure 5 for water plane 1 and in Figure 6 for water plane 2. The diagonal plot is a parity plot matrix showing the number of points of binding energies that are in the specified range for energies calculated by the ab initio method and by the site−site potential method. The x-axis on the diagonal plot is the site−site binding energies, and the y-axis is the ab initio binding energies.

charge

0.652 −0.326 0.000 0.5897

Figure 5. Parity plot for water plane-1 showing the number of binding energy points. 582

dx.doi.org/10.1021/jp410306v | J. Phys. Chem. B 2014, 118, 577−589

The Journal of Physical Chemistry B

Article

Figure 6. Parity plot for water plane-2 showing the number of binding energy points.

The binding energies ranging from −2 to 1.25 kcal/mol have been considered in the diagonal plot, because the attractive region is important as the energies are Boltzmann weighted when optimizing the potential parameters; see eq 20. The diagonal in Figures 5 and 6, highlighted in color, represents the number of points that have the same range of binding energies calculated by the ab initio method and by site−site potentials. Many-Body Effects. Klauda and Sandler24 showed that many-body effects can significantly change the total interaction energy between the guest molecule and the clathrate cage. They also showed for the methane hydrate that the two half cell calculations closely resemble the calculations of a full cage. To evaluate the many-body effects in the carbon dioxide hydrate system, initially, the half pentagonal dodecahedron of structure I with more than half water molecules, 15 water molecules, with a single guest carbon dioxide molecule is optimized for the minimum energy at the MP2/6-31G level. The 15 water molecules and guest carbon dioxide system are shown in Figure 7. The guest molecule inside the half cage is moved in different configurations, and the interaction energy was calculated for these 15 water

Figure 8. Parity plot of corrected site−site predicted 15 water molecule−carbon dioxide interaction energies.

molecules and single guest CO2 molecule. Six different configurations were obtained by moving the guest CO2 molecule toward the cage and also by rotating the CO2 molecule w.r.t. the 15 water molecule cell. Preliminary calculations were carried out at the MP2/aug-cc-pVTZ basis level, similar to the basis set used for PES calculations, but the computational time required for the interaction energy calculation for the 16 molecule system is more than a month with the available resources. Due to the computational limitations, the interaction energies were calculated at the MP2/6-31++G (2d, 2p) level for different configurations of guest in the 16 water molecule cell. The computational time required at the MP2/6-31++G (2d, 2p) level basis set is around 12 h. The site−site model was used to calculate the total interaction energy of the many-body system. The atomic site−site potentials obtained by optimizing the 6000 point ab initio potential energy surface reproduce the many-body interaction energies very well. There was no change in the site−site potential parameters when we tried to optimized them, such that the errors of the prediction of the site−site model w.r.t. the ab initio half cell calculations were minimized using the Boltzmann factor-weighted objective function χ given in eq 20. Figure 8 shows the results of the binding energies calculated on the 15 water molecules−CO2 system. We denote the site−site ab initio intermolecular potential model as the “VAS” model. The VAS model parameters for CO2−H2O are listed in Table 2. The pure CO2 intermolecular potentials were obtained from pair CO2−H2O potentials given in Table 2 by using the Lorentz−Berthelot41 combining rules. The force field used for water was the TIP4P/ice model.42 The pure

Figure 7. Single guest CO2 and 15 water molecules of the pentagonal dodecahedron of the structure I hydrate. Red color, oxygen atoms; white color, hydrogen atoms; black color, carbon atom. The dotted lines represent hydrogen bonding. 583

dx.doi.org/10.1021/jp410306v | J. Phys. Chem. B 2014, 118, 577−589

The Journal of Physical Chemistry B

Article

Table 3. CO2 Intermolecular Potential Parametera VAS model

a

EPM2 model46

L-J 6-12

ε/k (K)

σ (Å)

ε/k (K)

σ (Å)

O2C−CO2 OCO−OCO CO2 CO2

73.351 57.677

3.493 2.994

28.129 80.507

2.757 3.033

charge

0.652 −0.326

Bold atoms indicate specific atom−atom interactions.

carbon dioxide intermolecular potential parameters are given in Table 3. Reference Parameters. Holder et al.45 rearranged and integrated eqs 1, 2, and 9−12 and gave the following equation Δμwβ − H (T ,

P)

+

ΔCpβ − L, α(T0)

Figure 9. Thermodynamic reference parameters for structure I CO2 hydrate. The source of experimental data: Miller and Smythe,58 Falabella,59 Larson,60 Robinson and Mehta,61 Deaton and Frost,62 Ng and Robinson,63 Unruh and Katz,64 Adisasmito et al.,65 and Ohgaki et al.66

− bT0 ln(T /T0)

RT R β − L, α β − L, α 2 2ΔCp (T0) − b T0 ⎡ 1 1⎤ − ⎢ − ⎥ 2R T⎦ ⎣ T0

To account for possible lattice constant distortion when hydrates are formed from CO2, we have calculated the reference parameters using the EPM2 model by allowing a variable lattice constant of the hydrate. These calculations were done to know whether the EPM2 model could improve the hydrate properties if the flexibility of the water lattice of the hydrate phase48 is incorporated into the van der Waals and Platteeuw model. The increase in lattice constant will create an expanded water lattice simulating flexibility in the crystal structure to accommodate the CO2 molecule. The lattice constant was varied in the range 11.95−12.7 Å, and the reference chemical potential and reference enthalpy difference calculated are summarized in Table 4. There is certainly an increase in the reference chemical

β − L, α

P ΔVw b β − L, α (T − T0) − + ln a w 2R RT Δμw0 ΔH wβ − L, α(T0) ⎡ 1 1⎤ = + ⎢ − ⎥ RT R T0 ⎦ ⎣T +

(21)

If the left-hand side of eq 21 is defined as the variable Y, then Y=

Δμw0 RT

+

ΔH w0 + ΔH wL − α(T0) ⎛ 1 1⎞ ⎜ − ⎟ R T0 ⎠ ⎝T

(22)

ΔHL−α w (T0)

where is the latent heat of ice, Y is a function of the experimental conditions (T, P, composition), and other constants are given in Table 1. Using a linear regression analysis, the reference thermodynamic parameters obtained are Δμ0w = 1206 ± 3 J/mol and ΔH0w = 1260 ± 12 J/mol. The plot for the reference parameters is shown in Figure 9. The reference parameters were also calculated by the above method using Harris and Yung46 intermolecular potentials (EPM2 model) and using Potoff and Siepmann47 carbon dioxide and water intermolecular potentials. The intermolecular potentials for the carbon dioxide and water system are calculated using the combining rules, i.e., the Lorentz−Berthelot41 combining rules given in eqs 23 and 24, and the potentials for water are from the TIP4P-ice model. σii + σjj σij = (23) 2 εij =

εiiεii

Table 4. Reference Chemical Potential Difference and Reference Enthalpy Difference Calculated at Different Lattice Constants of Hydrate Using the EPM2 Model46 lattice constant (Å)

Δμ0w (J/mol)

ΔH0w (J/mol)

11.95 12.00 12.03 12.07 12.20 12.30 12.50 12.70

304.50 ± 3.6 318.56 ± 3.5 330.71 ± 3.2 336.64 ± 5.0 366.41 ± 2.6 387.92 ± 2.6 421.81 ± 2.6 447.04 ± 2.4

515.46 ± 21 482.20 ± 20 440.60 ± 18 387.15 ± 22 346.48 ± 15 273.13 ± 13 149.65 ± 11 45.74 ± 5

potential value, but it still falls well outside the values reported by earlier researchers for empty hydrate lattice reference parameters. Further, the reference enthalpy difference decreases with an increase in the lattice constant. The EMP2 model results using a flexible lattice still result in reference parameters that lie far outside those determined in previous studies. Phase Equilibrium. By solving eqs 1−12 using the reference parameters obtained by this model, the dissociation pressure is calculated at the desired temperature. Figures 10 and 11 compare the equilibrium pressure of CO2 hydrate at various temperatures ranging from 155 to 283.3 K with the experimental data. The absolute average deviation is less than 3.2% from the experimental data. Cage Occupancies. Cage occupancies, the fraction of each cage occupied by a guest molecule, are important, as they tell the amount of gas stored in the hydrate or the amount of gas that can be stored in the hydrate. The hydration number, n, can be

(24)

The reference parameters obtained by using the EPM2 model are Δμ0w = 331 ± 3 J/mol and ΔH0w = 440 ± 18 J/mol, and those for Potoff and Siepmann potentials are Δμ0w = 392.4 ± 3 J/mol and ΔH0w = 490 ± 12 J/mol. The reference parameters were obtained well outside the range obtained by earlier researchers either numerically or experimentally, as given in Table 1 of Anderson et al.27 for structure I hydrate. This shows the inability of the Harris and Yung potentials and Potoff and Siepmann potentials to accurately model carbon dioxide hydrates using the van der Waals and Platteeuw5 model framework. This also would call into question its applicability for molecular dynamic simulations of gas hydrates. 584

dx.doi.org/10.1021/jp410306v | J. Phys. Chem. B 2014, 118, 577−589

The Journal of Physical Chemistry B

Article

determined from the cage occupancies, as the hydration number is the average number of water molecules per guest molecule in the hydrate. For structure I hydrate, the hydration number can be calculated using eq 25. For fully occupied large cages, θL = 1, and small cages, θS = 1, of structure I gas hydrate, the hydration number is 5.75. 23 n= 3θL + θS (25)

temperatures between 263 and 278 K with pressures well above the equilibrium pressures, around 60 atm. The cage occupancies reported by Udachin et al.21 from the single crystal X-ray diffraction studies were 100% for the large cage (θL) and 71% for the small cage (θS); this yields the hydration number of 6.20. They prepared the crystal at a temperature of 276 K in the presence of excess liquid CO2 and a pressure of almost twice that of the equilibrium condition, at 38 atm. In all of the above experiments, the small cage occupancy reported for carbon dioxide hydrate is in the range 60−80% and the large cage occupancy is more than 98%. The CO2 hydrate samples prepared for determining the cage occupancies and hydration numbers in the above experiments were well above the equilibrium pressures, and these higher pressures during the synthesis can produce higher occupancies. Ripmeester and Ractliff23 prepared a sample under equilibrium conditions, at a temperature of 268 K and a pressure of 9.9 bar, and gave a lower limit to the hydration number of 7.0 for CO2 hydrate. They used solid state NMR to measure the relative cage occupancy, θS/θL, of 0.32 and assumed a Δμ0w value of 1297 J/m for a hydration number calculation, so the small cage occupancy was nearly 0.3136 assuming the 98% occupancy for a large cage. This means that the hydration number should be higher than 7.0 and the small cage occupancy should be in the range 25−40%. Cage occupancies were calculated at a particular temperature from eq 4 using the Langmuir constant obtained from the VAS model. CSMGEM is a thermodynamic code developed by Sloan and co-workers1 (Colorado School of Mines) to predict the phase equilibrium of the hydrate, and it uses the fitted Kihara potential parameters in predicting the occupancies and phase equilibria.1 The cage occupancy predicted by CSMGEM for a small cage is in between 47 and 40% in the temperature range between 256 and 283.3 K and almost fully occupied for large cages, 97% occupancy for large cage. CSMGEM predicted the phase equilibrium of carbon dioxide hydrate accurately, but it overestimates the cage occupancies. Klauda and Sandler24 predicted the small cage occupancy in between 54 and 90% in the temperature between 243.1 and 290 K. Sun and Duan25 using the site−site ab initio model had reported the hydration number for only two temperatures, at equilibrium conditions, at 273.1 and 274.5 K. The small cage occupancy was calculated for Sun and Duan data from hydration number, assuming 99% occupancy for large cage and obtained 55 and 60% occupancy at 273.15 and 274.5 K. In Figure 12, the predictions for the cage occupancy ratios (θS/θL) for the carbon dioxide hydrates obtained by using the

Spectroscopic measurements such as NMR and Raman have been used by different researchers to calculate the cage occupancy in which the integrated signal intensity ratios of the guests in the two hydrate cavities are measured.49 The signal intensity ratios between peaks for guests in each cage type reproduce the ratios of the cage occupancies (θS/θL, small cage to large cage) of the guest in the lattice cages. Sum et al.49 using Raman spectroscopy reported that the guest carbon dioxide appears to occupy only the large cavities of the sI hydrate, indicating 0% occupancy in the small cage. Later, NMR studies23 and neutron diffraction studies20 of CO2 hydrates show that the carbon dioxide molecules occupy both types of cavities in the clathrates. The cage occupancies determined by Henning et al.20 from neutron diffraction studies for the CO2 guest were more than 95% for the large cavities (51262), and for the small cages (512), they were in the range of 60−80%. This gives hydration numbers between 6.05 and 6.67. They prepared the sample at

Figure 12. Cage occupancy of CO2 hydrate at temperatures ranging from 155 to 283 K: (black dotted line) Klauda and Sandler;24 (yellow dashed line) CSMGEM;1 (blue triangles) Sun and Duan;25 (black square) Ripmeester and Ratcliffe;23 (red solid line) VAS model.

Figure 10. Calculation of CO2 hydrate equilibrium dissociation pressure using ab initio site−site potentials and regressed reference parameters for CO2. The source of experimental data: Miller and Smythe,58 Falabella,59 Larson,60 Robinson and Mehta,61 Deaton and Frost,62 Ng and Robinson,63 Unruh and Katz,64 Adisasmito et al.,65 and Ohgaki et al.66

Figure 11. Calculation of CO2 hydrate equilibrium dissociation pressure for T > 260 K, using ab initio site−site potentials and regressed reference parameters for CO2. The source of experimental data: Miller and Smythe,58 Falabella,59 Larson,60 Robinson and Mehta,61 Deaton and Frost,62 Ng and Robinson,63 Unruh and Katz,64 Adisasmito et al.,65 and Ohgaki et al.66

585

dx.doi.org/10.1021/jp410306v | J. Phys. Chem. B 2014, 118, 577−589

The Journal of Physical Chemistry B

Article

VAS model and by other researchers are compared. The small cage occupancies predicted by the VAS model for carbon dioxide structure I hydrate are in the range 32−51% for temperatures ranging from 155.5 to 283.3 K, whereas the large cage is more than 98% occupied. The cage occupancies predicted by Sloan/ CSMGEM, Klauda and Sandler, and Sun and Duan overestimated the small cage occupancies. Figure 13 compares the

Figure 14. TIP4P-ice water model.42

Figure 13. Hydration number for CO2 hydrate at different temperatures: (black dotted line) Klauda and Sandler;24 (yellow dashed line) CSMGEM;1 (blue triangles) Sun and Duan;25 (black square) Ripmeester and Ratcliffe;23 (red solid line) VAS model.

hydration number predicted by this model and by other researchers.1,23−25

■

DENSITY AND DIFFUSION CALCULATIONS In order to test the VAS model, molecular dynamics (MD) simulations were carried out calculating the density of pure CO2 and diffusivity of CO2 in water. The force field used for CO2 in the simulations is given in Table 3. All simulations were performed in a constant isothermal−isobaric (NPT) ensemble. The temperature was controlled by a Nosé−Hoover thermostat23,50,51 and the pressure by a Parrinello−Rahman52 barostat. Coulombic long-range interactions were calculated using the Ewald summation method, and non-electrostatic L−J interactions were calculated within a cutoff distance of 9 Å. SHAKE53 was used to maintain the geometry of CO2 and water molecules with a relative tolerance of 10−4 nm. Periodic boundary conditions were used in all the directions for all the simulations. The GROMACS32 package was used to perform the MD simulations. The simulation box consists of a cubic box with 1000 pure CO2 (4.4 Å × 4.4 Å × 4.4 Å) molecules for the density calculations, as given in Figure 15. The diffusivity calculations were performed on a system that contained a total of 510 molecules. The first step of each MD run consists of energy minimization with steepestdescent; this is to reduce the thermal noise in the structure and potential energies. The second step consists of a production run. MD simulations were run at different temperatures and pressures for density and at different temperatures and 1 bar pressure for diffusivity calculations. The average density was calculated after the equilibrium was reached: the calculated average temperature and pressure were the same as the imposed values, and no drift in the conserved quantity was observed. The densities simulated are given in Figure 16 using the VAS model. The densities are compared to the data obtained from fluid properties from the Span and Wagner equation of state (EOS) model.54 The densities obtained from the VAS model are in good agreement with data using the Span and Wagner EOS. In all the density calculation simulations, the initial configuration was in

Figure 15. Simulation box of pure CO2. Ball-and-stick model of CO2 with carbon and oxygen atoms shown in cyan and red, respectively.

Figure 16. Densities of pure CO2. Solid lines, Span and Wagner EOS;54 points, the VAS model. The vertical dashed line indicates the phase change from liquid to gas.

the condensed phase. For the simulation at 300 K and 50 bar, the density obtained was in the liquid phase, whereas it was supposed to be in the vapor phase. This is due to the metastable states and barriers for nucleation, so the vapor phase density was simulated from an initial configuration which is in the vapor phase. The density given in Figure 16 at 300 K and 50 bar is from the starting configuration which is in the vapor phase. The VAS model is not able to predict the phase change. MD simulations were also 586

dx.doi.org/10.1021/jp410306v | J. Phys. Chem. B 2014, 118, 577−589

The Journal of Physical Chemistry B

Article

water, which predicts the lesser diffusivity of CO2 in water compared with the TIP4P model. Each intermolecular potential model has a range of applicability. The VAS model well predicts the liquid CO2 densities, the CO2 hydrate cage occupancies, and phase equilibria and diffusivity of CO2 in water but fails when used for the prediction of vapor phase densities. The VAS model can be suitable to use in molecular simulations for liquid CO2 and CO2 hydrates.

■

CONCLUSIONS A site−site method has been applied to develop the CO2−H2O intermolecular potentials that characterizes the five-dimensional potential energy surface. The ab initio intermolecular potentials obtained from the 6000 point hyperspace energy surface were corrected for many-body effects. The corrections were employed by fitting the intermolecular potentials to quantum mechanical calculations on a system with 15 water molecules interacting with one carbon dioxide molecule. The reference thermodynamic parameters were calculated for structure I carbon dioxide hydrate using the VAS model as Δμ0w = 1206 ± 2 J/mol and ΔH0w = 1260 ± 12 J/mol. The EPM2 model and Potoff and Siepmenn model for CO2 failed to predict the cage occupancies and phase equilibrium when applied to hydrates using the vdWP framework. The reference parameters obtained by using the EPM2 model and Potoff and Siepmenn model are well outside the range obtained by earlier researchers using numerical and experimental methods. Using the VAS model, intermolecular potentials, and the reference parameters calculated from our potentials, the phase equilibrium pressure was computed with less than 3% absolute average deviation from the experimental data. The small cage occupancy predicted by this model for structure I CO2 hydrate is in the range 32−51% for temperatures ranging from 155.5 to 283.3 K, and the large cage is more than 98.5% occupied in the temperature range. The intermolecular potentials were also tested by calculating densities of pure CO2 and diffusion coefficients of CO2 in water using MD simulations. The densities and diffusivities of CO2 in water simulated are in good agreement with the literature data for pure liquid CO2 but overestimate the density of gaseous pure CO2. There is no perfect model which predicts all the properties for CO2. The EPM2 model is good for VLE predictions, but it inaccurately predicts the condensed phases, as seen in liquid density calculations and for CO2 hydrate reference property calculations. The VAS model is able to predict the liquid CO2 densities and CO2 hydrate cage occupancy and phase equilibria well, but it fails to predict them for the vapor phase. A more complex model is needed for CO2 which can predict all the phase properties.

Figure 17. Densities of pure CO2. Solid lines, Span and Wagner EOS;54 points, the EPM246 model. The vertical dashed line indicates the phase change from liquid to gas.

performed to calculate the densities using the EPM2 model and are given in Figure 17. All the density simulations with the EPM2 model were performed from an initial configuration in the liquid phase. It is known that the EPM2 model predicts the VLE data accurately, as the parameters of the EMP2 model were obtained by optimizing to the VLE data. From Figure 17, we can see that the EPM2 model underpredicts the liquid densities at higher pressures compared with the values obtained by using the Span and Wagner EOS. The diffusion coefficients for infinitely diluted species I in J were calculated by the Einstein equation DI = lim t →∞

⟨(ri(t ) − ri(0))2 ⟩ 6t

(26)

where ri(t) is the position of particle i at time t. The simulations were performed with two water models, TIP4P55 and TIP4Pice,42 at different temperatures and 1 bar pressure. The diffusion coefficient of CO2 in water is more dependent on temperature than on pressure.56,57 The simulated CO2 diffusion coefficients in water are listed in Figure 18. The simulated CO2 diffusion

■

AUTHOR INFORMATION

Corresponding Author

Figure 18. Temperature dependence of the diffusion coefficient of dissolved CO2 in water. Points: using the VAS model. Solid line: experimental data56 (diffusivity derived from the experimental data).

*Phone: +1 304 293 9334. Fax: +1 304 293 4139. E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■ ■

coefficients in water are in good agreement with the experimental values56 in the range of temperatures investigated for the TIP4P water model, but when used with the TIP4P-ice model for water, the simulated diffusion coefficients underpredict the diffusion coefficients compared with the experimental values. The TIP4Pice model was specifically designed to cope with the ice properties, and it predicts the melting point of ice as 272.2 K.42 The strong interaction of the water molecules in the TIP4P-ice model makes it difficult for the CO2 molecules to diffuse in the liquid

ACKNOWLEDGMENTS This technical effort was performed in support of the National Energy Technology Laboratory’s research in methane hydrates. REFERENCES

(1) Sloan, E. D.; Koh, C. A. Clathrate Hydrates of Natural Gases; CRC Press: Boca Raton, FL, 2007.

587

dx.doi.org/10.1021/jp410306v | J. Phys. Chem. B 2014, 118, 577−589

The Journal of Physical Chemistry B

Article

(2) Ohgaki, K.; Inoue, Y. A Proposal for Gas Storage on the Bottom of the Ocean, Using Gas Hydrates. Int. Chem. Eng. 1994, 34, 417−419. (3) Ohgaki, K.; Takano, K.; Sangawa, H.; Matsubara, T.; Nakano, S. Methane Exploitation by Carbon Dioxide from Gas Hydrates-Phase Equilibria for CO2-CH4 Mixed Hydrate System. J. Chem. Eng. Jpn. 1996, 29, 478−483. (4) Ota, M.; Abe, Y.; Watanabe, M.; Smith, R. Methane Recovery from Methane Hydrate Using Pressurized CO2. Fluid Phase Equilib. 2005, 228, 553−559. (5) Van der Waals, J.; Platteeuw, J. Clathrate Solutions. Adv. Chem. Phys 1959, 2, 1−57. (6) Lennard-Jones, J.; Devonshire, A. Critical Phenomena in Gases. II. Vapour Pressures and Boiling Points. Proc. R. Soc. London, Ser. A 1938, 165, 1−11. (7) Lennard-Jones, J.; Devonshire, A. Critical Phenomena in Gases. I. Proc. R. Soc. London, Ser. A 1937, 163, 53−70. (8) Parrish, W. R.; Prausnitz, J. M. Dissociation Pressures of Gas Hydrates Formed by Gas Mixtures. Ind. Eng. Chem. Process Des. Dev. 1972, 11, 26−35. (9) John, V. T.; Holder, G. Choice of Cell Size in the Cell Theory of Hydrate Phase Gas-Water Interactions. J. Phys. Chem. 1981, 85, 1811− 1814. (10) John, V. T.; Holder, G. Contribution of Second and Subsequent Water Shells to the Potential Energy of Guest-Host Interactions in Clathrate Hydrates. J. Phys. Chem. 1982, 86, 455−459. (11) John, V. T.; Holder, G. D. Langmuir Constants for Spherical and Linear Molecules in Clathrate Hydrates. Validity of the Cell Theory. J. Phys. Chem. 1985, 89, 3279−3285. (12) Sparks, K. A.; Tester, J. W. Intermolecular Potential Energy of Water Clathrates: The Inadequacy of the Nearest-Neighbor Approximation. J. Phys. Chem. 1992, 96, 11022−11029. (13) Sparks, K. A.; Jefferson, W.; Cao, Z.; Trout, B. L. Configurational Properties of Water Clathrates: Monte Carlo and Multidimensional Integration Versus the Lennard-Jones and Devonshire Approximation. J. Phys. Chem. B 1999, 103, 6300−6308. (14) Cao, Z.; Tester, J.; Sparks, K.; Trout, B. Molecular Computations Using Robust Hydrocarbon Water Potentials for Predicting Gas Hydrate Phase Equilibria. J. Phys. Chem. B 2001, 105, 10950−10960. (15) John, V. T.; Papadopoulos, K.; Holder, G. A Generalized Model for Predicting Equilibrium Conditions for Gas Hydrates. AlChE J. 1985, 31, 252−259. (16) Bazant, M.; Trout, B. A Method to Extract Potentials from the Temperature Dependence of Langmuir Constants for ClathrateHydrates. Physica A 2001, 300, 139−173. (17) Klauda, J.; Sandler, S. A Fugacity Model for Gas Hydrate Phase Equilibria. Ind. Eng. Chem. Res. 2000, 39, 3377−3386. (18) Pimpalgaonkar, H.; Veesam, S. K.; Punnathanam, S. N. Theory of Gas Hydrates: Effect of the Approximation of Rigid Water Lattice. J. Phys. Chem. B 2011, 115, 10018−10026. (19) Ravipati, S.; Punnathanam, S. N. Analysis of Parameter Values in the Van Der Waals and Platteeuw Theory for Methane Hydrates Using Monte Carlo Molecular Simulations. Ind. Eng. Chem. Res. 2012, 51, 9419−9426. (20) Henning, R.; Schultz, A.; Thieu, V.; Halpern, Y. Neutron Diffraction Studies of CO2 Clathrate Hydrate: Formation from Deuterated Ice. J. Phys. Chem. A 2000, 104, 5066−5071. (21) Udachin, K.; Ratcliffe, C.; Ripmeester, J. Structure, Composition, and Thermal Expansion of CO2 Hydrate from Single Crystal X-Ray Diffraction Measurements. J. Phys. Chem. B 2001, 105, 4200−4204. (22) Ikeda, T.; Yamamuro, O.; Matsuo, T.; Mori, K.; Torii, S.; Kamiyama, T.; Izumi, F.; Ikeda, S.; Mae, S. Neutron Diffraction Study of Carbon Dioxide Clathrate Hydrate. J. Phys. Chem. Solids 1999, 60, 1527−1529. (23) Ripmeester, J.; Ratcliffe, C. The Diverse Nature of Dodecahedral Cages in Clathrate Hydrates as Revealed by 129Xe and 13C Nmr Spectroscopy: CO2 as a Small-Cage Guest. Energy Fuels 1998, 12, 197− 200. (24) Klauda, J.; Sandler, S. Ab Initio Intermolecular Potentials for Gas Hydrates and Their Predictions. J. Phys. Chem. B 2002, 106, 5722−5732.

(25) Sun, R.; Duan, Z. Prediction of CH4 and CO2 Hydrate Phase Equilibrium and Cage Occupancy from Ab Initio Intermolecular Potentials. Geochim. Cosmochim. Acta 2005, 69, 4411−4424. (26) Cao, Z.; Tester, J.; Trout, B. Computation of the Methane-Water Potential Energy Hypersurface Via Ab Initio Methods. J. Chem. Phys. 2001, 115, 2550. (27) Anderson, B.; Tester, J.; Trout, B. Accurate Potentials for Argon Water and Methane Water Interactions Via Ab Initio Methods and Their Application to Clathrate Hydrates. J. Phys. Chem. B 2004, 108, 18705−18715. (28) Peng, D.; Robinson, D. A New Two-Constant Equation of State. Ind. Eng. Chem. Fundam. 1976, 15, 59−64. (29) McMullan, R. K.; Jeffrey, G. Polyhedral Clathrate Hydrates. Ix. Structure of Ethylene Oxide Hydrate. J. Chem. Phys. 1965, 42, 2725. (30) Mak, T. C.; McMullan, R. K. Polyhedral Clathrate Hydrates. X. Structure of the Double Hydrate of Tetrahydrofuran and Hydrogen Sulfide. J. Chem. Phys. 1965, 42, 2732−2737. (31) B. Klauda, J.; I. Sandler, S. Phase Behavior of Clathrate Hydrates: A Model for Single and Multiple Gas Component Hydrates. Chem. Eng. Sci. 2003, 58, 27−41. (32) van der Spoel, D.; Lindahl, E.; Hess, B.; Van Buuren, A.; Apol, E.; Meulenhoff, P.; Tieleman, D.; Sijbers, A.; Feenstra, K.; van Drunen, R. Gromacs User Manual Version 3.3; www.gromacs.org, 2008. (33) Carnahan, B.; Luther, H.; Wilkes, J. Applied Numerical Methods; Wiley: New York, 1969. (34) Holder, G.; Corbin, G.; Papadopoulos, K. Thermodynamic and Molecular Properties of Gas Hydrates from Mixtures Containing Methane, Argon, and Krypton. Ind. Eng. Chem. Fundam. 1980, 19, 282− 286. (35) Holder, G.; Zetts, S.; Pradhan, N. Phase Behavior in Systems Containing Clathrate Hydrates. Rev. Chem. Eng. 1988, 5, 1−70. (36) Zele, S.; Lee, S.; Holder, G. A Theory of Lattice Distortion in Gas Hydrates. J. Phys. Chem. B 1999, 103, 10250−10257. (37) Pradhan, N. Prediction of Multi-Phase Equilibria in Gas Hydrates. M.S. Thesis, University of Pittsburgh, Pittsburgh, PA, 1985. (38) Boys, S.; Bernardi, F. The Calculation of Small Molecular Interactions by the Differences of Separate Total Energies. Some Procedures with Reduced Errors. (Reprinted from Mol. Phys. 1970, 19, 553−566). Mol. Phys. 2002, 100, 65−73. (39) Møller, C.; Plesset, M. Note on an Approximation Treatment for Many-Electron Systems. Phys. Rev. 1934, 46, 618−622. (40) Peterson, K.; Klemperer, W. Structure and Internal Rotation of H2O-CO2, HDO-CO2, and D2O-CO2 Van Der Waals Complexes. J. Chem. Phys. 1984, 80, 2439−2445. (41) Tester, J.; Modell, M. Thermodynamics and Its Applications; Prentice Hall PTR: Upper Saddle River, NJ, 1997. (42) Abascal, J.; Sanz, E.; Fernández, R. G.; Vega, C. A Potential Model for the Study of Ices and Amorphous Water: TIP4P/Ice. J. Chem. Phys. 2005, 122, 234511−234520. (43) Conde, M.; Vega, C. Note: A Simple Correlation to Locate the Three Phase Coexistence Line in Methane-Hydrate Simulations. J. Chem. Phys. 2013, 138, 056101−056102. (44) Conde, M.; Vega, C. Determining the Three-Phase Coexistence Line in Methane Hydrates Using Computer Simulations. J. Chem. Phys. 2010, 133, 064507−064512. (45) Holder, G. D.; Malekar, S. T.; Sloan, E. D. Determination of Hydrate Thermodynamic Reference Properties from Experimental Hydrate Composition Data. Ind. Eng. Chem. Fundam. 1984, 23, 123− 126. (46) Harris, J.; Yung, K. Carbon Dioxide’s Liquid-Vapor Coexistence Curve and Critical Properties as Predicted by a Simple Molecular Model. J. Phys. Chem. 1995, 99, 12021−12024. (47) Potoff, J.; Siepmann, J. Vapor-Liquid Equilibria of Mixtures Containing Alkanes, Carbon Dioxide, and Nitrogen. AIChE J. 2001, 47, 1676−1682. (48) Ravipati, S.; Punnathanam, S. N. Calculation of Chemical Potentials and Occupancies in Clathrate Hydrates through Monte Carlo Molecular Simulations. J. Phys. Chem. C 2013, 117, 18549−18555. 588

dx.doi.org/10.1021/jp410306v | J. Phys. Chem. B 2014, 118, 577−589

The Journal of Physical Chemistry B

Article

(49) Sum, A. K; Burruss, R. C.; Sloan, E. D., Jr. Measurement of Clathrate Hydrates Via Raman Spectroscopy. J. Phys. Chem. B 1997, 101, 7371−7377. (50) Nosé, S. A Molecular Dynamics Method for Simulations in the Canonical Ensemble. Mol. Phys. 2002, 100, 191−198. (51) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A 1985, 31, 1695−1697. (52) Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52, 7182−7190. (53) Ciccotti, G.; Ryckaert, J. P. Molecular Dynamics Simulation of Rigid Molecules. Comput. Phys. Rep. 1986, 4, 346−392. (54) Span, R.; Wagner, W. A New Equation of State for Carbon Dioxide Covering the Fluid Region from the Triple Point Temperature to 1100 K at Pressures up to 800 MPa. J. Phys. Chem. Ref. Data 1996, 25, 1509−1596. (55) Jorgensen, W.; Madura, J. Temperature and Size Dependence for Monte Carlo Simulations of Tip4p Water. Mol. Phys. 1985, 56, 1381− 1392. (56) Lu, W.; Guo, H.; Chou, I.; Burruss, R.; Li, L. Determination of Diffusion Coefficients of Carbon Dioxide in Water between 268 and 473 K in a High-Pressure Capillary Optical Cell with in-Situ Raman Spectroscopic Measurements. Geochim. Cosmochim. Acta 2013, 115, 183−204. (57) Vlcek, L.; Chialvo, A. A.; Cole, D. R. Optimized Unlike-Pair Interactions for Water−Carbon Dioxide Mixtures Described by the SPC/E and EPM2 Models. J. Phys. Chem. B 2011, 115, 8775−8784. (58) Miller, S.; Smythe, W. Carbon Dioxide Clathrate in the Martian Ice Cap. Science 1970, 170, 531−533. (59) Falabella, B. J. A Study of Natural Gas Hydrates. Ph.D. Thesis, University of Massachusetts, University Microfilims, Ann Arbor, MI, 1975. (60) Larson, M.; Garside, J. Solute Clustering in Supersaturated Solutions. Chem. Eng. Sci. 1986, 41, 1285−1289. (61) Robinson, D.; Mehta, B. Hydrates in the Propane−Carbon Dioxide−Water System. J. Can. Pet. Technol. 1971, 10, 33−35. (62) Deaton, W.; Frost, E., Jr. Gas Hydrates and Their Relation to the Operation of Natural-Gas Pipe Lines; Bureau of Mines, Helium Research Center: Amarillo, TX, 1946. (63) Ng, H.; Robinson, D. Hydrate Formation in Systems Containing Methane, Ethane, Propane, Carbon Dioxide or Hydrogen Sulfide in the Presence of Methanol. Fluid Phase Equilib. 1985, 21, 145−155. (64) Unruh, C.; Katz, D. Gas Hydrates of Carbon Dioxide-Methane Mixtures. Trans. AIME 1949, 186, 83−86. (65) Adisasmito, S.; Frank, R.; Sloan, E. Hydrates of Carbon Dioxide and Methane Mixtures. J. Chem. Eng. Data 1991, 36, 68−71. (66) Ohgaki, K.; Makihara, Y.; Takano, K. Formation of CO2 Hydrate in Pure and Sea Waters. J. Chem. Eng. Jpn. 1993, 26, 558−564.

589

dx.doi.org/10.1021/jp410306v | J. Phys. Chem. B 2014, 118, 577−589