Configurational Properties of Water Clathrates: Monte Carlo and

To quantify the inadequacies associated with the use of the spherically symmetric Lennard-Jones Devonshire smooth cell approximation, Monte Carlo and ...
0 downloads 0 Views 140KB Size
Download PDF
6300

J. Phys. Chem. B 1999, 103, 6300-6308

Configurational Properties of Water Clathrates: Monte Carlo and Multidimensional Integration versus the Lennard-Jones and Devonshire Approximation Kevin A. Sparks, Jefferson W. Tester,* Zhitao Cao, and Bernhardt L. Trout Energy Laboratory and Department of Chemical Engineering, Room E40-455, Massachusetts Institute of Technology, 77 Massachusetts AVenue, Cambridge, MA 02139 ReceiVed: January 27, 1999; In Final Form: May 10, 1999

To quantify the inadequacies associated with the use of the spherically symmetric Lennard-Jones Devonshire smooth cell approximation, Monte Carlo and multidimensional quadrature integrations were performed for various water clathrate cavity geometries while accounting for the asymmetries of the host lattice using complete crystallographic structural data with effects included for multiple coordination shells. Using the van der WaalsPlatteeuw model as a basis, the configurational partition function was evaluated for typical hydrate hostlattice systems containing small (Ar, CH4, C2H6, CO2, ...; structure I type) and large (C3H6, C4H8, CHCl3, ...; structure II type) guest molecules. A range of intermolecular potential energy and size parameters were used to characterize configurational effects for structure I and II guest-host interactions. Additionally, angleaveraged potential energy profiles were calculated in order to compare with those determined using the LennardJones Devonshire approximation. It was found that the Lennard-Jones Devonshire approximation gives quantitatively correct results for smaller guest molecules (σ > 3.0 Å) but completely incorrect results for larger guest molecules (σ > 3.0 Å). In addition, quadrature calculations must include four coordination shells or more to provide an accurate estimation of physical properties.

Introduction and Motivation Given the practical importance of gas hydrates to the oil and gas industry and their suitability for tractable molecular simulations, these water clathrates have been the subject of several studies in physical chemistry and chemical engineering.1,2 With a well-defined host lattice of water molecules entrapping guest molecules such as CO2 and CH4, they provide an excellent framework for estimating configurational effects using Monte Carlo, molecular dynamics, and even numerical integration (see earlier studies by Tester et al.,3 John and Holder,4-6 Holder and co-workers,7-9 Rodger,10 and Tse and Davidson11). A key feature of clathrate compounds is that relatively weak van der Waals guest-host interactions stabilize the hydrate host lattice. Gas hydrates commonly form in two different structures (I and II) that have the following ideal compositions with fully occupied water cages (see Figure 1 and Table 1).

structure I: M1‚3M2‚23H2O structure II: 2M1‚M2‚17H2O where M1 ) guest molecules occupying small dodecahedral cavities and M2 ) larger guest molecules occupying larger tetrakaidecahedral or hexakaidecahedral cavities. Clathrates are generally regarded as nonstoichiometric compounds since all available cages within the host lattice structure need not be occupied. Given these structural characteristics, classical statistical mechanics provides a means for quantifying the binary interactions between host and guest molecules with high accuracy using the traditional pairwise addivity approximation. * Corresponding author.

This paper reports an improved approach to calculate configurational properties that highlights several limitations observed in earlier work. Theoretical Approach The starting point for our study is the van der WaalsPlatteeuw model.14-16 This model corresponds to a threedimensional representation of localized adsorption of guest molecules into cages formed by the host molecules. The difference between the chemical potential of water in the occupied clathrate lattice (µHw) and empty lattice (µβw) is expressed by the following equation:15,16

∑i νi ln(1 + ∑J CJifˆJ)

(µHw) - (µβw) ) - kT

(1)

where νi is the number of type i cavities per water molecule in the host lattice, fˆj is the fugacity of guest molecule J, and CJi is the Langmuir constant defined for each type J guest component within a type i cavity as

CJi t ZJi//kT

(2)

where ZJi/ is the classical configuration integral

ZJi* t

1 8π2

∫‚‚‚∫ exp[-U(r,θ,φ,R,β,γ)/kT] r 2 × sin θ sin β dθ dφ dr dR dβ dγ (3)

and U ) U(r,θ,φ,R,β,γ) is the total interaction potential between the guest molecule and all host molecules. The integration is 10.1021/jp9903108 CCC: $18.00 © 1999 American Chemical Society Published on Web 07/09/1999

Configurational Properties of Water Clathrates

J. Phys. Chem. B, Vol. 103, No. 30, 1999 6301

