Ab Initio Intermolecular Potentials for Gas Hydrates and Their

5722

J. Phys. Chem. B 2002, 106, 5722-5732

Ab Initio Intermolecular Potentials for Gas Hydrates and Their Predictions Jeffery B. Klauda and Stanley I. Sandler* Center for Molecular and Engineering Thermodynamics, Department of Chemical Engineering, UniVersity of Delaware, Newark, Delaware 19716 ReceiVed: September 20, 2001; In Final Form: February 21, 2002

Interaction energies between the guests CH4, C2H6, C3H8, and CO2 with hydrate cages are calculated using ab initio quantum mechanics (QM), and then fit to a Lennard-Jones (LJ) potential plus an electrostatic term to predict guest occupancies in the hydrate cavities. This procedure differs from that commonly used with the van der Waals and Platteeuw (vdWP) model in which the intermolecular potential parameters are fit to hydrate equilibrium pressures. The procedure here is seen to improve upon the vdWP model with Sloan’s Kihara potential parameters when compared to Raman spectroscopy and NMR data for guest compositions in the hydrate cavities. The LJ potential fit to QM interaction energies is also used in our recently developed fugacitybased model to predict hydrate equilibrium pressures. For the ethane hydrate below 287.4 K and the ethanecarbon dioxide mixed hydrate, the percent absolute average deviation (%AAD) from experimental hydrate equilibrium pressures is only 3.4 and 3.5, respectively. This is an improvement upon Sloan’s vdWP model that has a 10.8 and 4.8%AAD, respectively for these systems.

Introduction Gas hydrates or clathrates are a class of compounds in which water forms a continuous and known crystal structure with small cavities or cages. Compounds such as methane and ethane fill the hydrate cages and stabilize the water lattice. Sir Humphrey Davy1 in 1811 was the first to document the existence of a gas hydrate, specifically a chlorine hydrate. Since Davy’s initial discovery, there have been over 130 hydrophobic compounds found as guests in hydrates.2 Hydrates can form at conditions above the normal freezing point of water at sufficiently high pressure, e.g., methane hydrate is known to be stable up to 330.5 K and 1.5 GPa, above which ice VI is more stable.3 The focus of this research is on the hydrophobic guests CH4, C2H6, C3H8, and CO2 that are commonly found in nature and industry, and form two distinct hydrate crystal structures. Structure I (sI) is a body-centered cubic structure4 with a lattice parameter of 12.03 Å when ethylene oxide is the guest. The cubic cell contains 46 H2O molecules, two 12-hedra (512), and 6 large voids (51262), where 512 is used to indicate that the polyhedron contains 12 five-member ring faces.2 Structure II (sII) is a body-centered cubic structure5 with a lattice parameter of 17.31 Å for a double hydrate when tetrahydrofuran and hydrogen sulfide are guests. Each cubic cell contains 136 water molecules, eight large voids (51264), and 16 12-hedra. Small variations of the lattice parameter occur depending on the guests. Recently there has been considerable interest in naturally forming gas hydrates and methods to harvest their potential methane reserves. Naturally forming hydrates are found in the seafloor and permafrost regions, with a majority in the seafloor.6 It has been estimated that there may be 2 orders of magnitude more methane gas (expanded to atmospheric conditions) in ocean hydrates than in onshore reserves.2 The melting and dissociation of these oceanic and permafrost hydrates have been * Author to whom correspondence should be addressed. Phone: (302) 831-2945. Fax: (302) 831-4466. E-mail: [email protected].

proposed as contributors to global warming.7,8 Alternatively, it has been suggested that hydrates can be used to sequester carbon dioxide in the ocean to reduce its impact on global warming.9,10 Also the petroleum industry is concerned with the formation of large hydrated masses in natural gas pipelines that can slow or completely obstruct flow. We have previously developed a model to predict three-phase equilibrium in single component gas hydrates11 that eliminates the need to empirically fit intermolecular potential parameters as is commonly done in the van der Waals and Platteeuw (vdWP) hydrate model.12 We have also demonstrated that the stability of the various guests in the hydrate cavities can be determined using ab initio quantum mechanics (QM). Here we present a systematic method for calculating the interaction energies between guests and the hydrate lattice at varying guest positions and orientations using ab initio QM. We then fit these interaction energies to the Lennard-Jones 6-12 + Coulomb potential to predict each guest cage occupancy as a function of temperature and pressure, and to predict hydrate equilibrium pressures. Hydrate Models and Intermolecular Potentials between Guest(s) and the Water Lattice The disadvantage of the vdWP model is the required fitting of intermolecular potential parameters to hydrate equilibrium pressures of both pure and mixed gas hydrates. The ability of the model to accurately predict the phase behavior beyond the range of the fit can be poor, e.g., structural changes in mixture hydrates not previously known (as in the methane-ethane hydrate)13 cannot be predicted unless new potential parameters are used. Different sets of potential parameters lead to comparable correlations of measured hydrate phase behavior, but different extrapolations.2 In addition, while the correlation of equilibrium pressures may be satisfactory, the prediction of cage occupancies of guests based on these same parameters can be inaccurate.

10.1021/jp0135914 CCC: $22.00 © 2002 American Chemical Society Published on Web 05/10/2002

Intermolecular Potentials for Gas Hydrates

J. Phys. Chem. B, Vol. 106, No. 22, 2002 5723

The Langmuir absorption model is used to compute the fraction of cages occupied by a guest

fHw (T,P)

Cml(T) fl(T,P)

Θml(T,P) )

1+

∫∫

∑j Cmj(T) fj(T,P)

[

]

(2)

where r is the position vector and Ω is the orientation vector of the guest in the cage, and wgl(r,Ω) and wgg(r,Ω) are the guest-hydrate lattice (gl) and guest-guest (gg) interaction energies, respectively. The assumption contained in eq 2 is that the partition function of the internal modes of the guest in the hydrate phase is same as that in an ideal gas. In the vdWP model eq 2 is simplified to a one-dimensional integral by spherically averaging the water locations, and then a three-parameter Kihara potential is used to represent each guest-cage interaction ignoring guest-guest interactions.2 These potential parameters were fit to measured hydrate equilibrium pressures on the basis of the equality of the chemical potential of water in the hydrate phase (H) and in the other equilibrium phase(s), π, which depending on the conditions might be liquid water, ice, or both, that is

µβw - µHw ) ∆µHw ) ∆µπw ) µβw - µπw

(3)

Here µβw is the chemical potential of the empty hydrate, which is a metastable or hypothetical phase since a hydrate requires a guest to be stable. The difference between the chemical potential of water in the hypothetical and real (filled) hydrate phases is given by

νm ln(1 - ∑ Θmj) ∑ m j

(4)

where νm is the number of cages of type m per water molecule in the hydrate lattice. The chemical potentials of water in the empty hydrate and in the bulk liquid or ice in the vdWP model are related by

∆µπw(T,P) ∆µw(To,Po) ) RT RTo

∆Hπw(T) dT + o RT2 π P ∆Vw dP - ln(γwxw) (5) Po RT

∫TT

∫

where ∆µw(To,Po) is the reference chemical potential at the reference temperature, To, and pressure, Po, usually taken to be the freezing point of water and zero pressure, respectively. The other terms on the right-hand side of eq 5 correct for temperature and pressure differences from the reference state, and for the activity of impure water. To obtain the equilibrium pressure for the hydrate phase with liquid water or ice at a given temperature, one solves eqs 4 and 5 for P by iteration until the equality condition of eq 3 is satisfied. Instead of using eq 5, our single guest hydrate fugacity-based approach11 uses the following equilibrium condition:

fHw (T,P) ) fπw(T,P)

