J . Phys. Chem. 1990, 94, 6080-6089
6080
Stabillty of Gas Hydrates P.Mark Rodger' British Gas plc, London Research Station, Michael Road, London S W6 2AD, U . K . (Received: January IO, 1989: In Final Form: November 22, 1989)
Molecular dynamics computer simulations have been used to examine the validity of the assumptions upon which the cell theory of gas hydrates is based. It was found that one of these approximations-that the free energy of the water crystal is unaffected by the presence of guest molecules in the interstices of its lattice structure-is inappropriate for gas hydrates. On the basis of our calculations it is conjectured that kinetic effects may be more important than thermodynamic considerations in predicting the conditions under which gas hydrates may be expected to form.
Unfortunately, the fact that the link between theory and experiment is itself approximate makes it difficult to assess the theory unambiguously: if calculated and observed properties do not agree, then it might be either because the theory is invalid or because the cavity potential used was inadequate. Thus, by themselves, experimental results are insufficient to determine whether or not the cell theory is the appropriate way of describing the thermodynamic properties of gas hydrates. One practical consequence of this is that previous attempts to improve the cell theory have not addressed themselves to the fundamental assumptions but rather to better ways of describing the cavity p ~ t e n t i a l , l *and -~~ as such, they have not determined whether the model is a good one, or whether it is merely providing a flexible empirical correlation for the observed behavior. A convenient way to resolve this question is with the aid of computer simulations. Since such calculations involve an exactly known model, there is no uncertainty in the description of the guest-host interactions, and so it is possible to examine the validity of assumptions a-c directly. Then, if the properties of the model system are similar to those of the real system (as is the case with simulations of gas one can expect any qualitative conclusions drawn from the simulations (such as the validity of these assumptions) to apply to real hydrates as well. We note that molecular dynamics (MD) and Monte Carlo (MC) computer simulation schemes are usually based on classical mechanics and
1. Introduction Gas hydrates (water clathrates of small gas molecules) form an interesting class of solids that continue to attract considerable attention, particularly with respect to their thermodynamic properties and the conditions under which they form. At present the only practical way of describing these properties is with the cell theory developed by van der Waals and Platteeuwl and since extended to cover multiphase systems.24 According to this model the water molecules form a well-defined crystal lattice (the host lattice) containing cavities into which small gas molecules (guests) may be adsorbed; under appropriate conditions the adsorption energy may then reduce the free energy of the hydrate sufficiently to make the hydrate phase more stable than either pure water or ice.5 However, while this theory has been able to describe the properties of many gas hydrate^,^,^* the parametrization is not always c o n s i ~ t e n t , ~ ?and + ~ there ~ are still other examples of gas hydrates for which the model does not appear to w ~ r k . ~ , ' ~In- ' ~ view of this it is useful to reexamine the assumptions and approximations on which the cell theory is based. Van der Waals and Platteeuw enumerated four basic assumptions in their development of the cell theory, and these have not been removed in subsequent extensions of the model. They may be stated as follows. (a) The free energy of the host lattice is independent of which molecules, if any, are occupying the cavities. Thus, for example, the contribution of the host lattice to the system's free energy must be the same when all the cavities in the lattice are empty as it is when all the cavities are occupied by a guest molecule. An important corollary to this is that the empty hydrate lattice must be at least metastable: if it were not then it would spontaneously rearrange to give a different structure, thus ensuring that it would sample a different region of phase space from the occupied host lattice and so would also have a different free energy. This hypothetical metastable empty hydrate is referred to as the phydrate, or the @phase of water. (b) The encaged molecules are confined within cavities that never contain more than one guest molecule. (c) Guest-guest interactions are negligible. (d) Quantum effects may be ignored. In other words, the host lattice is considered merely to give rise to an environment (often termed the cavity potential) in which the guest molecules evolve, and the thermodynamic properties of the hydrate result from the classical behavior of individual guest molecules within this cavity potential. In order to make the link between theory and experiment, one then needs a description of the cavity potential, which will necessarily be approximate; this is because the cavity potential arises from the intermolecular interactions between the guest and the host lattice, and these are not known exactly for real guest and host molecules. Further, the many-bodied nature of such interactions makes a simplified (and therefore approximate) description of the cavity potential desirable.
1972,I I. 26. (3) Ng, H.-J.; Robinson, D. B. AfChE J . 1977,23,477. (4) Menten, P. D.; Parrish, W. R.: Sloan, E. D. fnd. Eng. Chem. Process Des. Dev. 1981,20, 399. ( 5 ) Throughout this paper the term "hydrate" will be used to refer to clathrate hydrates rather than to hydrated crystalline salts. (6) Hollander, F.: Jeffrey, G. A. J . Chem. Phys. 1977,66, 4699. (7) Sloan, E. D.; Khoury, F. M.; Kobayashi, R. Ind. Eng. Chem. Fundam. 1976,IS, 318. ( 8 ) Dharmawardhana, P. B.; Parrish, W. R.; Sloan, E. D. fnd. Eng. Chem. Fundam. 1980,19,410. (9) Anderson, F. E.; Prausnitz, J. M. AfChE J . 1986,32, 1321 (10) Holder. G . D.: Hand. J. H. AfChE J . 1982,28. 440. ( I I ) Holder, G. D.; Godbole, S. AfChE J . 1982,28, 930. (I 2) Holder, G. D.: Corbin, G.; Papadopoulos, K. fnd. Eng. Chem. Fundam. 1980. 19, 282. (13) Cady, G. H. J . Phys. Chem. 1981,85,3225. (14) Ng,H.-J.; Robinson, D. B. Ind. Eng. Chem. Fundam. 1976,15,293. ( 1 5 ) McKoy, V.;Sinanoglu, 0. J . Chem. Phys. 1963,38,2946. (16) John, V. T.; Holder, G. D. J . Phys. Chem. 1981,85, 181 I . (17) John, V. T.; Holder, G. D. J . Phys. Chem. 1982,86,455. (18) John, V. T.; Holder, G. D. J . Phys. Chem. 1985,89,3279. (19) John, V. T.; Papadopoulos, K. D.; Holder, G. D. AfChE J . 1985,31, 252. (20) Tse, J. S . ; Klein, M. L.: McDonald, 1. R. J . Phys. Chem. 1983,87, 4198. (21) Tse, J . S.; Klein, M. L.; McDonald, I . R. J . Chem. Phys. 1983,78, 2096. (22) Tse. J. S . ; Klein, M. L.: McDonald, 1. R. J . Chem. Phys. 1984,81,
'Present address: Department of Chemistry, University of Reading, Reading RG6 2AD, UK.
(23) Marchi, M.; Mountain, R. D. J . Chem. Phys. 1987,86,6454. (24) Tse, J S.: McKinnon, W . R.; Marchi, M. J . Phys. Chem. 1987,91, 4188
0022-3654/90/2094-6080$02.50/0
( I ) van der Waals, J. H.; Platteeuw, J. C. Adu. Phys. Chem. 1959,2, I . (2) Parrish, W. R.: Prausnitz, J. M. Ind. Eng. Chem. Process Des. Deu.
6146.
0 I990 American Chemical Society
Stability of Gas Hydrates therefore can not be used to determine the applicability of assumption d; however the wide ranging success of MD and MC in reproducing the observed properties of molecular systems such as water and ice suggests that, to a very good approximation, quantum effects can be neglected. It is also worth mentioning that double occupancy of any cavity in a type I hydrate is exceedingly unlikely. Current estimates of molecular size25indicate that even helium would substantially fill the large cavities and so make double occupany very energetically unfavorable. In this paper we present the results of an MD study of the validity of the assumptions that underlie the cell theory. It is expected that this work, together with our recent investigation into the factors influencing the cavity potential,26 will provide a sound basis for future analytical developments in theoretical models of gas hydrates. A sketch of the simulation procedure will be presented in section 11. In section I11 the properties of an empty hydrate are compared with those of partially and completely filled methane hydrate in order to determine the effect of guest molecules on the host lattice and the importance of guest-guest interactions, while in section 1V we consider a system in which the guest molecules are purely repulsive soft-spheres (and hence have no stabilizing adsorption energy). Some of the ramifications of this work will be discussed in section V.
The Journal of Physical Chemistry, Vol. 94, No. 15, 1990 6081 TABLE I: Intermolecular Potential Parameters Lennard-Jones Interactions, V(r) = 4c[(u/r)l2 - ( ~ / r ) ~ ]
atom pair 0-0 CH,-0 CH4-H C H4-C H4
R-Rb
elk, K 78.18 91.71 50.8 1 164.2 164.2 164.2 164.2 164.2 164.2
Electrostatic Parameters for Water 0-H bond length = 1.0 A charge on 0 atom = -0.82 au charge on H atoms = +0.41 au LH-0-H = 1 0 9 O 28’ 4Cross parameters for soft-sphere interactions were obtained by maintaining the ratio Yg-x/Yg3 used for methane, where Y = u or c, X = 0 or H, and g denotes the guest molecule. Equivalent parameters based on Lorenz-Berthelot mixing rules (with uo = 3.16 A and uH = 3.0 A) are given in parentheses. bSoft-sphereguest-guest interactions. TABLE 11: Some Thermodynamic Properties of the Methane Hvdrates‘
11. Computational Details
a,b A
The simulation procedure is the same as that used in our study of the cavity potentialz6and so we will give just a brief outline of it here. We performed constant temperature-constant volume MD simulations of 46 water molecules and up to 8 guest molecules which were initially arranged to form a single unit cell of a type I gas hydrate; positions of the oxygen atoms were obtained from neutron scattering experiments6and the hydrogen atoms were then added in a manner that was consistent with the observed 2-fold disorder in the proton sites. Periodic boundary conditions were used to mimic an infinite crystal. We note that assumption a of the cell theory requires that the host lattice should not relax according to the nature or number of guests present, and so the use of isochoric conditions is appropriate in this study. The equations of motion were integrated with the velocityVerlet alg~rithm,~’ and rigid water molecule constraints were implemented with the RATTLE algorithm;z8these calculations used a timestep of 5 fs, which gave energy fluctuations of -0.01% of kT in constant energy simulations. Unless otherwise specified, the system was allowed to equilibrate for 10 ps before we accumulated a further 20 ps of hydrate dynamics for analysis. Intermolecular interactions were the same as those used in earlier studies,zo*z6 with water-water interactions being described by the simple point charge model (SPC)29and Lennard-Jones potentials modeling the guest-host and guest-guest interactions. As with our earlier work, we have used just a single site for the methane molecule and have included the guest-hydrogen interactions explicitly. Long range electrostatic effects were described with a Ewald summation, while Lennard-Jones interactions were set to zero for intermolecular separations greater than half the length of the unit cell. Soft-sphere guests were obtained by removing the attractive r4 term from the guest-host LennardJones potentials, and the potential parameters were then obtained from the guest values by maintaining the ratio ug-x/ug-sused for methanem (g denotes guest, while x is either hydrogen or oxygen); note that since these were used only in simulations of singly occupied hydrates (see section IV), it is possible to define an equivalent set of u parameters that obey the more commonly used Lorenz-Berthelot mixing rules. A summary of the intermolecular
12.03
(25) Wilson, J. W.; Heinbockel, J. H.; Outlaw, R. A. J. Chem. Phys. 1988, 89, 929. (26) Rodger, P. M. J . Phys. Chem. 1989, 93,6850. (27) Swope, W. C.; Anderson, H. C.; Berens, P. H.; Wilson, K. R. J . Chem. Phys. 1982, 76, 637. (28) Andersen, H. C. J . Comput. Phys. 1983, 52, 24. (29) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermens, J. In Inrermolecular Forces; Pullman, B., Ed.; Reidel: Dordrecht, Holland, 198 I .
u,o A
3.16 3.33 3.25 3.64 (3.50) 3.64 (3.50) 3.30 (2.88) 3.00 (2.33) 2.80 (1.96) 2.50 (1.41)
T, K 145
nt 0
260
8 0 1 8 0 1 0 1 0
270
0
200 220 240
1
1
8 11.88 11.92 11.94
200 240 270
0 0 0
1
-Uww,d P, MPa -388 (5) -457 (5) -285 (6) -270 (6) -286 (7) -256 (6) -236 (7) -193 (6) -160 (7) -118 (7) -140 ( 8 ) -127 (7) -95 (7) -80 (8) -7 (6) -18 (7) -13 (7) +20 (7)
kJ mol-’ 51.21 (2) 5 1 . 1 3 (2j 49.62 (3) 49.34 (3) 49.57 (4) 49.16 (4) 48.67 (4) 48.39 (4) 48.07 (4) 47.92 (5) 47.33 (4) 46.92 (9) 46.98 (5) 47.56 (8) 49.96 ( 3 ) 48.78 (4) 47.62 (5) 47.22 (5)
-uwg.d -uBo,d”
kJ mol-’
kJ mol-’
4 . 1 3( I )
0.165 ( I )
0.49 ( I ) 4.02 (1)
0.166 ( I )
0.49 (1) 0.49 ( I ) 0.47 ( I ) 0.48 ( I ) 3.86 (2)
0.166 ( I )
0.49 (1)
‘Numbers in parentheses are the statistical uncertainty in the final digit. b u n i t cell length of the hydrate crystal. CNumber of guest molecules in each unit cell. d~ denotes water and g denotes guest. Energies are expressed per mole of water molecules present in the appropriate h drate. CCalculatedby using all images of guest molecules within 12
1.
potential parameters is given in Table I. In order to assess the sensitivity of our conclusions to the model potentials used, we have also performed some simulations of methane hydrate by using the RWK intermolecular potential for ~ a t e r . ~ Although ~ , ~ ’ this led to some quantitative changes, the qualitative behavior was very similar to that of the SPC model and, in particular, resulted in the same conclusions about the validity of the cell theory. Since the RWK potential was computationally more expensive, we have concentrated on the SPC model and report only the latter calculations in this paper. 111. Methane Hydrates
From a comparison of the thermodynamic properties of hydrate lattices with different proportions of the cavities occupied (Le. different guest occupancies), it should be possible to gain some information about the effect of the guest on the host lattice and the significance of guest-guest interactions. Accordingly we have calculated the pressure and various contributions to the potential energy of the hydrate for a range of temperatures and guest occupancies at a fixed unit cell length of 12.03 A;32these are ~
_
_
_
______ _ ~
~
~~
~
~
(30) Reimers, J. R.; Watts, R. 0.;Klein, M. L. Chem. Phys. 1982, 64,95. (31) Coker, D. F.; Watts, R. 0. J . Phys. Chem. 1987, 91, 2513. (32) McMullan, R. K.; Jeffrey, G. A. J . Chem. Phys. 1965, 42, 2725.
6082 The Journal of Physical Chemistry, Vol. 94, No. 15, 1990
Rodger
reported in Table 11. In view of the long range spacial order present in gas hydrates, no attempt has been made to include the usual isotropic corrections for the spherical truncation of the Lennard-Jones potential^,^^ although a long-range correction for U, has been estimated by calculating the interactions with all guest images less than one unit cell length away (- 12 A); this special treatment was needed because the spherical cutoff distance of 6.015A was comparable with the shortest guest-guest distance. Three guest occupancies have been considered: the empty hydrate, a hydrate with one large cavity per unit cell occupied by a methane molecule, and fully occupied methane hydrate, which will be denoted HO,H 1, and HS,respectively. All energies have been expressed in units of kJ/mol of water molecules in order to facilitate a direct comparison between systems that differ in the number of guest molecules in the unit cell and also to emphasize the relative magnitudes of the various contributions to the energy. The negative pressures in Table 11 indicate that the length adopted for the unit cell is slightly too large for the model potentials used and so some further simulations have been performed with the cell length adjusted to ensure that the pressure of the empty hydrate was -0. We note, however, that all of the pressures are small with respect to the bulk modulus of ice (- IO4 MPa) and so are unlikely to cause any mechanical instabilities in the hydrate crystals. A useful guide to the effect of guest molecules on the host lattice is provided by the water-water contribution to the potential energy (Uww). If, as asserted in assumption a, the host lattice is unperturbed by the introduction of guest molecules into the cavities, then U,, should be independent of the extent to which these cavities are occupied. From Table I1 it can be seen that this is not the case, with U,, being significantly lower for H8 than it is for HI (0.23 and 0.58 kJ/mol at 200 and 270 K, respectively); note that a shift in the potential of -0.6 kJ/mol would change estimates of the ice/hydrate coexistence pressure at 2 7 0 K by about 30%," and so these differences are substantial. Thus, contrary to assumption a, the successive addition of guest molecules appears to stabilize the arrangement of water molecules in the host lattice. Although average energies for HO have been quoted in Table 11, there was a significant drift in the total energy during each of these trajectories, and so a direct comparison with HI and H8 is not meaningful. The extent of this drift is shown in Table 111, where we have listed the difference between the average energies calculated from successive 5-ps sections of the analysis portion of each trajectory; negative values indicate that the energy of the system decreases as the trajectory progresses. In comparison, the maximum cumulative drift for the HI and H8 calculations was only 0.1 kJ/mol-about the same as the statistical error in these energy differences-and the direction of the drift was not constant within a given trajectory. Thus the energy drift is unique to the HO system. Since all of the simulations used exactly the same program and conditions, the lack of any significant drift in the HI and H8 systems suggests that the drift in HO is not due to
numerical error, and this is supported by the direction of the drift-to lower energies (numerical errors usually result in an increasing energy). It is also unlikely that the drift results from the negative pressures, since it is just as significant in the simulations using a smaller unit cell length. We therefore conclude that this energy drift is a real effect and reflects the fact that the empty hydrate does not equilibrate on a 30-ps time scale. Indeed the empty hydrate must be unstable on a picosecond time scale, since it spontaneously rearranges to adopt a lower energy configuration than that of the /3-hydrate. This was found to be the case for the whole temperature range considered. If, as seems likely from the above discussion, the hydrate crystal structure collapses in the absence of guest molecules, then one would expect to see some accompanying changes in its intermolecular structure, and these should be reflected in the characteristic radial distribution functions (rdfs) for the host lattice. In fact, when one considers the various atom-atom rdf s describing the distribution of water molecules (gxu(r),X,Y E (O,H]), then one sees very little difference between the different systems. It appears that the structure is more clearly defined for H8,but this difference is very small, and HO and H1 are essentially identical. In retrospect this should not be surprising. The three water rdf s-goo, goH, and gHH-primarily reflect the short range interctions between water molecules, and these are dominated by the tetrahedral hydrogen bonding network found in all condensed phases of water. This tetrahedral coordination is particularly clear in the hydrates considered by the present study, with the first peak in goo containing 4.00 oxygen atoms under all the conditions considered, but it is just as evident in, for example, both ice ICand low-density amorphous ice; indeed the water rdfs for ice Iczoand low-density amorphous ice34are remarkably similar to those for the hydrate lattice for atom-atom separations up to -5 A, with all three showing a clear, double-shelled structure. Thus it is unlikely that these three functions will provide a sensitive test for the presence or absence of a specific type of water structure, viz. the hydrate crystal. To obtain information about changes in the host lattice structure it is better to consider properties that are unique to the gas hydrates, for example the presence of the cavities. Since the only physical forces maintaining the cavities in the empty hydrate must arise from the long range order of the water lattice, any change in this long range order is likely to alter the cavity structure, and so the cavities should be good indicators of structural changes in the hydrate lattice. Consequently, we have calculated the radial distribution of oxygen and hydrogen atoms about the center of the two types of cavity found in type I gas hydrates. If the host lattice is stable, then these should reveal a void at small r followed by a well-defined cavity wall at - 4 A; however, if the host lattice is unstable, then the cavity centers will loose their significance and these cavity rdfs will change, becoming completely structureless for amorphous or liquidlike arrangements of the water molecules. Representative examples of the cavity rdfs are depicted in Figures 2-5. In every case it was found that well-defined cavities persisted in the HI and H8 models, but that this structure was absent in HO. Thus, in the absence of any guest molecules, the water molecules readily migrate into the cavities in the @hydrate lattice, providing clear evidence that the &phase is indeed unstable on the picosecond time scale. The structure of the cavity wall also reveals that there are some small differences between the two occupied hydrates, with the cavity wall in H8 being narrower and more densely packed than that in H 1; it is probably this sharper definition of the host structure in H8 that is responsible for its lower value of U,,. Further evidence of the decay of the empty hydrate lattice can be found from the mean square displacement of the water molecules in the hydrate lattice as a function of time (MSD). For stable solids, the constituent molecules will vibrate around their well-defined lattice sites but will not diffuse away from these positions (at least not on a picosecond time scale); such behavior will be indicated by the appearance of a plateau in the MSD, with
(33) Singer, K.;Taylor, A.; Singer, J. V . L. Mol. Phys. 1977,33, 1757.
(34) Bosio, L.: Johari. G. P.: Teixeira, J. Phys. Reu. Left. 1986, 56, 460.
TABLE 111: Energy Drift during the HO Simul.tio& E2 - E l ,
E3 - E*,
E,
- EJ,
cumulative drift:
T, K
a,bA
kJ mol-'
kJ mol-'
kJ mol-'
kJ mol-'
145 200
12.03 12.03 11.88 12.03 12.03 1 1.92 12.03 12.03 1 I .94
-0.09 -0.08
-0.14 -0.13 -0.14 -0.20 -0.20 -0.26 -0.23 -0.38 -0.20
-0.18 -0.23 -0.26 -0.32 -0.25 -0.32 -0.42 -0.22 -0.42
-0.41 -0.44 -0.50 -0.58 -0.59 -0.77 -0.83 -2.34 -0.9 1
220 240 260 270
-0.10
-0.06 -0.14 -0.19 -0.18 -I .74 -0.29
E, is the average energy calculated during the time interval ( 5 ( i 1),5i)ps, taken from the analysis portion of each trajectory. b u n i t cell length of the hydrate crystal. CThisis equivalent to E4 - E , .
The Journal of Physical Chemistry, Vol. 94, No. 15. 1990 6083
Stability of Gas Hydrates
I
0
I 2.0
1.0
I 3.0
4.0
5.0
rtA
5r-7r 4
2t I \ I I
il-
\
I
2'o
riA 3'0
W
3
I
k 4.0 5.0 0 ''O
2.0 1.5
0
tc
I\
I \
1.0
2.0
r,A
Figure 2. Radial distribution of oxygen atoms around the center of the large cavities at temperatures of (a) 270 K and (b) 200 K: -*-, HO; -, HI;---, H8.
I
3.0
4.0
5.0
Figure 1. Atom-atom radial distribution functions for the H1 (-)
and
H8 (---) systems, with unit cell length a = 12.03 A and temperature T = 270 K: (a) oxygen-oxygen; (b) oxygen-hydrogen; (c) hydrogenhydrogen. The corresponding functions for HO were indistinguishable from those for H I .
the height of this plateau reflecting the amplitude of the vibrations about the lattice sites. Equlibrium MSDs for the oxygen atoms in the various simulated methane hydrates are depicted in Figure 6,and a plateau region is clearly visible for both the HI and H8 systems. For HO, however, the MSD reveals that even at the lowest temperature studied (145 K) the motion of the water molecules is diffusive rather than bound: there are no well-defined lattice sites in the empty hydrate, rather the water molecules rearrange themselves freely. Note that the plateau for H I is slightly higher than that for H8, indicating that the amplitude of lattice vibrations is greater in HI; this is consistent with the better definition of the cavity wall in H8 mentioned above. Mean square displacements of the oxygen atoms from their initial configuration (Le. from their positions at the beginning of the equilibration phase) have also been calculated, since these should indicate the time scale on which the host lattice begins to degrade. The curves for the HO system at various temperatures are shown in Figure 7, and it is clear that the lattice begins to decay very quickly: even for the coldest system considered (1 45 K) it takes only about 3 ps for the water molecules to begin t o diffuse away from their lattice sites, while at the other two temperatures there is clear evidence of degredation after just 1 ps. These results have some important consequences for the cell theory. If the structure of the host lattice depends on the number of guest molecules trapped within it, then it is clear that the free energy of the water lattice cannot be independent of whether or not the cavities are occupied. In particular, since the empty hydrate spontaneously rearranges to a lower energy configuration, the mechanical stability of a gas hydrate is not a property of the
Figure 3. Radial distribution of hydrogen atoms around the center of the large cavities at temperatures of (a) 270 K and (b) 200 K: HO; -, HI;---, H8. -e-,
host lattice, but rather arises from t t - way the guest molecules perturb the host lattice. This is contrary to assumption a and its corollary, and so means that the cell theory of gas hydrates requires some modification. In order to achieve such modification, it is first necessary to understand the mechanism by which the guests perturb the host lattice, and a clue to this is contained in the small differences that
Rodger
6084 The Journal of Physical Chemistry, Vol. 94, No. 15, 1990
/
/' /'
6-
/'
/'
g4-
2,A+--
4.q,JL . .I --.v.-.,./.I 1.o 2.0
0
./ / ..
riA
3.0
4.0
5.0
Figure 4. Radial distribution of oxygen atoms around the center of the small cavities at temperatures of (a) 270 K and (b) 200 K: -.-, HO
-, HI;---, H8.
t/P
Figure 6. Mean square displacement (in
-~
A2)of oxygen atoms in the
hydrate lattice, calculated from the analysis portion of the trajectories at temperatures of (a) 145 K and (b) 270 K: -.-, HO; -, HI; ---, H8.
7
6 5-
-
c4cn
321
-
r'pi^s'h?h^V-.C-/"
,A'#'1
,
1
Y'
I U'
1
¶
8
t/P
Figure 5. Radial distribution of hydrogen atoms around the center of the small cavities at temperatures of (a) 270 K and (b) 200 K: HO -, HI;---, H8.
Figure 7. Mean square displacement (in A2)of oxygen atoms from their starting configuration in the HO simulations: (a) T = 145 K; (b) ---, T = 200 K;-, T = 270 K.
do exist between the singly and fully occupied hydrates (HI and H8).We noted above that in €38 the water molecules in the walls of the cavities are more tighly bound to their lattice sites and so give rise to a more sharply defined hydrate crystal structure than in H 1. Since the only difference between these two systems is that more of the cavities in H8 are occupied, this suggests that the guest molecules tend to confine the water molecules and,
thereby, restrict their ability to move. If this is the case, then the mechanism by which the guest molecules stabilize the hydrate structure is not via a favorable adsorption energy, but by being large enough to get in the way of the rearranging water molecules; i.e. they provide an excluded volume that damps out those critical modes of lattice vibration that lead to rearrangement of the host lattice. The instability of the empty hydrate-and the freedom
-e-,
Stability of Gas Hydrates with which the water molecules in the empty hydrate can migrate (see Figure 6)-would then be a necessary consequence of the lack of guest molecules: in the absence of any excluded volume there is nothing to prevent the water molecules from entering the cavities and so obtaining enough room in which to move around each other. In view of these conclusions, the sharp distinction between the properties of the empty and occupied hydrates is interesting, and it conveys some important information about the critical lattice vibrations. If the particular motion of the water molecules that leads to degradation of the host lattice were to involve deformations of just one or two cavities, then one would expect to see degradation of the hydrate crystal whenever a significant proportion of the cavities were unoccupied. However, despite the relatively large free volume in the H I system, with 88% of the cavities being unoccupied, all of its cavities remain well defined and have walls that are only slightly broader than those in H8;this indicates that in H1 the motion of the water molecules is still strongly restricted. Thus the critical lattice modes must require a large-scale cooperative motion of the water molecules involving the simultaneous deformation of at least several of the cavities. As such, they are likely to be low-frequency modes and may therefore interact strongly with the low-frequency translational motion of the guest molecules within the cavities. We note that such large scale cooperative motions are not favored by the small simulation cell size used in this work, and so if anything, our simulations are likely to underestimate the extent of degradation of the hydrate crystal (and hence, may also underestimate the significance of the breakdown of assumption a). The simulations described in this section may also be used to test assumptions b and c, and in this case the calculations support their validity. It is quite apparent from Table I1 that the guest-guest interactions are negligible. Even when the hydrate is fully occupied, and hence the import of the guest-guest term is maximized, it still represents less than 0.5% of the total energy; this figure is about the same for all conditions studied. Also, the guest molecules were found to be trapped within the cavities, with no evidence of cavity-to-cavity hopping. Calculations of the distribution of guest molecules around the center of the cavities revealed that each guest molecule was confined to a sphere of radius 1 A around the center of a cavity, which is much smaller than the shortest cavity center-cavity center separation of -6
-
A.
To conclude this section it is interesting to speculate about why the decay of the 0-hydrate lattice has not been observed in previous simulation^.^^-^^*^^ Certainly one reason is that the present work is the only one to have examined either the cavity radial distribution functions or the mean square displacement of atoms in the host lattice, and these provide the clearest evidence that the 0hydrate is unstable. However, the energy has been considered in other simulations, and one must look more closely to explain why other workers have not found evidence of the energy drift we observed in the HO system. A possible explanation for this may lie in the different ensembles probed by the various simulations: all the previous simulations have kept the total energy constant (coupled either with constant volume2*22 or constant p r e s s ~ r e ~whereas ~ , ~ ~ ) we have used constant temperature (i.e. constant kinetic energy) conditions. Since the constant temperature simulations are closely related to simulated annealing, which is a common technique for locating minima in multidimensional surfaces, it is quite likely that constant temperature techniques will be more efficient than the constant energy methods at escaping from unstable or metastable regions of the potential energy surface, and so will reveal any rearrangement of the host lattice more clearly. It is also possible that differences in statistical noise might make a significant drift toward higher temperatures-which is the constant energy analogue of our drift to lower energies-harder to detect than the energy drift. IV. Soft-Sphere Hydrates If our conjecture that bulk gas hydrates are stabilized by excluded volume effects is correct, then hydrate structural properties
The Journal of Physical Chemistry, Vol. 94, No. 15, 1990 6085 TABLE I V Properties of the Soft-Sphere Hydrates at 270 K O a,b A
12.03
11.94
-L, +L., guest HI ~ 3 6 R33 R30 R28 R25 HlA R33A
f, MPa
kJ mol-l
-95 -11 -93 -65 -74 -90 +20 -66
46.98 46.98 46.94 46.83 46.89 46.84 47.22 47.20
(7) i7j (8) (7) (9) (7) (7) (7)
(5) (5j (5)
(5) (5) (6) (5) (5)
kJ md-l
-0.48 (1) 0.267’(i) 0.134 (1) 0.065 (3) 0.042 (4) 0.042 ( 5 ) -0.49 (1) 0.142 (2)
AA,b
kJ mol-l -0.156 0.642 0.428 0.330
(4) (8j (2) (6)
0.261 (17) -0.155 (3) 0.447 (2)
a Numbers in parentheses are the statistical uncertainty in the final digit(s). All energies are expressed per mole of water molecules in the appropriate hydrate. Helmholtz free energy relative to a hypothetical hydrate water lattice (see text).
should not depend on the nature of any attractive forces between guest and host molecules, and in particular, the structure should be the same for a purely repulsive guest as it is for a LennardJones guest. (We stress that this conjecture does not relate to how the hydrate forms, but to how stable it is once it has formed; growth of a gas hydrate and the structure resulting from this growth will be discussed in section V.) Further, it should be possible to gain some idea of how much excluded volume is needed to stabilize the host lattice by varying the radius of the repulsive guest molecule. In order to test this we have performed a number of simulations in which a single repulsive soft-sphere was placed in one of the large cavities of an otherwise empty unit cell of a type I gas hydrate, with the radius of the guest particles ranging from 55% to 80% of the cavity radius; these simulations will be denoted Rnm, where the u parameter for the soft-sphere is approximately n.m A (see Table I). As with the HI and H8 calculations, we did not observe any significant or consistent drift of the energy in any of the trajectories, indicating that the softsphere hydrates may indeed be stable on a picosecond time scale. Table IV lists some of the thermodynamic properties of the repulsive soft-sphere hydrates. We particularly note that there is no significant difference between the various host-host contributions to the potential energy, U,,, and that these agree with the corresponding values for the HI system to within statistical uncertainty. This suggests that although the behavior of the water lattice depends on the number of guest molecules present, it is quite insensitive to the exact nature of the guest.3s Confirmation of this insensitivity can be obtained from the various radial distribution functions arising in the soft-sphere hydrates. As would be expected from the discussion in section 111, goo,&H, and gHH were indistinguishable for the various singly occupied hydrates. More surprising, however, was the fact that the structure of the large cavity walls revealed very little difference between H1 and the various R systems, as can be seen from Figure 8. Indeed the only significant difference was found in the distribution of water molecules around the small cavities, which is depicted in Figure 9. In this case the wall of the cavity does appear to spread out more for the smaller guest molecules, but the effect is still minor and is clearly apparent only for the smallest guest (R25). Hence we conclude that the attractive guest-host interactions are not essential to maintaining the structure of bulk gas hydrates and that it is the excluded volume of the guest that prevents rearrangement of the host lattice. The soft-sphere simulations also indicate that the excluded volume required to stabilize the host lattice is minimal, with a single small guest molecule restricting the motion of the water molecules almost as much as a guest with a radius that is 50% larger. Thus it appears that the large scale cooperative lattice modes that are responsible for the decay of the hydrate lattice must involve large amplitude deformations of the cavities, so that a small guest will disrupt this motion just as effectively as a large (35) Greater sensitivity to the nature of the guest might be introduced by allowing the host lattice to relax about the guest molecule, e.g. by changing the unit cell length; however, such effects are beyond the scope of the present (isochoric) study, and in any case would only represent a further breakdown in assumption a.
6086
The Journal of Physical Chemistry, Vol. 94, No. 15, I990
r/A
Radial distribution of (a) oxygen and (b) hydrogen atoms around the center of the large cavities in various singly occupied hydrates: -, HI; R36; * * * , R30; -*-, R25. Figure 8.
---1
Rodger barriers. From our simulations it is not possible to say more than that the barrier for singly occupied systems is sufficient to prevent passage on a time scale of -20 ps; however, from the nature of the interactions that give rise to it (i.e., substantial overlap between two or more molecules), it is likely that the time scale for crossing the potential barrier is macroscopically long, so that gas hydrates may be kinetically stable even when they are thermodynamically unstable. In principle, the question of kinetic or thermodynamic stability could be addressed by comparing the free energies of our simulated systems with that of ice; unfortunately, because the rearrangement of the water molecules in these two systems is so different, such a calculation is likely to be a major undertaking and is beyond the scope of this work. An alternative is to use the @hydrate as a reference system instead of ice, as is done in the cell theory; the free energy differences between the various hydrate systems can then be compared to see whether stability of the host lattice is consistent with thermodynamic stability. Since the empty hydrate is unstable, it is not possible to use the HO simulation for this purpose, and so we have defined a hypothetical @-hydratethat is the water lattice obtained by removing the guest molecule from the singly occupied hydrates after the simulation has been performed; i.e. the 8-hydrate is taken to be the host lattice as seen in the HI and R simulations. In view of the similarity of the host lattice properties in the H1 and R systems this definition is reasonable-it is also consistent with the motivation behind assumption a of the cell theory and so provides a self-consistent test of the cell theory. We note, however, that this definition would not be so useful for comparing hydrates with different guest occupancies, as their lattice structures would differ (see, for example, the differences between H I and H8 noted in section 111). To calculate the change in free energy induced in the system by varying the guest molecule, we have used the equation A2
6t
r/A
Radial distribution of (a) oxygen and (b) hydrogen atoms around the center of the small cavities in various singly occupied hydrates: -, H1; ---, R36; R30; ---,R25. Figure 9.
..e,
guest. Conversely, this implies that the energy barrier to rearrangement in an occupied hydrate will be very high, since it requires a substantial overlap of the electron distributions of two or more molecules. Furthermore, since these critical vibrations involve the simultaneous entry of water molecules into several cavities, larger guest occupancies than that of H1 (such as those observed e ~ p e r i m e n t a l l y ' ~are ~ ~ likely ~ * ~ ~to ) lead to even higher
- A I = -kT In ((expt(U2 - Ul)/(kT)l)J
(1)
and allowed the reference system, 1, to range over all the singly occupied simulations. Averages and standard deviations from calculations with the different reference systems are reported in Table IV; the small magnitude of the standard deviations indicates that the different reference systems gave very similar results. In eq 1, k is Boltzmann's constant, T is the temperature, indicates and average over the equilibrium ensemble of system 1, and Ai and U, are the Helmholz free energy and potential energy of system i, respectively. Since the water-water interactions were the same in all the systems studied, the difference in potential energies is just the difference between the two guest-host energies, UwB.The free energy of trapping a single guest in the @-hydrate reference described above was obtained by letting Uwglbe identically zero in eq 1. From Table IV it is clear that all of the R systems are thermodynamically less stable than the 8-hydrate which, in turn, has already been shown to be unstable; because of the purely repulsive nature of the guest-host interactions in the R systems this result was expected, but it does confirm that the stability observed in the R simulations arises from kinetic rather than thermodynamic considerations. A more instructive use for the calculated free energies is to compare them with empirical estimates of the free energy difference between the &phase and ice. These two quantities are not directly comparable, since the calculations give the Helmholtz free energy whereas the experimental estimates are for Gibb's free energy. However, they are linked through the following thermodynamic identities: ~r=
( a G / d N i ) , , N , = (aA/aNi)v,T,N,
G = CpjNj
(2) (3)
where N k ( k = i, j ) is the number of guest or host molecules, V is the volume, G is the Gibb's free energy, and p is the chemical (36) Barrer, R. M.; Edge, A. V. J . Proc. R. Soc. London l967,300A, 1. (37) Cady, G. H. J . Fluorine Chem. 1978, I I , 2 2 5 .
Stability of Gas Hydrates potential. Note that eq 1 evaluates the change in A when one guest molecule is added to the hydrate at constant V,T, and N-; this is the same as the chemical potential of a guest molecule in a large cavity, and so may be used with eq 3 to give the change in Gibb’s free energy associated with encaging the guest molecules. Thus it can be seen from Table IV that the chemical potential of methane trapped in one of the large cavhies is quite small, being only -0.156 kJ/mol. We have also calculated the chemical potential of methane in one of the small cavities in the hypothetical fl-hydrate, this time using Widom’s test particle insertion method,38939and found it to be -0.091 f 0.002 kJ/mol. If, to be consistent with the cell theory, it is assumed that the chemical potential is independent of whether or not the other cavities are occupied, then the free energy difference between the fully occupied methane hydrate and the &hydrate (which will be denoted AGJ is found from eq 3 to be -1.12 f 0.01 kJ/mol of water. The magnitude of this value is actually smaller than the current empirical estimates of the free energy difference between the phydrate and ice at 273 K (AGO,), which are in the range 1.24-1.30 kJ/mo1,2.8J2*40 and so suggests that methane hydrate is thermodynamically less stable than ice at temperatures just below the ice point; we note that methane hydrate is known to form at these temperatures. There are several approximations in this argument, however, that may make the resultant sign of the free energy difference between methane hydrate and ice uncertain. In the first place, AGO, values have been obtained by empirically fitting the parameters in the cell theory to experimental data, and since we have already shown that one of the basic assumptions in the cell theory is not valid for gas hydrates, the physical interpretation of these empirically fitted parameters is not obvious. Similarly, the fact that the guest molecules alter the behavior of the host lattice also means that AGf will not be linear with respect to the number of guest molecules present. Thus both AGO, and AGf are likely to differ from the values quoted above, and in view of the size of the estimated difference in free energy between ice and methane hydrate (only 0.1-0.2 kJ/mol), this could easily change the conclusions. Further, the potentials used in this work are model potentials and hence do not give an exact description of methane hydrate; as a result, it is again unlikely that the magnitude of AGO, AGf is large enough to be conclusive. Finally, it should be remembered that the values in Table IV have not been corrected for the truncation of the Lennard-Jones potentials beyond a center of mass separation equal to half the unit cell length. Estimates of this truncation correction based on a spacially isotropic distribution of molecules beyond the cutoff distance33 suggest that this could increase lAGA by as much as SO%, which is enough to ensure that methane hydrate is stable; we note, however, that for systems like gas hydrates, which have long range spacial order and large voids in the host lattice, the isotropic approximation is not likely to be very accurate. This need to consider the truncation correction confirms the findings of John and Holder” about the importance of long-range guest-host interactions in determining the conditions under which a hydrate may form and suggests that the long range interactions may even be decisive in determining whether the hydrate will form a t all. One conclusion that can be drawn from the above discussion is that the difference between stability and instability in a gas hydrate can be very subtle, since the answer may depend on the cumulative effect of several small factors. As a result, it is essential that the theory give a good description of these subtleties if it is to have much predictive power.
+
V. Discussion From the results of our MD simulations it is apparent that one of the fundamental assumptions of the cell theory-that the free energy of the host lattice is unperturbed by the presence of guest (38) Widom, B. J . Chem. Phys. 1963.39, 2808. (39) Romano, S.; Singer, K. Mol. f h y s . 1979, 37, 1765. (40) Davidson, D. W.; Handa, Y.P.; Ripmeester, J. A. J . f h y s . Chem.
1986, 90,6549.
The Journal of Physical Chemistry, Vol. 94, No. 15, 1990 6087 molecules-is inappropriate for gas hydrates. Although the properties of the host lattice are remarkably insensitive to the exact nature of the guest molecule, they do depend appreciably on the number of guest molecules present. This is particularly obvious when occupied and unoccupied hydrates are compared, since the host lattice is mechanically stable in the former but spontaneously rearranges to a lower energy configuration in the latter (possibly amorphous ice, although this cannot be confirmed as the simulations were stopped before the rearrangement was complete); however, there are also significant differences between the various occupied hydrates, with the water molecules in singly occupied hydrates (Le. with one guest molecule per unit cell) having a higher potential energy and being more loosely bound to their lattice sites than those of the fully occupied hydrates. This means that the phonon modes of the host lattice must change in the presence of guest molecules, and hence so too will the host contribution to the free energy of the system. Thus there are some basic difficulties in the cell theory of gas hydrates as it now stands. One way to circumvent these difficulties might be to adopt some other reference system than the &hydrate, with an obvious choice being a hydrate with the appropriate proportion of its cavities occupied. For example, the host properties of all the singly occupied hydrates we studied were very similar, so that one of these might be a suitable reference system for calculating the properties of the others; indeed we found that each of the singly occupied hydrates did provide a good starting point from which to evaluate the free energy of each of the other singly occupied systems. This solution is not completely satifactory, however, as the host lattice is dependent on the percentage of cavities occupied, and so one would need to adopt a different reference system each time the guest occupancy changed. In practice, this means that the properties defining the reference system, namely the variation in its chemical potential and enthalpy from those of ice, would need to become functions of guest occupancy; since these properties are determined empirically, to adopt such a continuously variable reference system runs the risk of becoming merely a more flexible (and possibly more arbitrary) means of fitting experimental data. Another drawback with this method is that it will increase the complexity involved in calculations of hydrate properties: since properties such as hydrate pressures (i.e. the pressures at which the hydrates will just begin to form) depend on the percentage of cavities occupied and this percentage depends on the reference chemical potential which, in its turn, depends on the percentage of cavities occupied, it is probable that some sort of iterative cycle would have to be introduced into the calculation procedure. As a result of these factors, it is likely that a change in the reference system used with the cell theory will neither improve our understanding of the formation and properties of gas hydrates nor increase our confidence in the predictive power of the cell theory for gas hydrates. An alternative conclusion to be drawn from the inapplicability of assumption a to gas hydrate systems is that the cell theory is on the wrong track. This theory embodies a picture of gas hydrates in which a metastable (and hence mechanically stable) crystalline water lattice is made thermodynamically stable by the favorable effect it has on the guest molecules trapped within it. Thus their thermodynamic stability would arise not from any distortion, restriction, or relaxation that the guest molecules induce in the host lattice but from the way the attractive dispersion forces between the guest and water molecules lower the energy of the guest. This picture contradicts the results of our simulations in at least two ways: (i) on its own the host lattice is unstable, so that the mechanical stability of gas hydrates is not a property of the water lattice but is rather a consequence of the way the guest molecules restrict the motion of the host molecules, and (ii) the attractive guest-host forces do not appear to be essential to stabilizing gas hydrates in the bulk, with purely repulsive guests giving rise to precisely the same host behavior as a model methane guest. What our simulations do reveal is that hard core repulsive guest-host interactions play a central role in stabilizing gas hydrates. There are several observations that support this conclusion: (i) the host lattice collapses unless a guest is physically occupying
6088 The Journal of Physical Chemistry, Vol. 94, No. 15, 1990
at least one cavity in each unit cell, (ii) purely repulsive guests are just as effective a t stabilizing the host lattice as guests that may be attracted to the water molecules, and (iii) any small differences that do exist between the various singly occupied hydrates appear to be related to the size of the guest, with the smallest guest (R25) giving the water molecules the most room in which to move and so inducing the least structure in the host lattice. These factors have led us to conjecture that the role of the guest molecule is merely to prevent the water molecules from having enough room in which to rearrange; in other words, the guest molecules constitute an excluded volume which disrupts any critical lattice vibrations and thereby prevents rearrangement of the host lattice. In this case the characteristic crystal structures of the host lattice arise as a means of maximizing hydrogen bond formation around these excluded volumes. This interpretation is actually very similar to the current understanding of the anomalous behavior of supercooled water. A number of the thermodynamic properties of supercooled water appear to diverge at a temperature of about -45 O C , and this is usually considered to arise from an association between water molecules, probably due to hydrogen b 0 n d i n g . 4 ~ ~The ~ presence of a polar or ionic solute tends to suppress these anomalies (much as methanol inhibits hydrate formation9), but the introduction of hydrate-forming solutes appears to enhance them.z5-44.45There now appears to be a consensus that this behavior is due to the spontaneous formation of transient clathrate-like cages in water at low (especially supercooled) temperature^;^^^^^^ these cages give a reasonably favorable arrangement of hydrogen bonds but, being transient, must be mechanically unstable-in this respect they recall the behavior we observed in the @-hydrate. If polar or ionic solutes are present, then they will interact strongly with the water molecules and cause the water molecules to associate with the solute rather than with other water molecules; hence they will inhibit the formation of clathrate-like cages in water. On the other hand, hydrate-forming molecules do not interact strongly with water and so do not disrupt the hydrogen bonding networks set up in supercooled water, they can, however, occupy the cages and so inhibit their decay. Indeed, dissolving hydrate-forming substances in water actively promotes the formation of clathrate-like cages49*50 as a means of maximizing hydrogen bonding interactions around a body to which the water molecules are not strongly attracted. Thus the water structure is again determined by the interplay between the excluded volume of the guest and the formation of hydrogen bonds between the water molecules. Although excluded volume effects are sufficient to prevent the decay of gas hydrates in the bulk, they do not suffice to explain how the hydrates form in the first place. Clearly, if there were no attractive forces between the guest and host molecules, then they would not associate, and the hydrate would never have the opportunity to form; thus the attractive guest-host interactions must play a major role in the initial formation and growth of gas hydrates. At first sight this conclusion might appear contrary to the results of our simulations, but this is not the case. It is well-known that the formation and growth of gas hydrate crystals is primarily an interfacial p h e n o m e n ~ n , 'and ~ ~ ~as~such, the factors that influence formation and growth will differ from those affecting the bulk properties studied in our simulations. For example, many of the "cavities" on the surface of a hydrate crystal will not be fully formed and hence will not completely encage a guest molecule. (41) Stillinger, F. H. Science 1980, 209, 451. (42) Rice, S.A.; Sceats, M. G. J . Phys. Chem. 1981, 85, 1108. (43) Speedy,R. J. J . Phys. Chem. 1984,88, 3364. (44) Speedy, R. J.; Ballance, J. A,; Cornish, B. D. J . Phys. Chem. 1983, 87, 325. (45) Sorensen, C. M. J . Chem. Phys. 1983, 79, 1455. (46) Oguni, M.; Angell. C. A. J . Phys. Chem. 1983, 87, 1848. (47) Euliss, G. W.; Sorensen, C. M. J . Chem. Phys. 1984, 80. 4767. (48) D'Arrigo, G.; Paparelli, A. J . Chem. Phys. 1988, 88, 405. (49) Owicki, J. C.; Scheraga, H. A. J . Am. Chem. Sot. 1977, 99, 7413. (50) Swaminathan, S.; Harrison, S.W.; Beveridge. D. L. J . Am. Chem. SOC.1978, 100, 5705. ( 5 I ) Barrer, R. M.; Ruzicka, D. J. Trans. Furuday Sot. 1962, 58. 2262.
Rodger As a result, the host lattice will not trap guest molecules on its surface, and any association between guest molecules and the hydrate surface must arise from the attractive guest-host forces. Thus energetic considerations will be more important on the surface than they are in the bulk. In addition, since hydrate formation usually occurs at either a solid/gas or liquid/gas interface, there will be diffusion of the constituent molecules across the interface and this will tend to set up an equilibrium distribution of molecules in the interfacial region. Both of these factors mean that migration of water and guest molecules can occur at the hydrate surface and that this will be governed primarily by energetic considerations; hence the surface of a gas hydrate will be stable only if it is in equilibrium with its surroundings. We note that the possible importance of the guest in determining hydrate formation was first suggested more than 20 years ago,36but the idea does not seem to have been developed since then. The picture of gas hydrates that begins to emerge from this discussion is one in which a kinetically stable bulk is encased within a thermodynamically stable surface. Guest molecules adsorb to the water surface (be that solid or liquid) because it is energetically favorable for them to do so. If a few water molecules are then added to this surface, it is likely that they will simply orient themselves around the adsorbed molecules, particularly if the subsequent arrangement is favorable to the formation of hydrogen bonds among the water molecules. Eventually, when enough water molecules have been added to the surface, it may become energetically more favorable for the guest molecules to desorb and the water molecules to reorient themselves; however there will be an energy barrier to such a reaction, and if this is sufficiently high, then the rearrangement is unlikely to occur and the guest-water mixture will appear to be stable. Subsequent adsorption of guests and condensation of water will then continue this process and so give rise to a hydrate crystal. Thus the formation of a gas hydrate may be viewed as the growth of a lattice of water molecules on a contaminated water or ice surface. Since it is the contamination of the surface that disrupts the packing of the water molecules, it is likely that the height of this energy barrier will be related to the extent of the contamination, i.e. to the coverage of the surface by the adsorbent; this coverage is governed by the experimental conditions. If the coverage is very low, then the adsorbed molecules will not create enough gaps in the water lattice to force the formation of a hydrate crystal. On the other hand, if the coverage is extensive, then the guest molecules will reinforce their disruption of the surface water structure and so will lead to the open hydrogen bonding network that characterizes gas hydrates. Between these two extremes there will be a critical adsorption coverage that just ensures the formation of a hydrate water lattice on a macroscopic scale. The presure at which this coverage occurs will be the hydrate pressure. Another consequence of this model is that the type of hydrate formed should depend largely upon the adsorption structure at the water surface and be relatively insensitive to the guest-cavity interactions. Clearly a given type of hydrate will not form if the guests are too big to fit into its cavities, but beyond that it is likely to be the pattern of adsorption sites on the surface that determines which of the possible types of hydrate will form. One would therefore expect the pattern of adsorption sites occupied by type I hydrate-forming molecules to be similar to each other but to be distinct from the adsorption structure adopted by type I1 hydrate-forming molecules. Note that the relative occupancy of the various adsorption sites (which determines the adsorption structure) will depend on both the size of the guest relative to the adsorption site and the strength of the guest-host dispersion forces, so that the hydrate structure will arise from an interplay of dispersion forces and molecular size, but this will be determined at the surface rather than in the bulk. Although this model is not yet quantitative-such an implementation is currently being developed-it does explain some of the experimental and theoretical difficulties that have been encountered with gas hydrates. For example, if the bulk hydrate is kinetically stable, then it is likely that the distribution of guest molecules within it will be inhomogeneous. The resulting dis-
Stability of Gas Hydrates tribution will partly reflect the history of the conditions under which the hydrate formed, but may also depend upon the kinetics of formation, since if the hydrate grows sufficiently rapidly, then an equilibrium adsorption of guest molecules will not be maintained during hydrate growth. Since these factors may depend on the experimental conditions and apparatus used, it should not be surprising that experimental hydration numbers for gas hydrates show considerable ~ariati0ns.l~ From a theoretical point of view, it also explains why the cell theory is often able to reproduce the observed behavior, as the fundamental process in both models is adsorption. However, the cell theory does use the wrong model for the environment of the guest molecule when empirically fitting the intermolecular potentials-it considers a spherical cavity in the bulk instead of a partially formed cavity on the surface-and so one should expect it to yield potential parameters that differ from those obtained from, e.g., viscosity measurements (compare refs 9 and 52). Further, the variaton of the potential parameters with gas composition and experimental conditions (compare refs 9, 3, and 10) can be seen to arise from variations in the tendency of different molecules to adsorb onto different members of the ensemble of partially formed cavities found on the hydrate surface. We also note that, unlike the cell theory, our model is in line with chemical intuition about hydrophobic and hydrophilic interactions. These concepts express the conventional wisdom that water molecules, being polar, will be strong attracted to each other, and that this attraction will exclude other molecules unless those other molecules attract the water even more strongly, i.e. unless they themselves are either charged or very polar. While these considerations admit the possibility that nonpolar molecules may adsorb onto the water surface and then become kinetically trapped within a water lattice, they do suggest that the immersion of nonpolar molecules in a water lattice (which is fundamental to the cell theory) should not be thermodynamically more stable than separate water and guest phases. In view of the present lack of information about the surface of gas hydrates, however, the mechanism for hydrate formation that we have postulated in this section can only be conjecture. More work is certainly required to clarify the role of the surface in stabilizing gas hydrates, and both experimental and theoretical studies, especially of the possible relationship between hydrate formation and the nature of the surface adsorption, would be (52) Tee, L. S.;Gotoh, S.;Stewart, W. E. Ind. Eng. Chem. Fundam. 1966, 5, 363.
The Journal of Physical Chemistry, Vol. 94, No. 15, 1990 6089 helpful. The recent interest in quantitative kinetic studies of gas hydrate f ~ r m a t i o n ~is ~also - ~ ~welcome in this context.
VI. Conclusion During the course of this paper we have presented the results of a molecular dynamics computer simulation study of the factors influencing the stability and thermodynamic properties of gas hydrates. These results indicate that there are some fundamental difficulties with the conventional picture of gas hydrates (embodied in the cell theory), in which a metastable water lattice is made thermodynamically stable (but is mechanically unperturbed) by a favourable interaction energy with the guest molecules it encages. In fact, it is clear from our simulations that the water lattice is unstable rather than metastable, and as such, the guest molecules must serve to damp out those critical lattice vibrations that lead to rearrangement of the host matrix. Further, since repulsive soft-sphere hydrates were found to exhibit the same host structure as methane hydrate, it appears that it is the hard core repulsions and not the attractive dispersion forces that are important. Thus the excluded volume of the guest molecules seems to mechanically alter the behavior of the host lattice, and this perturbation plays an important part in stabilizing gas hydrates in the bulk. On the basis of our calculations we have conjectured a new model to explain the formation and stability of gas hydrates that includes these excluded volume effects. In this model the attractive dispersion forces cause a layer of guest molecules to adsorb to a water or ice surface; this layer of adsorbed "contaminants" then disrupts the structure adopted by any water molecules that are subsequently added onto this surface, causing them to form the open hydrogen bonded network that is characteristic of gas hydrates. If the adsorption coverage is sufficiently great and the potential energy barrier to rearrangement of this surface structure is sufficiently high, then the surface will be kinetically stabilized and the growth of a gas hydrate crystal will result. Work is in progress to quantity this model. Acknowledgment. The author thanks Dr. J. T. K. Tan for his helpful comments during the development of this work, and British Gas plc for permission to publish. (53) Vysnianskas, A.; Bishnoi, P. R. Chem. Eng. Sei. 1985, 40, 299. (54) Englezos, P.; Kalogerakis,N.; Dholabhai, P. D.; Bishnoi, P. R. Chem. E m . Sei. 1987. 42. 2647. 7 5 5 ) Englezk P.; Kalogerakis,N.; Dholabhai, P. D.; Bishnoi, P. R. Chem. Eng. Sei. 1987, 42, 2659.