Figure 1. (A) Structure I: water clathrate polyhedra space group Pm3n; lattice constant ) 12.03 ( 0.01 Å. (B) Structure II: water clathrate polyhedra space group Fd3m; lattice constant ) 17.31 ( 0.01 Å (see also refs 1 and 22-26).

TABLE 1: Lennard-Jones σ Values and Ideal Hydrate Compositions for Representative Guest Componentsa,b water clathrate guest

σ12 from gas viscosities (Å)

methane ethane ethylene carbon dioxide xenon cyclopropane hydrogen sulfide

3.1995 3.542 3.402 3.291 3.344 3.724 3.132

argon krypton nitrogen oxygen propane cyclopropane I-butane

3.0915 3.148 3.2195 3.054 3.8795 3.724 3.9595

σ12 from virial coefficients (Å)

Structure I 3.3255 3.9305 3.537 3.5285 3.3705

Structure II 3.07 3.2435 3.1675 3.1105 4.139

Thus, eq 1 can be rewritten as

ideal fully occupied hydrate composition 4CH4‚23H2O 3C2H6‚23H2O 4C2H4‚23H2O 4CO2‚23H2O 4Xe‚23H2O 3C3H6‚23H2O 4H2S‚23H2O

3Ar‚17H2O 3Kr‚17H2O 3N2‚17H2O 3O2‚17H2O C3H8‚17H2O C3H6‚17H2O C4H10‚17H2O

a σ ) (σ + σ )/2 (1 is the guest molecule; 2 is the water molecule 12 11 22 where σ22 ) 2.641 Å; σ22 is from gas viscosity data only as virial coefficient data for H2O are not available). b Virial coefficient σ12’s are provided by ref 12. Gas viscosity σ12’s are provided by ref 13.

carried out spatially over configuration space using spherical coordinates (r,θ,φ) and Euler orientation angles (R,β,γ) for the guest molecule. The fraction of available cages of type i occupied by guest molecule k is given by

Ckifˆk

yki ) 1+

∑J CJi fˆJ

∑i νi ln(1 - ∑J yJi)

(µHw) - (µβw) ) +kT

(5)

For various situations involving different phases in equilibrium, such as along three-phase monovariant lines or at four-phase invariant points in binary (single guest type-host) systems, eq 5 is used to constrain the system to determine equilibrium conditions (see, for example, Holder et al.17,18). Since analytical integration of eq 3 is intractable given the asymmetries of the host-guest potential, several approximations have been applied. Most commonly, a Lennard-Jones and Devonshire (LJD) liquid cell model is used to simplify the integration.15,16 The potential is averaged and uniformly distributed on a single spherical surface for each cage which is located at the average radius of that cavity. The averaging in this case is over both spherical angles (θ and φ) and effectively smoothes the potential function U(r,θ,φ) to a dependence on radial distance only (see also Hirschfelder, Curtis, and Bird12). The LJD model expresses the Langmuir constant as a onedimensional integral,

CJi(LJD) t C* )

∫

4π exp[-〈W(r)〉/kT]r 2 dr kT

(6)

where the integration is carried out over the cavity radius and 〈W(r)〉 ) 〈U(r)〉LJD is the angle-averaged, spherically symmetric cell potential defined at each r by

(4) 〈W(r)〉 t

1 4π

∫02π ∫0π U(r,θ,φ) sin θ dθ dφ

(7)

6302 J. Phys. Chem. B, Vol. 103, No. 30, 1999

Sparks et al.

TABLE 2: Cell Radii and Coordination Numbers for the LJD Model for Structure I and II Hydrates cell radius Ri, Å coordination numbers Zc,i John and Holder6 this study John and Holder6

cell

this study

structure I, small cavity structure I, large cavity structure II, small cavity structure II, large cavity

3.905

3.875

20

20

4.326

4.30a

24

21

3.902

3.870

20

20

4.682

4.703

28

28

a One reviewer noted that the actual R value used by John and Holder i was 4.182, not 4.30 as reported in ref 6.

where

[

〈W(r)〉 ) 2Ζc,i

(

)

(

)]

σ12 10 a 11 σ6 a δ + δ - 5 δ4 + δ5 11 R R Ri r Ri r i i

Calculational Results and Discussion We explored several methods to evaluate the guest-host configurational partition function. Initially, the Metropolis Monte Carlo sampling technique21 described earlier was implemented. However, due to the uncertainties of the “free volume” estimates, this method was quickly abandoned. This method did, however, provide an accurate estimate of the ensemble-averaged guesthost potential energy which can be rigorously calculated using the configuration integral given in eq 3. Therefore, the results of the Metropolis MC simulation could be used to compare directly with average potential values resulting from the direct integration of the configuration integral. Several standard integration techniques were used to evaluate guest-host configurational properties; these included (1) simple Monte Carlo integration, (2) composite trapezoidal rule, and (3) GaussLegendre quadrature. The simple Monte Carlo integration scheme is best described as an elementary application of the mean value theorem of integral calculus, which in this context is expressed as

Z/Ji )

with