(

)

- ∆µHw (T,P) exp RT

(7)

where ∆µHw (T,P) is given by eq 4 and fβw(T,P) is the fugacity of the hypothetical empty hydrate lattice given by

wgl(r,Ω) + wgg(r,Ω) exp dr dΩ kT

∆µHw (T,P) ) - RT

)

fβw(T,P)

(1)

where fl(T,P) is the fugacity of guest l in the vapor or liquid phase. The Langmuir constant, Cml, in the model can be computed from

1 Cml(T) ) RT

in which the fugacity of water in the hydrate phase is

(6)

fβw(T,P) ) Psat,β w (T) exp

(

)

Vβw(T,P)(P - Psat,β w (T)) RT

(8)

where the temperature dependence of the empty hydrate lattice molar volume, Vβw(T,P), is obtained from Tse14 and the effect of pressure is given by a correlation11 based on experimental X-ray diffraction data.15 The vapor pressure of the empty hydrate lattice was modeled as β ln(Psat,β w ) ) A ‚ln(T) +

Bβ + Cβ + Dβ‚T T

(9)

where the constants Aβ, Bβ, and Dβ in eq 9 were allowed to vary for each guest and hydrate structure, and Cβ is a constant. Further details of this fugacity-based model can be found in ref 11. Procedure for Evaluating Interaction Energies from QM. As a result of the development of quantum mechanical software and with current computer speed, it is feasible to calculate accurate intermolecular interaction energies. The most accurate method to evaluate the integral of eq 2 is to use a high level of QM theory and a large basis set to calculate the interaction energies of a guest at many locations and orientations in a cavity, but this is computationally intensive. Also the extension to differing cavity dimensions or to water molecules beyond the initial cage would require many additional and computationally intensive QM calculations. An alternative that reduces the number of calculations is to use an intermolecular potential model with parameters fit to the results of some QM calculations. Accurate potentials have been obtained by calculating many ab initio interaction energies and using a large number of potential parameters.16 However, Rowley and Pakkanen17 have shown that reasonable fits to ab initio interaction energies can be obtained using simple atom site potentials, for example, using the modified Morse potential. We will show that this approach can be used for guest-water interactions in hydrates without a serious loss of accuracy. In quantum mechanical calculations there are two general sources of error: error due to the level of electron correlation used, and the error in using an incomplete basis set. Since exactly solving the Schro¨dinger equation is possible only for very small molecules, approximate methods that include electron correlation have been developed for larger molecules. The couple cluster method, CCSD(T), developed by Raghavachari et al.18 is the most accurate of these approximate methods and is frequently used as a representation of the exact solution of the Schro¨dinger equation.19 However, as CCSD(T) calculations scale approximately as the number of basis sets to the sixth power, it is not possible to use this method for a large number of molecules as in a hydrate. For hydrogen-bonded and weakly bound molecules (binding energies greater than a few tenths of a kcal/mol), the results of second-order Møller-Plesset perturbation theory (MP2) match reasonably closely with the results obtained using CCSD(T) with the same basis set.19 The error due to an incomplete basis set can be reduced by including additional basis functions. A systematic method that reduces both of these errors will be presented.

5724 J. Phys. Chem. B, Vol. 106, No. 22, 2002

Klauda and Sandler

TABLE 1: Methane-Water Interaction Energies at Various Basis Sets wgw(rmin) [kJ/mol] basis set

HF

MP2

MP3

MP4(SDQ)

CCSD(T)

6-31 g(d) 6-311 g(d) 6-31+g(d) 6-31+g(d,p) 6-311+g(d)

-0.099 -0.134 -0.122 -0.111 -0.175

-0.525 -0.561 -0.651 -0.688 -0.680

-0.517 -0.571 -0.652 -0.673 -0.692

-0.488 -0.538 -0.630 -0.648 -0.668

-0.515 -0.576 -0.681 -0.705 -0.726

Computing the interaction energies using the complete hydrate crystal lattice is also not possible at present, because of the number of basis functions that would result. Therefore, we considered only the three basic cages (512, 51262, and 51264) formed in the sI and sII hydrates. The geometries were optimized to each guest-cavity pair since previous experimental20,21 and computational work22,23 have demonstrated that the crystal lattice size can depend on the guest. As an initial guess, the locations of the water cage molecules were taken from crystallography measurements. The guest-cavity geometry was then optimized at the Hartree-Fock (HF) level of theory with a double-ζ basis set including polarization functions on non-hydrogen atoms, 6-31 g(d), using the Gaussian 9824 suite of programs. This optimized geometry was then used for calculating guest-cavity interaction energies with a more accurate electronic correlation method. Since CCSD(T) calculations are very computationally intensive for the guest with the entire cavity, the interaction energy between one water molecule and a methane molecule at the optimized cavity minimum energy distances was used as a test of the accuracy of the MP2 level of theory. Since the interaction energy is the difference between the energy of the molecular complex and that of the individual guest and the water molecules, the problem of basis set superposition error (BSSE) arises.25 This is reduced in the energy calculations of the individual guest or water molecule by including the basis set of the other molecule in the calculation, but leaving “ghost” atoms in its place. A comparison between the HF, MP2, MP3, MP4(SDQ), and CCSD(T) levels of theory, correcting for the BSSE, is presented in Table 1 for various modest size basis sets. The HF method was the worst, with an 80% absolute average deviation (%AAD) from the more accurate CCSD(T) result. Including electron correlation using MP2, MP3, and MP4(SDQ) improved the predictions dramatically, reducing the errors to 3.5, 2.9, and 7.0% AAD, respectively. Although the MP3 results were slightly closer to the CCSD(T) result than MP2, the difference in accuracy was minimal, while the computational cost was much greater. Therefore, on the basis of this comparison and more generally with the results presented by Dunning, we concluded that MP2 was an adequate level of theory for the present application.19 The error due to an incomplete basis set must also be reduced. However, too much computation time and memory are required for interaction energy calculations for the guest and the entire cavity at the appropriate basis set size. For the case considered here, we have found that calculations of the interaction energies for the full cage to be computationally possible only at MP2/ 6-31 g(d), MP2/6-31+g(d), and MP2/6-31+g(d,p) levels. Alternatively, the cage of water molecules can be broken into sections and the interaction energies between the guest and each section (consisting of numerous water molecules) can be calculated. The interaction energies for each section are then summed to estimate the interaction energy between the guest and the complete cage. Because at fixed position, the orientation of the guest relative to each section is different, we did not

Figure 1. Effect of the basis set size on QM interaction energies with varying types of cage cutting of the 512 cavity and methane at T ) 273 K. All calculations are for the minimum energy structure with the full cage at the HF/6-31 g(d) level.

merely calculate the energy for just one section and multiply that by the number of sections. Initially, the cage was split into two parts (10 water molecules each), so that the basis set size could be increased to 6-31+g(2d,pd) while requiring only modest computational time. In subsequent calculations sectioning into 4 parts (5 water molecules) and 20 parts (single water molecule) were also examined, as was the use of larger basis sets. In Figure 1, a comparison is made of the calculated energies of CH4 in the 512 hydrate cavity for various numbers of cavity sections and different basis sets using the geometry of the minimum energy structure obtained with HF/6-31 g(d) and the complete cavity. There it is seen that the value for the interaction energy between the guest and the sum for the two halves of the cage match well to that of the full cage at MP2/6-31+g(d,p). However, when splitting the cage into two, calculations are possible only with the modest size basis set of 6-31+g(2d,pd), which is far from the basis set limit. For the various basis sets considered, splitting the cage into four sections leads to interaction energies that closely follow the results for both the whole cage and for the two-half cages. As shown in Figure 1, the calculated interaction energies of the 4-section cage appear to be approaching an asymptote for any basis set larger than MP2/6-31+g(3d,3p). However, the sum of the computed interaction energies of methane with each water molecule separately (that is splitting the cage into 20 parts) is far from the results for the full, half, and quarter cages for all but the inaccurate minimally sized basis set, indicating many-body effects are important in hydrates. Cao et al.26 calculated 18000 ab initio interaction energy points between methane and a single water molecule to model the interaction in hydrate cavities. However, on the basis of the results shown in Figure 1, it appears that the delocalization of electrons along a hydrogen bond that occurs in hydrates is important in calculating the interaction energies between the guest and the hydrate. Consequently, while the work of Cao et al.26 undoubtedly leads to an accurate methane-water pair potential, we believe its use for the interaction of a methane molecule in a hydrate cavity leads to an under-prediction of guest-host interaction energies. Although it is known that the computed total energy of a molecule is most accurate at quintuple-ζ basis sets or higher,19

Intermolecular Potentials for Gas Hydrates

J. Phys. Chem. B, Vol. 106, No. 22, 2002 5725

we expect that the interaction energies of guests in hydrates may be less dependent on the number of contracted basis functions in the valence shell. For example, the third and fourth basis sets in Figure 1 indicate that going from a double-ζ basis set (631++g(d)) to a triple-ζ basis set (6311++g(d)) results in only a 4% energy change. For the accurate calculation of the interaction energy between two molecules a large basis set is required,19,26 however we are calculating the interaction between methane and five water molecules. This provides a diffuse electron density since an electron from a water molecule has a finite probability to be located in the basis set of an adjacent water molecule. In a sense, this produces a result similar to the approach frequently used in two-molecule calculations in which basis functions are located between the centers of mass of the molecules to increase the electron density distribution without a large increase in computational cost. Cao et al.26 found that a quintuple basis set (cc-pV5Z) adds about 50% to the energy of the modest size double-ζ basis set (6-31++g(2d,2p)), and therefore a correction to the smaller double-ζ basis to approximate the quintuple basis set was used. However, since we include the five water molecules in our calculations, this results in more accurate interaction energies with a smaller basis set for each molecule. Therefore we believe using MP2/6-31+g(3d,3p) provides a reasonable compromise between accuracy and computational cost, and this is what was used here. The procedure we have used can be generalized to all guestcavity combinations using the following steps and criteria: (1) Optimize the geometry guest-cavity complex to determine the minimum energy configuration at the HF/6-31 g(d) level. (2) Evaluate the interaction energies between the guest and the complete cavity using the geometry in Step 1 at only the MP2/6-31 g(d) and MP2/6-31+g(d) levels, because the computational cost at MP2/6-31+g(d,p) is too high for larger cavities and guests. (3) Split the cavity into sections in such a way (perhaps by symmetry or by trial-and-error) that the sum of the interaction energies of the guest and the sections closely matches that in Step 2. (4) Then calculate the interaction energies of the guestcage at various positions and orientations using MP2/6-31+g(3d,3p). Fits to Potential Models. The procedure above was used to compute the interaction energy between the guest and the water molecules in the cavity at a single fixed location. However, to accurately represent the interaction energies between the guest and the entire crystal lattice, not only should the first shell (or cavity) be considered, but also water molecules beyond the first shell.27-29 Since it is not computationally feasible to calculate the interaction energy between the guest and the entire crystal lattice at all positions and locations, we used a potential model to represent the guest-host interaction. The interaction energy between the guest and the crystal lattice in eq 2 can be separated into two parts,

wgl(r,Ω) ) wgc(r,Ω) + wgl-gc(r,Ω)

(10)

where wgc(r,Ω) is the interaction energy between the guest and the first shell lattice, and wgl-gc(r,Ω) is the interaction energy between the guest and the lattice beyond this first shell. For the interaction energy between methane and the 512 cage, wgc(r,Ω), a total of 86 guest positions and orientations were calculated at the MP2/6-31+g(3d,3p) level. These data were then fit to several potentials using the following Boltzmann factor-weighted objective function,

[ (

∑i

)

- w(ri,Ωi)

# data points

χ)

exp

kT

- exp pred

(

)]

- w(ri,Ωi) kT

2

QM

(11) where χ is minimized by changing parameters in the potential model. Two potentials were considered: the Lennard-Jones (LJ) 6-12 potential

∑ i)1

∑ j)1

[( ) ( ) ] σij

# guest sites # water sites

w(rij) )

4ij

σij

12

-

rij

6

rij

(12)

in which the adjustable parameters are the characteristic energy ij, and the soft core radius σij; and the Kihara potential

w(rij) ) ∞, # guest sites # water sites

w(rij) )

∑ i)1

∑ j)1

[( ) ( )] σij

4ij

rij e 2aij

12

rij - 2aij σij

rij - 2aij

-

6

, rij > 2aij (13)

in which ij, σij, and aij are the adjustable parameters. The parameters that best fit the QM energies are given in Table 2. In this representation we have treated each water molecule as having a single interaction site located on the oxygen atom, and each CH2, CH3, or CH4 group as a single site located on the carbon atom. Even though the Kihara cell potential has an additional parameter, a somewhat better fit was obtained with the LJ potential. A comparison was made of the first water shell contribution to the Langmuir constant computed using the fitted potentials with that using the 86 QM-calculated energies. Most of the points used are in the attractive well region, with only about 16% of the points in the repulsive region. The integral in eq 2 was evaluated numerically using Simpson’s Rule in spherical coordinates. Small spherical shell volumes of 0.1 Å width, ∆Vi, were constructed and the number of interaction energy points that lie within each shell, NQM,i, were weighted as follows # spherical shells NQM,i

Cml(T) )

∑ i)1

( )

wij ∆Vi

exp ∑ kT N j)1

(14)

QM,i

This was then repeated using the LJ potential with the parameters in Table 2 for the positions and orientations of the guest in the QM optimized cage geometry. As shown in Table 4, the calculated Langmuir constants using both the simple LJ water molecule site-methane molecule site and the water atom site-methane molecule site potentials match the Langmuir constants obtained from the direct use of quantum-mechanically calculated energies to within 0.5 and 0.1% absolute deviation (%AD), respectively. The LJ potentials with interaction sites only on the carbon and oxygen atoms, and the LJ H2O site potential with additional sites on the hydrogen atoms of water both fit the QM interaction energies and Langmuir constant reasonably well, but the latter potential has two more parameters than the former; for this reason the simple LJ molecule site potential was used in all later calculations. Table 4 compares values of the first shell contribution to the Langmuir constant computed using various potentials and parameters. Sloan’s2 parameters for the Kihara cell potential were fit to hydrate equilibrium pressures using the vdWP model, but give a 440% AD in the Langmuir constant compared with

5726 J. Phys. Chem. B, Vol. 106, No. 22, 2002

Klauda and Sandler

TABLE 2: Potential Parameters (Interaction sites are labeled in bold.) potential

interaction sites

/k (K)

σ (Å)

LJ Kihara LJ CH4 atom site

CH4-H2O CH4-H2O CH4-H2O CH4-H2O H2O-CH4 H2O-CH4

123.062 126.900 123.976 1.494 113.574 33.014

3.501 3.119 3.451 3.207 3.525 2.236

LJ H2O atom site

a (Å) 0.383

TABLE 3: Potential Fits and Comparison to Langmuir Constant (C) at 273 K Considering Only the First Shell C (Pa-1)

potential LJ Kihara LJ CH4 atom site LJ H2O atom site from QM

no. of C %AD R2 fit of R2 fit of w/kT parameters from QM eq 11