δN ) [(1 - r/Ri - a/Ri)-N - (1 + r/Ri - a/Ri)-N]/N where Zc,i ) coordination number of cavity i, Ri ) effective cavity i radius, a ) core radius of interaction, σ ) the distance between molecular cores at which there is no interaction, and  ) the depth of the intermolecular potential well. Table 2 provides a listing of the Zc,i and Ri values used in this study. They were based directly on crystallographic data for structure I and II hydrates. Phase equilibrium calculations for most water clathrates or gas hydrates involve evaluations of eq 1 for two basic hydrate structures (I and II) with two types of cells each (see Figure 1 and Tables 1 and 2 for crystallographic information). The majority of research using the van der WaalsPlatteeuw model over the past 30 years has followed this approach employing either a Lennard-Jones 6-12 or a Kihara19 potential to capture specific guest-host interactions using appropriate size (σ and a) and energy () parameters (see also, Parrish and Prausnitz,20 Sloan,1 and Makogon2 for more detailed discussion). Commonly, σ, a, and  parameters are based on second virial coefficient or transport coefficient measurements, or they may be fit to gas hydrate phase equilibrium data. In all cases, they are highly correlated. In an earlier paper,14 we examined the inadequacy of the nearest-neighbor approximation commonly used to evaluate the thermodynamic properties of clathrates using the van der Waals-Platteeuw model.15,16 Specifically, we quantified the importance of including guest-host interactions beyond the first water shell using a lattice summation approach. In addition, interactions between guests in adjacent cages were examined. The magnitude of the errors introduced by including only nearest-neighbor (first-shell) interactions can be seen by comparing the Langmuir constant CJi that includes all water shell guest-host interactions and guest-guest interactions to one that includes only interactions in the first shell. At 273 K, even for small sized guests such as argon, the ratio (CJi (n shells)/CJi (first shell)) ranged from 2.68 to 3.71 for all four cavity types in structures I and II, and for larger guests, such as isobutane (i-C4H10), the ratio increased to 1772 for the hexakaidecahedron cage of a structure II hydrate. Table 1 provides estimates of Lennard-Jones σ values and ideal stoichiometric compositions for several representative guest molecules.

∫Ve-U/kT dV ≈ 〈e-U/kT〉V (

〈

〉

〈e-2U/kT〉 - 〈e-U/kT〉2 N

(8)

where N is the number of random guest configurational samples. The ensemble averages 〈e-U/kT〉 and 〈e-2U/kT〉 are approximated by -U/kT

〈e

〉t

N

1

∑e N k)1

-U/kT

〈e

-2U/kT

〉t

1

N

e-2U/kT ∑ N k)1

(9)

Applying simple Monte Carlo integration to the configuration integral of a spherically symmetric guest molecule results in a significant simplification of eq 2. Here we have chosen not to include the variance term as a matter of convenience.

Z/Ji ≈

2π2R N

N

∑ e-U(r ,θ ,φ )/kTrk2 sin θk k)1 k

k

k

(10)

Similarly, for an asymmetric guest molecule we can write

Z/Ji ≈

π3R N

N

∑ e-U(r ,θ ,φ,R ,β ,γ )/kTrk2 sin θk sin βk k)1 k

k

k

k

k

(11)

The major advantage associated with simple Monte Carlo integration is that it is extremely easy to implement, and its accuracy, proportional to N -1/2, is independent of the dimensionality of the integration. One disadvantage, however, is that it is ill-suited for efficiently estimating positionally averaged potential energy profiles within the different host lattice cavities. Since standard numerical integration methods are typically set up to use a one-, two-, or three-dimensional grid, they are better suited to yield positionally averaged information such as angle-averaged potential energy profiles. This structural characteristic is an advantage in that the resulting potential profiles can be directly compared to those derived from the LJD model (see eqs 6 and 7, for example). The composite trapezoidal rule was evaluated as a method for the estimation of the configurational partition function. However, given the geometric complexities of water clathrate systems, a multi-interval 10-point Gauss-Legendre quadrature formula was found to be a much more efficient integration technique in terms of the number of grid points required for a given level of accuracy. The flexibility associated with the use

Configurational Properties of Water Clathrates

J. Phys. Chem. B, Vol. 103, No. 30, 1999 6303

Figure 2. Angle-averaged potential energy for structure II pentagonal dodecahedron [note for a ) 0 the Kihara and Lennard-Jones 6-12 potentials are identical].

Figure 3. Angle-averaged potential energy for structure II hexakaidecahedron [note for a ) 0 the Kihara and Lennard-Jones 6-12 potentials are identical].

of such a multi-interval (n) integration method,

∫abf(x) dx ) ∫ab/nf(x) dx + ∫b/n2b/nf(x) dx + ... + ∫(nb - 1)b/nf(x) dx

(12)