6.994 × 10-7 6.768 × 10-7 7.155 × 10-7

0.5 2.7 2.9

0.778 0.736 0.780

0.914 0.848 0.884

2 3 4

6.963 × 10-7

0.1

0.776

0.914

4

6.957 × 10-7

Figure 2. Boltzmann averaged scaled interaction energies at T ) 273 K between water molecules and methane at varying radii from the center of the cavity.

TABLE 4: Langmuir Constant for Common Potentials (Pa-1) with % AD from QM Prediction empirical parameters

QM cage

% AD

Sloan 2nd virial Kihara 2nd virial LJ OPLS-SPC

3.775 × 10-6 2.087 × 10-6 7.240 × 10-7 6.687 × 10-7

442.6 200.1 4.1 3.9

that computed from our QM energies. However, since Sloan only considered the first hydrate shell in his calculations, the effect of water molecules beyond the first shell has been empirically included in his constants by fitting the potential parameters to measured hydrate equilibrium pressures.2 Also, the set of Kihara potential parameters obtained in this way are not unique, and other sets may fit the hydrate phase behavior equally well, but lead to different values for the Langmuir constant.2 For comparison the LJ and Kihara potential parameters obtained from second virial coefficients30,31 and the following combining rules

σgw )

σg + σw ag + aw ; gw ) xgw; agw ) (15) 2 2

have been used to obtain parameters for eqs 12 and 13, where the subscripts g and w represent the guest in the cavity and water, respectively. The Kihara potential parameters obtained by fitting the second virial predict a Langmuir constant that is too large by 200%, while the LJ potential obtained the same way matches the QM prediction reasonably with only a 4.1% AD. The OPLS32(for methane)-SPC33(for water) potentials commonly used in simulations were also used with the combining rule in eq 15 and produced only a 3.9% AD from the QM prediction. The HF/6-31 g(d) optimized geometry for the 512 cage with methane used to fit the LJ potential was similar to that obtained from crystallography measurements,34 but the average radius

of the cage was slightly larger than that measured from X-ray diffraction. This slight discrepancy may be due to the small basis set used or not including water molecules beyond the first cage. Therefore, in the calculations that follow we used oxygen placements obtained from X-ray diffraction. It has been claimed that the contribution to the guest interaction energy from water molecules beyond the first shell can be significant.27-29 Therefore, we also used the fitted potential to estimate contributions to total interaction energy from beyond the first shell; this differs from the assumption made in most vdWP calculations that only the first shell need be considered.2,12,35 Calculations for three shells with coordination numbers of 20, 20, and 60 water molecules and the average radii given in Table 5 were used with the oxygens placed in their measured positions from methane hydrate crystallography,34 and the effect on the guest-water lattice interaction energy is shown in Figure 2. These results are for the guesthost interaction energy at varying distances from the center of the 512 cavity in Figure 2 are a Boltzmann factor average of 10 positions at each distance from the center of the cavity. The effect of adding the second shell of water molecules, for a total of 40 water molecules, increases the dimensionless energy at the center of the cavity from -8.8 (for 20 water molecules) to -9.6 (40 water molecules), and further to -10.2 when the interaction of methane and 100 water molecules is included. This 14.6% increase in energy has a large effect on the Langmuir constant at 273 K, which is 2.93 × 10-7 Pa if only the first shell is considered, but is 1.12 × 10-6 Pa, a 280% increase, when the first three shells are considered. This shows that even a small change in the total interaction energy can have a large effect on the Langmuir constant. Therefore, the number of water molecules in all three shells given in Table 5 was used for other cages and hydrate structures to accurately represent the interaction energy between the guest and the crystal lattice.

TABLE 5: Shell Radii at T ) 248.15 K average shell radius in Å C2H4O/CH4/CO2 for sI

shell coordination number

structure

cage

1st

2nd

3rd

1st

2nd

3rd

I

512 51262 512 51264

3.906/3.887/3.893 4.326/4.282/4.297 3.902 4.682

6.672/6.609/6.632 7.104/7.044/7.066 6.672 7.472

8.789/8.699/8.724 8.908/8.828/8.853 8.669 9.311

20 24 20 28

20 24 20 28

60 60 60 60

II

Intermolecular Potentials for Gas Hydrates

J. Phys. Chem. B, Vol. 106, No. 22, 2002 5727

TABLE 6: LJ Potential Parameters (Interaction sites are labeled in bold) 512 guest

interaction sites

51262

51264

/k (K) σ (Å) /k (K) σ (Å) /k (K) σ (Å)

CH4 CH4-H2O 123.062 3.501 123.747 3.512 104.183 3.497 CO2 CO2-H2O 103.976 3.48 120.664 2.716 120.365 2.629 CO2-H2O 17.211 3.129 78.482 3.071 95.144 3.032 C2H6 (CH3)2-H2O 114.328 3.523 102.937 3.459 82.630 3.635 C3H8 (CH3)2CH-H2O 99.849 3.664 (CH3)2CH2-H2O 62.544 3.568

Another common approximation used in the vdWP models is that the water molecules are averaged or smeared over a sphere of a given radius, rather than using their known locations from X-ray diffraction.2,12,35 The Langmuir constant calculated for methane by spherically averaging the water molecules with the radii given in Table 5 are 3.46 × 10-7 Pa if only the first shell is considered and 1.21 × 10-6 Pa when the first three shells are considered, resulting in deviations of 18.1 and 8.0% AD, respectively, from the results if each water location is explicitly considered. Consequently, the effect of this assumption is less significant than that of neglecting higher shell interactions. For the 86 QM calculated interaction energies for methane in the 512 cavity, each position or orientation required an average CPU time of 1 day at the MP2/6-31+g(3d,3p) level on a Silicon Graphics Origin 2000 with 195 MHz MIPS R10000 processors, so the complete set required three months of CPU time. Since this is the smallest cage and guest, the computational time required for other guests or cages increases significantly as the CPU time scales approximately as O2N3 where O is the occupied orbitals and N is the number of basis functions. Therefore, it would be beneficial if the number of guest positions and orientations could be reduced significantly without reducing accuracy. To test this, a subset containing only 15 of the 86 QM energies were used to refit the potential parameters. This subset consisted of the energy minimum of the guest position in the cage and 14 points at random distances and orientations in the attractive part of the potential surface, since these give the largest contribution to the Langmuir constant and the objective function of eq 11. The LJ potential parameters fit to the subset of 15 data points has the same correlation coefficient, R2, as parameters fit to the larger number of QM energies. Also the first shell contribution to the Langmuir constant of the LJ fit to the 15 QM energies is smaller than the previous result by only 3.9%AD. Consequently only 15 energy points were used to fit the potential parameters for the larger guests and cavities. (However, we used the parameters the fit to the 86 points of the methane-512 cage interaction energy as these had already been computed and were slightly more accurate.) We used different potential parameters fit to QM calculations for the three types of cages formed in sI and sII; this also differs from previous models that assumed the same value for all cages and structures. This difference in potential parameters from the QM calculations can be attributed to the different environments in each cage type as a result of the different water orientations and the inability of the water molecules to rotate in their fixed crystal positions. For all guests except carbon dioxide, eq 12 was used to represent the intermolecular potential and the best-fit parameters are given in Table 6. However, for carbon dioxide, due to strong preferential orientations in the small and large cavities, as seen in the X-ray diffraction experiments of Udachin et al.,36 not only were the 15 points chosen a priori used, but additional calculations were preformed to better represent the attractive part of the intermolecular potential. Also, the Columbic

contribution to the interaction energy is significant for carbon dioxide and leads to unrealistic LJ parameter values if ignored. Therefore, the potential used for the carbon dioxide-water interaction was # guest sites # water sites

w(rij) )

∑ i)1

∑ j)1

{ [( ) ( ) ] } σij

4ij

12

-

rij

σij

6

+

rij

qiqj rij

(16)

where qi and qj is the charge on the atoms of the guest and water molecule, respectively, where the values used were the Mulliken charges obtained from the MP2/6-31 g(d) QM calculations. The objective function χ in eq 11 was used to find the best fit of the atom site LJ parameters in eq 16, and the values are given in Table 6. The charge contribution for the other guests was ignored because for the interaction sites used, the charges on the guest are essentially zero. Guest-Guest Interaction Energy. The contribution of the guest-guest interaction energy to the Langmuir constant in eq 2 has been ignored in previous thermodynamic hydrate models. However, both Roger37 and Sparks and Tester27 have shown that the guest-guest interaction energy can contribute up to 10% to the total interaction energy which, as shown above, can greatly affect the Langmuir constant. Since the separation distance between guests in a hydrate is significantly larger than location of the minimum in the two-body potential, the predominate contribution to the interaction energy for the nonpolar guest molecules is the dispersion energy. This dispersion energy is commonly modeled using ascending powers of r-n,

wgg(rij) ) -

C6,ij r6ij

-

C8,ij r8ij

-

C10,ij r10 ij

(17)

where here the values of the dispersion coefficients (C6, C8, and C10) were estimated using the London-Margenau38 formulas with the reported guest polarizability and ionization potentials.39 Since for the large separation distances of the guest-guest interactions, the movement of the guests in the cavities has a minimal effect on the configurational integral in eq 2, and this equation can be simplified to

( )∫∫ [ ] ∫∫ [ ]

Wgg 1 ml Cml(T) ) exp RT kT 1 ) Cgg ml(T) RT

exp -

exp -

wgl(r,Ω) dr dΩ kT

wgl(r,Ω) dr dΩ kT

(18)

where Wgg ml is the interaction energy between the guest in cage m and the guests in the surrounding cages assuming all guests are located at the center of their cages. For methane at the center of the small 512 cavity the dimensionless energy in Figure 2 is lowered by -0.4 to a value of -10.6, however this small change increases the Langmuir constant by 53%. Further, the effect of guest-guest interactions can be more important in mixed gas hydrates, because of stronger attractions. For example, the dimensionless energy for methane at the center of the small cavity of the methane-propane hydrate at 193 K is decreased by -0.97 as a result of guest-guest interactions, and the Langmuir constant is increased by 260%. Cage Occupancies Cage occupancies have been reported using Raman spectroscopy and nuclear magnetic resonance (NMR), in which the

5728 J. Phys. Chem. B, Vol. 106, No. 22, 2002 integrated signal intensity ratios of the guests in the two hydrate cavities are measured. Traditionally NMR had been used as an indirect method for measuring cage occupancies, although Sum et al.40 were the first to use Raman spectroscopy to measure cage occupancy and guest vibrational frequency changes in the cavities. Raman spectroscopy has an advantage that it is simpler to use and less resource intensive than NMR. The disadvantage is the inability to quantitatively analyze the Raman spectra because the scattering cross sections, which are proportional to the polarizability of the molecule, may differ for the guest in the gas and in the various hydrate cages.41 Furthermore, as only the ratio of guest occupancies in the cavities can be measured, and since the equality between eqs 4 and 5 is used to determine the absolute cage occupancy, this leads to a dependence of the absolute cage occupancy on whether the vdWP or our model is used, although the difference is small and within experimental error. Equation 1 with the guest fluid-phase fugacity obtained using the PRSV equation of state42 was used to predict the cage occupancies. For methane and carbon dioxide in the sI hydrate, the positions of the water molecules were obtained from X-ray diffraction spectra.34,43 Since there are no experimental data on the water positions for other guests, it was necessary to approximate their distortion of the hydrate lattice. For ethaneforming sI hydrates, the water placement of the ethylene oxide hydrate4 was used since these guests are similar in size and in both cases there is very little occupation of the small cavity. To our knowledge the only diffraction study of sII hydrates is the tetrahydrofuran and hydrogen sulfide double hydrate, so those structural parameters were used for all sII hydrate guests.5 For guests with multiple LJ sites, the orientation of the guest must be included in the evaluation of the Langmuir constant. For each guest center-of-mass position in the hydrate cage, 256 orientations of the guest were used in the integration. The integral in eq 18 was then evaluated using a procedure similar to eq 14, but also including an integral over orientation. The water and carbon dioxide Mulliken atomic charges were obtained from QM calculations of the cage + water complex at the MP2/6-31 g(d) level. Charges of -0.916 and +0.458 on the oxygen and hydrogen atoms, respectively, were the average of the QM-calculated charges for all the water molecules in the 512 cage, and were found to be independent of cage and guest occupation. If the charges on each water molecule are calculated individually (that is, removing the other 19 water molecules but keeping the geometry fixed), the average charge on the oxygen and hydrogen atoms for each of the twenty water molecules in the 512 cage would be -0.835 and +0.4175, respectively. Clearly, the interactions between nearby water molecules significantly affect the charges, and therefore the charges computed for the cage-guest complex, and not for the isolated water molecules, were used. The charge on the carbon atom of CO2 was calculated to be +0.649 and +0.652 in the sI small and large cages, respectively, and for the oxygen atoms on carbon dioxide the average charges were -0.3245 and -0.326 in the sI small and large cages, respectively. Single Gas Hydrates. Since the analysis of NMR spectra requires fewer approximations than Raman spectra, the cage occupancies or their ratios reported from NMR are thought to be more accurate. Ripmeester and Ratcliffe studied methane and other hydrate guest occupations at 193 K using NMR.44 The pure methane hydrate occupancy ratio, θS/θL, from NMR was 0.916 ( 0.01, where S is the small 512 cavity and L is the large 51262 cavity. Using the fitted LJ parameters in Table 6 in eqs 1 and 18 the ratio is 0.941, while using the Kihara potential

Klauda and Sandler

Figure 3. (a) Measured45 and predicted 512 cage occupation of methane at T ) 273.8 K. (b) Measured45 and predicted 51262 cage occupation of methane at T ) 273.8 K.

parameters empirically fit by Sloan2 to calculate the Langmuir constant, which only considers the first shell and averages the water molecules over a sphere, results in a prediction of 0.961 for this ratio. This latter prediction is well outside the reported experimental error bars, while the prediction using the LJ parameters deviate from the NMR results by only 2.7%. If the vdWP model is used to estimate absolute cage occupancies from the experimental data, the occupancy of the large and small cavities are 0.97 and 0.89, respectively. The predictions using our fitted LJ potential are 0.99 and 0.93, while the predictions using Sloan’s fit of the Kihara potential are 0.99 and 0.95. Sum et al.40 and Uchida et al.45 measured cage occupancies for methane hydrates using Raman spectroscopy. Sum et al.40 measured occupancies near the equilibrium pressure while Uchida et al.45 measured occupancies above the equilibrium pressure. From eq 1 as the pressure is increased the cage occupancy should also increase, but this was not what was observed, indicating an inconsistency between the two sets of results. For example, at 273.6 K, Sum et al.40 measured a small cage occupancy of 0.92 at 2.8 MPa but Uchida et al.45 reported 0.72-0.76 for pressures ranging from 3.0 to 6.4 MPa. Without a more accurate measuring technique, i.e., NMR, we cannot determine which of these two data sets is more accurate. Tulk et al.41 have suggested the calibration of Raman spectra with NMR should be used to ensure more accurate values for cage occupancies. In Figure 3 the predictions for the cage occupancies using Langmuir constants obtained from our LJ parameters in Table 6 and Sloan’s empirical fit of the Kihara potential2 are compared