exceeded the advantages of a single-interval, higher-order formula in that the accuracy of the integration was not restricted by the choice of the integration formula. The actual evaluation of the multidimensional configurational partition function involved the simple repeated application of the one-dimensional, 10-point Gauss-Legendre quadrature formula discussed in detail by Carnahan et al.22 A subdivision of each dimension into 4 intervals with 10 quadrature points each was usually found to return values of sufficient accuracy, usually on the order of 1 part in 10 000. However, for the fiveand six-dimensional integrations associated with asymmetric guest molecules, computational time increased significantly to a point where we commonly included only two or one subdivisions, respectively. For example, in the calculation of the configurational integral of cyclopropane guest-water system, the use of two intervals for each dimension, corresponding to 106 Gauss points, required approximately 360 min of Pentium II (330 MHz) CPU time. A sample angle-averaged potential energy profile of a spherical guest within the structure II pentagonal dodecahedron is shown in Figure 2. Similarly, the angle-averaged profiles of different sized spherical guest within the structure II hexakaidecahedron are shown in Figures 3 and 4. Due to the similar geometry of the various cavities, the structure I water clathrate cavity profiles are not presented here even though calculations were made (see Sparks23). The dashed curves in each of these figures represent the Lennard-Jones Devonshire spherical cell approximation obtained using eq 7 where potential effects of the first shell of water molecules were included. The family of solid curves represents the full three-dimensional integration over the rigid host lattice. The uppermost solid line includes the first-shell interaction only. The five additional solid lines represent the inclusion of subsequent water-shell interactions in the potential energy calculations. The convergence of the total potential energy as a function of these succeeding interactions is clearly illustrated in each of the figures. Again, the inclusion of these additional water shell interactions has a rather significant effect on the potential energy profile and must therefore be considered in the configurational partition function evaluation.

Figure 4. Angle-averaged potential energy for structure II hexakaidecahedron.

As is quite obvious from Figures 2 and 3, for the smaller guest molecules (σ e 3.0 Å), the Lennard-Jones Devonshire approximation does a remarkably good job in describing the potential energy profiles within the different water clathrate cavities. Of course, the addition of the subsequent water-shell interactions would be still be necessary in order to accurately capture their important contributions. For the larger guest molecules (σ > 3.0 Å), as depicted in Figure 4 for σ ) 3.5, the Lennard-Jones Devonshire (LJD) approximation is completely inadequate in describing the potential energies and is therefore incapable of providing reliable estimates of the configurational partition function. Using Table 1 as a guide, this result suggests that for any asymmetric guest molecules larger than oxygen or nitrogen serious errors using the LJD approximation will result. To explore further the inadequacies of the LJD approximation in regard to the asymmetries of the host lattice, we have performed calculations for a wide range of Kihara intermolecular potential parameters following the Q* convention first proposed by John and Holder.6 For a particular cage type J and guest molecule i,

Q* t

CJi(quadrature) CJi(LJD)

)

C C*

(13)

With their definition of Q* we have compared the Langmuir constant determined using the Lennard-Jones Devonshire approximation [CJi(LJD) ) C*] with those determined from the complete three-dimensional quadrature integration over the host

6304 J. Phys. Chem. B, Vol. 103, No. 30, 1999

Figure 5. Langmuir constant ratio Q* results from this study for the structure I pentagonal dodecahedron.

Figure 6. Langmuir constant ratio Q* results from this study for the structure I tetrakaidecahedron.

lattice [CJi(quadrature) ) C], which is expected to provide a substantially higher level approximation than the LJD model. To provide a direct comparison with John and Holder’s results, only first-shell interactions were included in the calculation of CJi (quadrature) ) C. The results for the structure I pentagonal dodecahedron are illustrated in Figure 5. The symbols represent the actual calculations, while the solid lines are simple cubic splines through the data points. The results for the structure I tetrakaidecahedron are shown in Figure 6. Additionally, the results for the structure II pentagonal dodecahedron and hexakaidecahedron are given in Figures 7 and 8. Detailed crystallographic information for gas hydrate structures is available in the literature (von Stackelberg and Muller,24,25 Jeffrey and co-workers,26,27 Mak and McMullan,28 Sloan, Happel, and Hnatow,29 and Sloan1). In our study, the crystallographic parameters developed by McMullan and Jeffrey30 were used for the structure I hydrate lattice and by Mak and McMullan28 for the structure II lattice. Details are provided in Table 3, and the oxygen atom coordinates used in our simulations are tabulated by Sparks.23 It is clearly evident from Figures 5-8 that there exists appreciable deviation between the LJD value and that determined from the precise, three-dimensional quadrature integrations. We know that for the larger, more asymmetric water

Sparks et al.

Figure 7. Langmuir constant ratio Q* results from this study for the structure II pentagonal dodecahedron.

Figure 8. Langmuir constant ratio Q* results from this study for the structure II hexakaidecahedron.

TABLE 3: Crystallographic Parameters Used in This Studya parameter

structure I

structure II

space group lattice constant oxygen fractional position generating functions

Pm3n 12.03 ( 0.01 Å y(k) ) 0.30710 z(k) ) 0.11819

Fd3n 17.31 ( 0.01 Å x(g) ) -0.05744 z(g) ) -0.24487 x(e) ) -0.09228

a Sources: McMullan and Jeffrey,27 Mak and McMullan,28 and Sparks.23