Intermolecular Potentials for Gas Hydrates to the measurements from Uchida et al.45 Our predictions, which include guest-guest (GG) interactions, are within experimental error as seen in Figure 3. The effect of including the guestguest interaction is an increase in the cage occupancy. Sloan’s empirical fit consistently over-predicts the small cage occupancy, as seen Figure 3a. The QM LJ fit has a 16.2% and 0.28% AAD from the measured small and large cage occupancies at the 16 temperatures and pressures, some of which are presented in Figure 3,45 respectively, while Sloan’s Kihara potential has 21.3 and 0.68% AAD, respectively. The reported experimental percent error for these 16 measuremens is on average 22.3 and 2.6 for the small and large cages, respectively, however half of predictions made using the Sloan Kihara potential are outside this experimental error. Mixed Gas Hydrates. Ripmeester and Ratcliffe44 also studied guest occupations for the sII methane-propane hydrate at 193 K using NMR, and reported the occupancy ratio of θL,C3H8/ θS,CH4 as 1.14 and of θL,CH4/θS,CH4 as 0.140. Our LJ fit to QM with guest-guest interactions results in a prediction of 1.07 for the cage occupancy ratio between propane in the large 51264 cage and methane in the small 512 cage, however the Sloan Kihara potential leads to a prediction of 7.54 for this ratio. This large deviation is because of the small value predicted for the Langmuir constant of methane in the small cavity, resulting in a predicted occupation of 0.133 instead of the measured value of 0.780. For the cage occupancy ratio of methane in the large and small cage, both models predict a value of less than 0.001, which deviates significantly from the measured value of 0.14. The measured NMR peak for methane in the large cavity is weak and only slightly distinguishable from the background noise.44 Therefore, there may be significant error in the cage occupancy ratio of methane in the small and large cavity; Ripmeester and Ratcliffe44 did not report an estimate for the experimental error. Sum et al.40 also studied the cage occupancies of the methane-carbon dioxide mixture gas hydrate using Raman spectroscopy. They reported that carbon dioxide appeared to occupy only the large 51262 cavities as the Raman bands did not split into two peaks representing the different environments of the guest in the 512 and 51262 cavities. Recent X-ray diffraction studies of pure carbon dioxide hydrates36,43 found carbon dioxide occupancy in the small 512 cavity. The Raman active intensity vibrational mode of carbon dioxide splits in the gas phase due to Fermi resonance;40 this may alter the interpretation of the splitting of the bands in the hydrate phase and result in the failure of Raman spectroscopy to capture the behavior of guest occupancy. Therefore, the results of Sum et al.40 are questionable because of the assumption that carbon dioxide only occupies the large cage. However, the measured overall composition of the guest in the hydrate does not require a knowledge of individual cage occupancies when only the ratio of overall guest intensities is used; this is what is presented in Figure 4a. For the lower vapor composition of methane, the predictions of both our model and the Sloan Kihara potential are within experimental error; however, for the system with a methane vapor composition of 0.66 the Sloan potential leads to an overprediction of the overall methane composition in the hydrate. Ohgaki et al.46 also have measured the overall composition of the guest in the methane-carbon dioxide gas hydrate. Instead of using spectroscopy, they determined the overall guest composition in the hydrate at a constant temperature of 280.3 K with varying vapor phase compositions (and therefore varying hydrate equilibrium pressures) using the volumetric properties of the vapor, liquid, and hydrate phases and mass balances. The

J. Phys. Chem. B, Vol. 106, No. 22, 2002 5729

Figure 4. (a) Overall composition of methane in the hydrate phase (ZCH4) for the methane-carbon dioxide mixed hydrate at varying temperatures between 273.2 and 278.2 K. (b) Overall composition of methane in the hydrate phase (ZCH4) for the methane-carbon dioxide mixed hydrate at T ) 280.3 K.

reported methane composition in the hydrate phase is shown in Figure 4b. It is evident that predictions using both our model and the Sloan model over-estimate the overall composition of methane in the hydrate phase, and are outside the reported experimental error of 5%. However, similar direct measurements for the hydration number (average number of water molecules per guest in the hydrate lattice) of methane hydrates by several authors have shown significant differences,44 and therefore the accuracy of this method of measurement is uncertain. Subramanian et al.13 measured cage occupancies and structural transitions for the methane-ethane mixture hydrate. They found a previously unknown phase change from sI to sII at methane vapor composition between 0.72 and 0.75, observed as a change in the Raman vibrational frequency and the relative areas under the peaks that correspond to a different environment of the guest molecules.13 Figure 5 displays some of the methane and ethane cage occupancies reported by Subramanian et al.13 In Figure 5a the sI large 51262 cage occupancies, which have been corrected for an error in ref 13, are shown for both methane and ethane.47 Clearly, predictions based on both our QM LJ parameters and the parameter set of Sloan match the experimental data well and similar results for the 51264 cavities in sII are obtained. The importance of guest-guest interactions for the occupancy of methane in the small cage is shown in Figure 5b. For the 512 cavity in both sI and sII, our QM LJ parameters without guest-

5730 J. Phys. Chem. B, Vol. 106, No. 22, 2002

Klauda and Sandler

Figure 6. Equilibrium pressures of the methane sI hydrate versus temperature.

I-H-V region and a small subset of data points in the LwH-V region which gave

ln(Psat,β w [Pa]) ) 4.64130‚ln(T) +

Figure 5. (a) The sI 51262 cage occupancy in the methane-ethane hydrate at T ) 274.2 K versus the composition of methane in the gas phase. (b) The sII methane 512 cage occupancy in the methane-ethane hydrate at T ) 274.2 K versus the composition of methane in the gas phase.

guest interactions led to an under-prediction of the occupancy of methane. However, when guest-guest interactions are included, the occupancy of methane in the small cage is significantly increased, and close to the measured values. Since the vdWP model of Sloan2 does not consider guest-guest interactions, the model may nonetheless predict the correct small cage occupancy of the methane-ethane hydrate because the potential parameters were fit to mixed guest hydrate equilibrium pressures. However, this effective potential leads to overpredictions of the small cage occupancy of methane in the pure hydrate, as the guest-guest interaction energy is less in the pure methane than in the methane-ethane hydrate. Hydrate Equilibrium Pressures Using either the vdWP model or our fugacity-based model, the equilibrium pressure is found by iteration of eqs 3 or 6, respectively, over pressure at fixed temperature until the equality is satisfied. In the vdWP model the Kihara cell potential parameters have been fit to both pure and mixed gas hydrate equilibrium pressures so that the result is a correlation of the available data. Using our fugacity-based model the vapor pressure constants in eq 9 were fit to only pure gas hydrate equilibrium pressures, so that the equilibrium pressures for mixed gas hydrates are predictions. The experimental equilibrium pressures for sI methane hydrate2 are plotted in Figure 6 using the empty hydrate lattice vapor pressure constants (Aβ, Bβ, and Dβ in eq 9) fit to hydrate equilibrium pressures in the

- 5366.10 + 2.77891 + T - 0.008332‚T (19)

Recently, Chou et al.48 using X-ray diffraction have found a structural change of methane hydrate from sI to a hydrate with sII characteristics at 100 MPa and approximately 304 K. For the sI hydrate from 149 to 304 K, the % AAD for the Sloan and our fugacity-based models are 3.95 and 3.28, respectively. The experimental data in Figure 7 for sI ethane hydrate were obtained from Sloan2 for T < 290 K, Nakano et al.49 for 290 K< T < 298 K, and Morita et al.50 for T > 298 K. The empty ethane hydrate lattice vapor pressure equation fitted to the hydrate equilibrium pressures in the I-H-V region and five data points in the Lw-H-V region for this system is