clathrate formers, the use of the Lennard-Jones Devonshire approximation results in estimates of the dissociation pressure which are far below those determined experimentally. The inverse relationship between the Langmuir constant of the guest molecule and its fugacity, as dictated by van der WaalsPlatteeuw model (see eqs 1-5) requires smaller Langmuir constants if we are to better predict phase equilibria such as three-phase (solid hydrate-vapor-ice) dissociation pressures. The decreasing nature of the “numerically approximated, quadrature-determined” Langmuir constants as illustrated in Figures 5-8 are therefore physically reasonable. Another key point illustrated in these figures is that our calculations exhibit completely different behavior than that reported earlier by John and Holder.6 The results of their threedimensional integrations as illustrated in Figures 9-12 instead

Configurational Properties of Water Clathrates

J. Phys. Chem. B, Vol. 103, No. 30, 1999 6305

Figure 12. Structure II hexakaidecahedron results for Q* from John and Holder.6 Figure 9. Structure I pentagonal dodecahedron results for Q* from John and Holder.6

Figure 10. Structure I tetrakaidecahedron results for Q* from John and Holder.6 Note that an “optimized” coordination number ) 21Zc,i and cell radius Ri ) 4.182 Å (incorrectly reported as 4.30 Å) was used by John and Holder.

Figure 11. Structure II pentagonal dodecahedron results for Q* from John and Holder.6

show a nearly opposite trend. For the same range of potential parameters, John and Holder reported an increasing Langmuir constant ratio for increasing values of σ. The magnitude of the changes was also smaller, usually on the order of a factor of 2 or 3. Using somewhat questionable reasoning, they apparently discounted their computed integration results in favor of an empirical correlation for Q*, which instead gave the theoretically correct downward trend similar to the calculated results of our

study. It is quite likely that John and Holder’s integrations are inconsistent. The authors did not document the numerical methods they used nor did they report any checks of their integration procedures. We note that the lattice parameters and oxygen atom positions used by John and Holder as cited by John31 are slightly different than the ones employed in our simulations. While this undoubtedly would have an effect on actual numerical values for CJi, it should not result in different derivative trends with increasing σ. In a later paper, Holder and co-workers32 discuss the possibility that lattice distortions could influence the chemical potential of the occupied lattice and the Langmuir constants in a manner which would show the correct downward trends for larger guest molecules. While this may be plausible, further experimental confirmation is needed. Our results show the correct theoretical dependency without invoking such ad hoc hypotheses. Table 4 provides a direct comparison of results from this study with those reported by John and Holder.6 Here we show the calculated Langmuir constants and Q* ratios in both structure I and II, large and small cavities [Figure 1]. Two quadrature methods for the integration were used: (1) trapezoidal approximation with 10 000 equally spaced points radially outward from the center of the cavity to the outer radius; (2) GaussLegendre approximation employing 10 points per interval with 4 intervals (40 points total) over the same radial dimensions. As can be seen in Table 4, both approximations give the same result. For the LJD approximation in Table 4 only, we used values for the cavity coordination numbers (Zc,i) and optimal cell radii (Ri) recommended by John and Holder6 to provide a meaningful comparison (see Table 2). The Zc,i values used by John and Holder6 were the same as those used in our study except for the larger, tetrakiadecahedral cavity of the structure I system where an optimal Zc,i of 21 was used; however, as indicated, the Ri values were different. In general, John and Holder obtained Zc,i and Ri values for the LJD model by fitting the average 〈W(r)〉 values to those obtained by discrete summation. We used Monte Carlo integrations to test the results of our multidimensional Gaussian-Legendre quadrature and found excellent agreement. The additional potential energy profiles reported previously also provided us with a quantitative check of our methods. In particular, the similarity between the LJD model results and the exact integrations for the smaller guests provided us with a positive verification of our integration methods. Furthermore, by setting the intermolecular potential

6306 J. Phys. Chem. B, Vol. 103, No. 30, 1999

Sparks et al.

TABLE 4: Comparison of Calculated Langmuir Cosntants of Spherical Guest Molecules as a Function of Molecular Parameters for First-Shell Interactions Onlya-c

σ, Å R, Å /k, K

C, atm-1 (John and Holder6)

C*, atm-1 Q* ) C/C*, atm-1 (John and (John and C*(trap),d atm-1 C*(Gauss),e atm-1 C(first shell), atm-1 Q* ) C/C*, atm-1 6 6 Holder ) Holder ) (this study) (this study) (this study) (this study)f

2.50 3.00 3.50 2.50 3.00 3.50 2.50 3.00 3.50

0.20 0.20 0.20 0.30 0.30 0.30 0.30 0.30 0.30

130 130 130 130 130 130 300 300 300

1.11E-02 5.94E-02 1.87E-03 1.25E-02 6.53E-02 5.13E-05 2.28 1.61E+03 4.70E-03

1.14E-02 6.95E-02 1.53E-03 1.33E-02 7.70E-02 3.24E-05 2.87 3.58E+03 1.83E-03

0.968 0.855 1.227 0.939 0.848 1.583 0.792 0.624 2.583

2.50 3.00 3.50 3.75 2.50 3.00 3.50 3.75 2.50 3.00 3.25 3.50 3.75 2.50 2.75 3.00 3.25 2.50 2.75 3.00 3.25

0.25 0.25 0.25 0.25 0.40 0.40 0.40 0.40 0.40 0.40 0.40 0.40 0.40 0.88 0.88 0.88 0.88 0.88 0.88 0.88 0.88

130 130 130 130 130 130 130 130 300 300 300 300 300 300 300 300 300 130 130 130 130

1.40E-02 7.22E-02 1.41E-01 2.68E-03 1.62E-02 9.67E-02 2.55E-02 0 2.08 1.05E+03 1.95E+04 1.74E+03 0 8.04E+01 3.25E+03 2.03E+03 2.13E-05 3.41E-02 7.91E-02 2.56E-02 3.60E-06

1.05E-02 5.98E-02 9.22E-02 8.32E-04 1.25E-02 8.50E-02 1.35E-02 0 1.4 1.07E+03 1.68E+04 5.35E+02 0 1.02E+02 5.43E+03 2.33E+03 3.10E-06 3.36E-02 8.76E-02 2.44E-02 1.40E-06

1.326 1.208 1.524 3.215 1.293 1.138 1.89

2.50 3.00 3.50 3.75 2.50 3.00 3.50 3.75 2.50 3.00 3.25 3.50 3.75

0.20 0.20 0.20 0.20 0.30 0.30 0.30 0.30 0.30 0.30 0.30 0.30 0.30

130 130 130 130 130 130 130 130 300 300 300 300 300

1.19E-02 6.65E-02 1.57E-03 0 1.36E-02 7.31E-02 3.57E-05 0 2.78 2.16E+03 1.80E+03 2.14E-03 0

1.15E-02 7.00E-02 1.31E-03 0 1.34E-02 7.71E-02 2.52E-05 0 2.98 2.69E+03 1.90E+03 1.07E-03 0

1.036 0.949 1.2

2.50 3.00 3.50 4.00 4.25 2.50 3.00 3.50 4.00 4.25 2.50 3.00 3.50 4.00 4.25

0.30 0.30 0.30 0.30 0.30 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50

130 130 130 130 130 130 130 130 130 130 300 300 300 300 300

1.81E-02 8.51E-02 1.21 1.64 4.41E-03 2.10E-02 1.47E-01 2.38 3.28E-02 0 1.52 5.27E+02 4.51E+06 7.02E+03 0

1.61E-02 8.63E-02 1.29 1.62 3.47E-03 1.95E-02 1.52E-01 2.55 2.78E-02 0 1.32 5.95E+02 5.54E+06 5.06E+03 0

1.123 0.992 0.939 1.016 1.272 1.078 0.967 0.935 1.181

1.491 0.981 1.161 3.245 0.79 0.599 0.871 6.87 1.015 0.908 1.049 2.57

1.012 0.948 1.417 0.933 0.805 0.948 2

1.126 0.886 0.814 1.385

Structure I, Small Cavity 1.10E-02 6.90E-02 1.50E-03 1.30E-02 7.70E-02 3.20E-05 2.9 2.60E+03 1.80E-03 Structure I, Large Cavity 9.40E-03 4.30E-02 0.17 2.80E-02 1.10E-02 6.00E-02 9.00E-02 7.70E-04 0.63 2.20E+02 6.70E+03 1.40E+04 1.3 20 1.10E+03 1.40E+04 16 2.30E-02 6.90E-02 9.10E-02 2.00E-03 Structure II, Small Cavity 1.10E-02 7.00E-02 1.30E-03 6.90E-10 1.30E-02 7.70E-02 2.50E-05 3.50E-15 3 2.70E+03 1.90E+03 1.10E-03 1.10E-25 Structure II, Large Cavity 1.60E-02 8.60E-02 1.3 1.6 3.50E-03 1.90E-02 0.15 2.5 2.80E-02 6.00E-09 1.3 5.90E+02 5.50E+06 5.00E+03 1.00E-11

1.10E-02 6.90E-02 1.50E-03 1.30E-02 7.70E-02 3.20E-05 2.9 2.60E+03 1.80E-03

1.50E-02 6.10E-02 1.50E-07 2.00E-02 2.20E-02 3.00E-18 35 9.70E+02 5.00E-32

1.4 0.88 1.0E-04 1.5 0.28 9.40E-14 12 0.37 2.80E-29

9.40E-03 4.30E-02 0.17 2.80E-02 1.10E-02 6.00E-02 9.00E-02 7.70E-04 0.63 2.20E+02 6.70E+03 1.40E+04 1.3 20 1.10E+03 1.40E+04 16 2.30E-02 6.90E-02 9.10E-02 2.00E-03

1.70E-02 0.12 7.30E-03 1.30E-08 2.80E-02 8.70E-02 1.70E-11 1.30E-32 31 1.70E+04 3.3 1.20E-16 0 9.30E-20 0 0 0 6.20E-13 0 0 0