ln(Psat,β w [Pa]) ) 4.76152‚ln(T) +

- 5419.38 + 2.77891 + T - 0.0097740‚T (20)

Our fugacity-based model closely matches the experimental data as seen in Figure 7a; however, the vdWP model with the Sloan parameter set under-predicts the ethane hydrate equilibrium pressures, which may be a result of fitting the parameters to a data set containing both mixed and single gas hydrates. For the three phase equilibrium pressures in the I-H-V and Lw-H-V regions, Sloan’s vdWP model has a 10.8% AAD, while our fugacity-based model has a 3.35% AAD (which is a slight improvement to our earlier fugacity-based model parameter set11 that led to a 3.60% AAD). Also our model matches the measured phase change of ethane from vapor to liquid at 287.4 K, while the prediction of the vdWP model for this phase change occurs st 288.5 K. Both models have errors of about 20% AAD for the slightly scattered high-pressure experimental data in the LwH-LCO2 region. To verify that our fugacity-based model can be extended to mixed hydrates, the equilibrium pressures for the ethane-carbon dioxide hydrate were used. The vapor pressure constants fit to the pure CO2 hydrate equilibrium pressures were found to be

ln(Psat,β w [Pa]) ) 4.59071‚ln(T) +

- 5345.28 + 2.77891 + T - 0.00752178‚T (21)

Intermolecular Potentials for Gas Hydrates

J. Phys. Chem. B, Vol. 106, No. 22, 2002 5731 Because fewer assumptions are made in interpreting NMR spectroscopy data, those results are believed to be more accurate than those from Raman spectroscopy. For the NMR-measured guest occupations in pure methane and methane-propane hydrates, predictions using our QM-based LJ potential are superior to those resulting from parameters that have been empirically fit to the Kihara potential by others. Even considering the lower accuracy of the results obtained from Raman spectroscopy, predictions from our QM-based LJ potential are in reasonable agreement with the experimental measurements. The QM-based LJ potentials for ethane and carbon dioxide were used to predict hydrate equilibrium pressures, and we found that with our parameters, the fugacity-based model leads to more accurate predictions of simple and mixed hydrate equilibrium pressures that the vdWP model with the parameter set of Sloan. Acknowledgment. Financial support of this research was provided from Contract DE-FG02-85ER13436 from the Division of Basic Energy Sciences of the U.S. Department of Energy. The authors also thank the University of Delaware and the National Center for Supercomputing Applications (NCSA) for the use of computers for the quantum mechanical calculations. Nomenclature Nomenclature Roman Characters

Figure 7. Equilibrium pressures of the ethane sI hydrate versus temperature.

For the mixed hydrate the following mixing rule was used for the constants in eqs 18 and 19, # of guests

Xmix )

∑ i)1

Zi X i

(22)

where X is the Aβ, Bβ, or Dβ parameter in eq 9, and Zi is the overall guest composition in the hydrate phase predicted by the Langmuir adsorption isotherm. Using this mixing rule with the parameter values in eqs 20 and 21, our fugacity-base model produces a 3.5% AAD from the 40 ethane-carbon dioxide mixed hydrate equilibrium pressure data points2 compared to a 4.8% AAD using the Sloan vdWP model and parameter set.

a ) Kihara hard core radius C(T) ) Langmuir constant f ) fugacity GG ) guest-guest H ) enthalpy k ) Boltzmann constant LJ ) Lennard-Jones P ) pressure q ) atomic charge QM ) quantum mechanics r ) radial distance r ) positional vector R ) gas constant T ) temperature vdWP ) Van der Waals and Platteeuw w(r) ) intermolecular potential V ) volume x ) molar composition in liquid phase Z ) molar composition in hydrate phase Greek Characters

Conclusions We have presented a new, more accurate method to obtain intermolecular potentials of guest molecules in gas hydrates. Instead of fitting potential parameters to hydrate equilibrium pressures, the potential parameters were obtained from a small set of ab initio calculated intermolecular energies. The LJ potential plus electrostatic term was used to accurately fit the results of the QM calculations, and led to a good prediction of the Langmuir constant. Potential parameters used in most vdWP type models ignore the effect of water molecules beyond the first shell, but we have shown that including them can increase the value of the Langmuir constant by 280%. Also, guest-guest interactions have been shown to make a significant contribution, but these have been ignored in vdWP type models. NMR and Raman spectroscopy methods are currently used to experimentally quantify guest occupations in hydrate cavities.

χ ) objective function ∆ ) change in a property  ) characteristic energy γ ) activity coefficient µ ) chemical potential Θ ) fraction of cage occupation σ ) soft core radius Ω ) orientational vector ν ) number of cages per water molecule in the hydrate lattice Superscripts and Subscripts g ) guest gc ) guest-cage gl ) guest-lattice H ) hydrate phase L ) liquid water

5732 J. Phys. Chem. B, Vol. 106, No. 22, 2002 o ) reference state w ) water β ) empty hydrate phase π ) ice or liquid water References and Notes (1) Davy, H. Philos. Trans. R. Soc. London 1811, 101, 1. (2) Sloan, E. D. Clathrate Hydrates of Natural Gases, 2nd ed.; Marcel Decker, Inc.: New York, 1998. (3) Dyadin, Y. A.; Aladko, E. Y.; Larionov, E. G. Decomposition of methane hydrates up to 15 kbar. MendeleeV Commun. 1997, 34. (4) McMullen, R. K.; Jeffrey, G. A. Polyhedral Clathrate Hydrates. IX. Structure of Ethylene Oxide Hydrate. J. Chem. Phys. 1965, 42, 2725. (5) Mak, T. C. W.; McMullen, R. K. Polyhedral Clathrate Hydrates. X. Structure of the Double Hydrate of Tetrahydrofuran and Hydrogen Sulfide. J. Chem. Phys. 1965, 42, 2732. (6) Kvenvolden, K. A. Potential effects of gas hydrate on human welfare. Proc. Natl. Acad. Sci. U.S.A. 1999, 96, 3420. (7) Hatzikiriakos, S. G.; Englezos, P. The Relationship between Global Warming and Methane Gas Hydrates in the Earth. Chem. Eng. Sci. 1993, 48, 3963. (8) Hatzikiriakos, S. G.; Englezos, P. Permafrost Melting and Stability of Offshore Methane Hydrates Subject to Global Warming. Int. J. Offshore Polar Eng. 1994, 4, 162. (9) Brewer, P. G.; Orr, F. M.; Friederich, G.; Kvenvolden, K. A.; Orange, D. L. Gas hydrate formation in the deep sea: In situ experiments with controlled release of methane, natural gas, and carbon dioxide. Energy Fuels 1998, 12, 183. (10) Brewer, P. G.; Friederich, C.; Peltzer, E. T.; Orr, F. M. Direct experiments on the ocean disposal of fossil fuel CO2. Science 1999, 284, 943. (11) Klauda, J. B.; Sandler, S. I. A Fugacity Model for Gas Hydrate Phase Equilibria. Ind. Eng. Chem. Res. 2000, 39, 3377. (12) van der Waals, J. H.; Platteeuw, J. C. Clathrate Solutions. AdV. Chem. Phys. 1959, 2, 1. (13) Subramanian, S.; Kini, R. A.; Dec, S. F.; Sloan, E. D. Evidence of structure II hydrate formation from methane plus ethane mixtures. Chem. Eng. Sci. 2000, 55, 1981. (14) Tse, J. S. Thermal-Expansion of the Clathrate Hydrates of EthyleneOxide and Tetrahydrofuran. J. Physique 1987, 48, 543. (15) Hirai, H.; Kondo, T.; Hasegawa, M.; Yagi, T.; Yamamoto, Y.; Komai, T.; Nagashima, K.; Sakashita, M.; Fujihisa, H.; Aoki, K. Methane hydrate behavior under high pressure. J. Phys. Chem. B. 2000, 104, 1429. (16) Bukowski, R.; Sadlej, J.; Jeziorski, B.; Jankowski, P.; Szalewicz, K. Intermolecular potential of carbon dioxide dimer from symmetry-adapted perturbation theory. J. Chem. Phys. 1999, 110, 3785. (17) Rowley, R. L.; Pakkanen, T. Determination of a methane intermolecular potential model for use in molecular simulations from ab initio calculations. J. Chem. Phys. 1999, 110, 3368. (18) Raghavachari, K.; Trucks, G. W.; Pople, J. A.; Headgordon, M. A 5th-Order Perturbation Comparison of Electron Correlation Theories. Chem. Phys. Lett. 1989, 157, 479. (19) Dunning, T. H. A road map for the calculation of molecular binding energies. J. Phys. Chem. A 2000, 104, 9062. (20) Berecz, E.; Balla-Achs, M. Gas Hydrate; New York, 1983. (21) Von Stachelberg, M.; Mu¨ller, H. R. Z. Electrochem. 1954, 58, 25. (22) Hwang, M. J.; Holder, G. D.; Zele, S. R. Lattice Distortion by Guest Molecules in Gas-Hydrates. Fluid Phase Equilib. 1993, 83, 437. (23) Zele, S. R.; Lee, S. Y.; Holder, G. D. A theory of lattice distortion in gas hydrates. J. Phys. Chem. B 1999, 103, 10250. (24) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Zakrzewski, V. G.; Montgomery, J. A.; Stratmann, R. E.; Burant, J. C.; Dapprich, S.; Millam, J. M.; Daniels, A. D.; Kudin, K. N.; Strain, M. C.; Farkas, O.; Tomasi, J.; Barone, V.; Cossi, M.; Cammi, R.; Mennucci, B.; Pomelli, C.; Adamo, C.; Clifford, S.; Ochterski, J.; Petersson, G. A.; Ayala, P. Y.; Cui, Q.; Morokuma, K.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Cioslowski, J.; Ortiz, J. V.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Gonzalez, C.; Challacombe, M.; Gill, P. M. W.; Johnson, B. G.; Chen, W.; Wong, M. W.; Andres, J. L.; Head-Gordon, M.; Replogle, E. S.; Pople, J. A. Gaussian 98, Revision A.9; Gaussian, Inc: Pittsburgh, PA, 1998.