1.8 2.8 4.30E-02 4.6E-07 2.6 1.5 1.90E-10 1.70E-29 49 77 4.9E-04 8.6E-21

1.10E-02 7.00E-02 1.30E-03 6.90E-10 1.30E-02 7.70E-02 2.50E-05 3.40E-15 3 2.70E+03 1.90E+03 1.10E-03 1.10E-25

1.60E-02 6.30E-02 9.90E-08 1.20E-21 2.10E-02 2.00E-02 8.80E-19 0 38 8.20E+02 1.40E-05 2.90E-33 0

1.5 0.9 7.6E-05 1.7E-12 1.6 0.26 3.5E-14

1.60E-02 8.60E-02 1.3 1.6 3.40E-03 1.90E-02 0.15 2.5 2.80E-02 6.00E-09 1.3 5.90E+02 5.50E+06 5.00E+03 1.00E-11

2.70E-02 0.24 2.8 7.80E-05 4.60E-16 5.60E-02 1.3 1.70E-02 6.40E-33 0 67 8.70E+05 1.60E+03 0 0

1.7 2.8 2.2 4.9E-05 1.3E-13 2.9 8.7 6.80E-03 2.3E-31

4.7E-21

2.7E-11

13 0.3 7.4E-09 2.70E-30

52 1500 2.9E-04

a 0 means the value is less than 1 × 10-50. b The C* values reported by John and Holder6 were incorrectly cited as R ) 4.3 Å for the structure I i tetradecahedral cell with Zc,i ) 21. The Ri value actually used to compute C* was apparently 4.182 Å, as reported by one reviewer. c All C* values given in this table used Ri and Zc,i parameters reported by John and Holder.6 d C*(trap) ) Langmuir constants calculated from LJD model using trapezoidal numerical integration. e C*(Gauss) ) Langmuir constants calculated from LJD model using Gauss-Legendre quadrature numerical integration based on the crystallographic structure of water clathrates in the first shell. f Q* ) Ratio of the two Langmuir constants C and C*.

Configurational Properties of Water Clathrates

J. Phys. Chem. B, Vol. 103, No. 30, 1999 6307 constant was calculated as a simple function of the number of water shells in the integration. Again, as explained in an earlier paper by Sparks and Tester,14 residual contributions beyond firstnearest-neighbor interactions can be quite large. In fact, in Figure 13, we observe a change of more than an order of magnitude in the value of the Langmuir constant as we go from a firstnearest-neighbor, single-shell treatment with 28 water molecules included to an asymptotically convergent value with more than 4 shells and greater than 1000 water molecules included. Conclusions

Figure 13. Study for subsequent water shell contributions using multidimensional quadrature integration for the structure II hexakaidecahedron cavity.

energy parameter, , to zero, we were able to equate the configurational partition function to the known integration volume. Note that in Table 4 the C* value for the LJD model reported by John and Holder6 and the C* values calculated in this study using the LJD parameters reported in ref 6 agree reasonably well except for the structure I large cavity for all values of σ, a, and /k examined. Given the specific mathematical form for 〈W(r)〉, one would expect identical values for C*. This discrepancy has been clarified by one of the reviewers, who noted that there was an error in the optimal cell radius, Ri, of 4.3 Å for the structure I large cavity reported in ref 6. The value actually used by John and Holder6 was Ri ) 4.182 Å as stated in two earlier papers by John and Holder.4,5 We believe that other discrepancies for the C* values reported in Table 4 are due to the slightly different oxygen position lattice coordinates used by each group for structure I and structure II hydrates (see John31 and Sparks23). As mentioned earlier, another internal check on the validity of our approach was provided by an independent calculation of 〈U/kT〉, the ensemble-averaged interaction potential which is given by: 〈U/kT〉 )

∫ ... ∫[U/kT] exp[-U/kT]r sin θ sin β dθ dφ dr dR dβ dγ ∫ ... ∫ exp[-U/kT]r sin θ sin β dθ dφ dr dR dβ dγ 2

2

(14) where again U ) U(r,θ,φ,R,β,γ). With the Metropolis Monte Carlo algorithm, an estimate of the average potential is given directly by

〈U/kT〉MC )

1

M

∑ Uj /kT M j

(15)

where M ) total number of configurations sampled (Tester et al.3). Close agreement was obtained between Monte Carlo results from eq 10 and multidimensional quadrature integrations of eq 9 for the structure I hydrate of cyclopropane (see Sparks23). In addition to the importance of including the host lattice asymmetries in the evaluation of the configurational partition function, the inclusion of the subsequent water-shell interactions is also important. A good example of this is illustrated in Figure 13. Using complete three-dimensional integrations, the Langmuir