Klauda and Sandler (25) Boys, S. F.; Bernardi, F. The Calculation of Small Molecular Interactions by the Differences of Separate Total Energies. Some Procedures with Reduced Errors. Mol. Phys. 1970, 19, 553. (26) Cao, Z. T.; Tester, J. W.; Trout, B. L. Computation of the methanewater potential energy hypersurface via ab initio methods. J. Chem. Phys. 2001, 115, 2550. (27) 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. (28) Sparks, K. A.; Tester, J. W.; Cao, Z. T.; 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. (29) John, V. T.; Holder, G. D. Contribution of 2nd and Subsequent Water Shells to the Potential-Energy of Guest-Host Interactions in Clathrate Hydrates. J. Phys. Chem. 1982, 86, 455. (30) Tee, L. S.; Gotoh, S.; Stewart, W. E. Normal Fluids: The LennardJones 12-6 Potential. Ind. Eng. Chem. Fundam. 1966, 5, 356. (31) Tee, L. S.; Gotoh, S.; Stewart, W. E. Molecular Parameters of Normal Fluids: The Kihara Potential with Spherical Core. Ind. Eng. Chem. Fundam. 1966, 5, 363. (32) Jorgensen, W. L.; Maxwell, D. S.; TiradoRives, J. Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids. J. Am. Chem. Soc. 1996, 118, 11225. (33) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The Missing Term in Effective Pair Potentials. J. Phys. Chem. 1987, 91, 6269. (34) Gutt, C.; Asmussen, B.; Press, W.; Johnson, M. R.; Handa, Y. P.; Tse, J. S. The structure of deuterated methane-hydrate. J. Chem. Phys. 2000, 113, 4713. (35) Parrish, W. R.; Prausnitz, J. M. Dissociation Pressures of Gas Hydrates Formed by Gas Mixtures. Ind. Eng. Chem. Process Des. DeV. 1972, 11, 26. (36) Udachin, K. A.; Ratcliffe, C. I.; Ripmeester, J. A. Structure, composition, and thermal expansion of CO2 hydrate from single-crystal X-ray diffraction measurements. J. Phys. Chem. B 2001, 105, 4200. (37) Rodger, P. M. Stability of Gas Hydrates. J. Phys. Chem. 1990, 94, 6080. (38) Reed, T. M.; Gubbins, K. E. Applied Statistical Mechanics: Thermodynamic and Transport Properties of Fluids, 1st ed.; ButterworthHeinemann: Boston, 1973. (39) Lide, D. R.; Frederiskse, H. P. R. CRC Handbook, 76th ed.; CRC Press: Boca Raton, FL, 1995. (40) Sum, A. K.; Burruss, R. C.; Sloan, E. D. Measurement of clathrate hydrates via Raman spectroscopy. J. Phys. Chem. B 1997, 101, 7371. (41) Tulk, C. A.; Ripmeester, J. A.; Klug, D. D. The application of Raman spectroscopy to the study of gas hydrates. Ann. N.Y. Acad. Sci. 2000, 912, 859. (42) Stryjek, R.; Vera, J. H. PRSV - an Improved Peng-Robinson Equation of State for Pure Compounds and Mixtures. Can. J. Chem. Eng. 1986, 64, 323. (43) Ikeda, T.; Mae, S.; Yamamuro, O.; Matsuo, T.; Ikeda, S.; Ibberson, R. M. Distortion of host lattice in clathrate hydrate as a function of guest molecule and temperature. J. Phys. Chem. A 2000, 104, 10623. (44) Ripmeester, J. A.; Ratcliffe, C. I. Low-Temperature CrossPolarization Magic Angle Spinning C-13 NMR of Solid Methane Hydrates - Structure, Cage Occupancy, and Hydration Number. J. Phys. Chem. 1988, 92, 337. (45) Uchida, T.; Hirano, T.; Ebinuma, T.; Narita, H.; Gohara, K.; Mae, S.; Matsumoto, R. Raman spectroscopic determination of hydration number of methane hydrates. AICHE J. 1999, 45, 2641. (46) 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. (47) Subramanian, S.; Sloan, E. D. Correct Values for Table 5 in Chem. Eng. Sci. 55: 1981 2000. Personal Communication, 2001. (48) Chou, I. M.; Sharma, A.; Burruss, R. C.; Shu, J. F.; Mao, H. K.; Hemley, R. J.; Goncharov, A. F.; Stern, L. A.; Kirby, S. H. Transformations in methane hydrates. Proc. Natl. Acad. Sci. U.S.A. 2000, 97, 13484. (49) Nakano, S.; Yamamoto, K.; Ohgaki, K. Natural gas exploitation by carbon dioxide from gas hydrate fields - high-pressure phase equilibrium for an ethane hydrate system. Proc. Inst. Mech. Eng. Part A-J. Power Energy 1998, 212, 159. (50) Morita, K.; Nakano, S.; Ohgaki, K. Structure and stability of ethane hydrate crystal. Fluid Phase Equilib. 2000, 169, 167.