This investigation quantitatively demonstrates the inaccuracies introduced by using the spherically symmetric Lennard-Jones Devonshire approximation to model water clathrates. For large guest molecules whose σ values are greater than 3.0 Å, significant errors in estimating the three-dimensional Langmuir constant (CJi) are observed. In addition, with guest-host asymmetries included, multiple water host shell interactions out to at least the fourth nearest neighbor are required to obtain consistently uniform values of CJi. Acknowledgment. The authors thank Norsk Hydro, the Idaho National Environmental Engineering Laboratory (INEEL), and the MIT Energy Laboratory for their partial support of this research. We also gratefully acknowledge insightful discussions with Prof. Irwin Oppenheim at MIT regarding the statistical mechanics of these systems. Bonnie Caputo, Gillian Kiley, and Loretta May did very commendable work in successfully producing countless electronic versions of the text, equations, and figures. References and Notes (1) Sloan, E. D., Jr. Clathrate Hydrates of Natural Gases; Marcel Dekker: New York and Basel, 1990. (2) Makogon, Y. F. Hydrates of Hydrocarbons; Pennwell Publishing Co.: Tulsa, OK, 1997. (3) Tester, J. W.; Bivins, R. L.; Herrick, C. C. AIChE J. 1972, 18, 1220. (4) John, V. T.; Holder, G. D. J. Phys. Chem. 1981, 85, 1811. (5) John, V. T.; Holder, G. D. J. Phys. Chem. 1982, 86, 455. (6) John, V. T.; Holder, G. D. J. Phys. Chem. 1985, 89, 3279. (7) Holder, G. D.; Hand, J. H. AIChE J. 1982, 28, 440. (8) Holder, G. D.; Malekar, S. T.; Sloan, E. D. Ind. Eng. Chem. Fundam. 1984, 23, 123. (9) Holder, G. D.; Hwang, M. J. Adsorption Theory in Clathrate Hydrates. In Fundamentals of Adsorption; Liapis, A. I., Ed.; Engineering Foundation: New York, 1987. (10) Rodger, P. M. J. Phys. Chem. 1989, 93, 6850. (11) Tse, J. S.; Davidson, D. W. Proc. 4th Can. Permafrost Conf. 1982, 329. (12) Hirschfelder, J. O.; Curtis, C. F.; Bird, R. B. Molecular Theory of Gases and Liquids; Wiley & Sons: New York, 1954; p 286. (13) Reid, R. C.; Prausnitz, J. M.; Poling, B. The Properties of Gases and Liquids, 4th ed.; McGraw- Hill: New York, 1987. (14) Sparks, K. A.; Tester, J. W. J. Phys. Chem. 1992, 96, 11022. (15) van der Waals, J. H. The Statistical Mechanics of Clathrate Compounds. Trans. Faraday Soc.: London 1956, 52, 184. (16) van der Waals, J. H.; Platteeuw, J. C. AdV. Chem. Phys. 1959, 2, 1. (17) Holder, G. D.; Corbin, G.; Papadopoulos, K. D. Ind. Eng. Chem. Fundam. 1980, 19, 282. (18) Holder, G. D.; John, V. T.; Yen, S. Geological Implications of Gas Production from In-situ Gas Hydrates; SPE/DOE Symposium on Unconventional Gas Recovery; SPE/DOE, 1980; p 8929. (19) Kihara, T. ReV. Mod. Phys. 1953, 25, 831. (20) Parrish, W. R.; Prausnitz, J. M. Ind. Eng. Chem. Process Des. DeV. 1972, 11, 26. (21) Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. J. Chem. Phys. 1953, 21, 1087.

6308 J. Phys. Chem. B, Vol. 103, No. 30, 1999 (22) Carnahan, B.; Luther, H. A.; Wilkes, J. O. Applied Numerical Methods; Wiley: New York, 1969. (23) Sparks, K. A. Configurational Properties of Water Clathrates Through Molecular Simulation. Ph.D. Thesis, Department of Chemical Engineering, Massachusetts Institute of Technology, Cambridge, MA, 1991. (24) von Stackelberg, M.; Muller, H. R. J. Chem. Phys. 1951, 19, 1319. (25) von Stackelberg, M.; Muller, H. R. Z. Elektrochem. 1954, 58, 25. (26) Jeffrey, G. A.; McMullan, R. K. Prog. Inorg. Chem. 1967, 8, 43. (27) Jeffrey, G. A. J. Inclusion Phenom. 1984, 1, 211. (28) Mak, T. C. W.; McMullan, R. K. J. Chem. Phys. 1965, 42, 2732.

Sparks et al. (29) Sloan, E. D., Jr.; Happel, J.; Hnatow, M. A. International Conference on Natural Gas Hydrates, The New York Academy of Sciences, New York, 1994; p 715. (30) McMullan, J. R. K.; Jeffrey, G. A. J. Chem. Phys. 1965, 42, 2725. (31) John, V. T. Compared Predictions of Gas Hydrate Phase Equilibria. Ph.D. Thesis, School of Engineering and Applied Science, Columbia University, New York, 1982. (32) Holder, G. D.; Zele, S.; Emile, R.; leBlond, C. Modeling Thermodynamics and Kinetics of Hydrate Formation. Ann. N.Y. Acad. Sci. 1994, 715, 344.