Heterogeneous Crystal Growth of Methane Hydrate on Its sII [001

Feb 5, 2008 - This paper presents a systematic molecular simulation study of the heterogeneous crystal growth of methane hydrate sII from supersaturat...
1 downloads 0 Views 716KB Size
Download PDF
J. Phys. Chem. B 2008, 112, 2399-2404

2399

Heterogeneous Crystal Growth of Methane Hydrate on Its sII [001] Crystallographic Face Jenel Vatamanu and Peter G. Kusalik* Department of Chemistry, UniVersity of Calgary, 2500 UniVersity DriVe NW, Calgary, Alberta, Canada T2N 1N4 ReceiVed: September 20, 2007; In Final Form: December 7, 2007

This paper presents a systematic molecular simulation study of the heterogeneous crystal growth of methane hydrate sII from supersaturated aqueous methane solutions. The growth of sII hydrate on the [001] crystallographic face is achieved through utilization of a recently proposed methodology, and rates of crystal growth of 1 Å/ns were sustained for the molecular models and specific conditions employed in this work. Characteristics of the crystals grown as well as properties and structure of the interface are examined. Water cages with a 51263 arrangement, which are improper to both sI and sII structures, are identified during the heterogeneous growth of sII methane hydrate. We show that the growth of a [001] face of sII hydrate can produce an sI crystalline structure, confirming that cross-nucleation of methane hydrate structures is possible. Defects consisting of two methane molecules trapped in large 51264 cages and water molecules trapped in small and large cages are observed, where in one instance we have found a large 51264 cage containing three water molecules.

1. Introduction Methane hydrates (MH) are crystalline inclusion compounds of methane and water, in which water molecules form cages around methane.1 Great interest in methane hydrates is being stimulated by several important considerations, including the fact that they contain the most abundant source of carbon-based energy on Earth.2,3 Additionally, these compounds are of relevance in global warming,4 geology,5 astronomy6 and marine systems.7 The most common gas hydrate structures are the so-called sI, formed by dodecahedra (512) and tetradecahedra (51262) cages, and sII, containing 512 and larger hexadecahedra (51264) water cages.1 Methane hydrates are stable at relatively high pressure.8 The most stable form of methane hydrate at low to moderate pressures is sI, where the relatively small sized methane molecules stabilize the 512 and 51262 water cages. The larger 51264 cages of sII gas hydrates are normally stabilized by larger inclusion molecules (e.g., tetrahydrofurane).9 However, it has been experimentally shown that at relatively high-pressure structure sII can form for smaller guest molecules like methane.10 Also, in situ transformations between sI and sII have been reported for methane hydrates.11 In addition, one of our recent studies12 presented findings of an interconversion from sI to sII during the heterogeneous crystal growth of the [001] face of MH-sI. Previous simulation studies exploring molecular level properties of methane hydrates have addressed issues such as the interfacial properties or interfacial relaxation of methane hydrate crystals,13 melting,14 properties of the crystalline hydrate phase,15 hydrate stability and their unusual self-preservation,16 the influence of inhibitors on the methane hydrates,17 and the nucleation of methane hydrate cages.18 As well, larger scale and more phenomenological models based on phase field theory19 have been used to predict the (macroscopic) kinetics of CO2 and CH4 gas-hydrate crystal growth,20 interconversion * Corresponding author. E-mail: [email protected].

of CH4 to CO2 hydrate,21 and the morphological structure of CO2 hydrate.22 In a previous article,12 we had presented a simulation framework that was effective in achieving steady-state growth of a methane hydrate. In that work, Newtonian dynamics (i.e., with no thermostat forces) were used to sample the molecular behavior at a growing interface by controlling the temperature in the system with two local thermostats acting far away from the crystal/liquid interface (see ref 23 for further details). Heterogeneous crystal growth was driven at steady-state conditions via an imposed temperature gradient within the system. Although the interface was evolved with real (Newtonian) dynamics, this simulation setup has the disadvantage of relying upon relatively large temperature gradients that, in principle, can affect the quality of the measured properties across the interface and can enhance the demixing of methane from water.12 In this paper, we will employ a methodology24 that performs a canonical sampling at and near the interface, while maintaining a constant temperature along the direction of asymmetry in the system. In this paper, we report the results of a set of molecular dynamic simulations of the [001] crystallographic face of MHsII where steady-state crystal growth from supersaturated solutions of methane and water is achieved. Both constant pressure and constant volume conditions are examined, and we are successful in observing steady-state growth at two different temperatures and various pressures. For the methane supersaturation levels employed, the kinetics of growth are essentially independent of the applied pressure over the range 0 and 2000 atm. The MH-sII structure is found to be stable at these pressures over timescales of several tens of nanoseconds. Several interesting structural features, including [110] microfacets, are found. Structural defects consisting of large cages containing two methane molecules are identified; an instance of three water molecules trapped in one large 51264 cage and another of a single water molecule in a small cage are also noted. Additionally,

10.1021/jp077583k CCC: $40.75 © 2008 American Chemical Society Published on Web 02/05/2008

2400 J. Phys. Chem. B, Vol. 112, No. 8, 2008

Vatamanu and Kusalik TABLE 1: Summary of the MD Trajectories Performeda

Figure 1. MD simulation setup employed in this study. The red area indicates the region where the temperature pulse is applied, and the horizontal lines indicate where the active transport is operative.

one trajectory exhibits growth of sI methane hydrate (which is expected to be more stable) on a [001] MH-sII face. The remainder of this paper is organized as follows. In section 2, a brief description of the simulation methodology is given. Our results and discussions will be presented in section 3, with concluding remarks given in section 4. 2. Methodology In contrast with the setup employed in our earlier work,12 where the heterogeneous crystal growth was controlled using temperature gradients in the system, the molecular dynamics (MD) simulations reported in this paper were performed using the methodology described in ref 24. A schematic of the present system setup is shown in Figure 1. Periodic boundary conditions are applied in all three Cartesian directions, where the system consists of two phases, liquid and crystal, separated by two interfaces. Along the direction of asymmetry (assigned to be our z direction) the temperature profile was constant across the entire system except for a small region near one interface (the melting interface) where a temperature pulse, with a maximum of 355 K, was applied. While melting occurs at the interface with the temperature pulse present, the constant undercooling of the remainder of the system induces crystallization at the other interface. If the rate of crystallization at the one interface equals the rate of melting at the other, then steady-state conditions (crystallization) can be achieved. In the present study, the extent of methane supersaturation and temperature undercooling is used to provide the thermodynamic driving force for crystal growth. To help sustain the former, a small force (less than 2% of the average force experienced by a methane molecule) called “active transport”, was applied near the melting interface to shuttle methane molecules beyond a small region near this interface (and toward the crystallizing interface) in a fashion similar to ref 12. The values of the width and cutoff distance of the Gaussian weighting function (see eq 7 of ref 24) were 16 Å and 15 Å, respectively. We emphasize that both the active transport and the temperature pulse are applied far from the interface of interest at which the crystallization is occurring; therefore, the impact of any nonstandard terms due to active transport is localized to regions of the system that are not the focus of detailed investigation or analysis. We remark that we have carefully verified that the instantaneous linear momentum is rigorously conserved using procedures similar to those in refs 12 and 24. The initial structure of the methane hydrate sII was taken from the X-ray crystallographic data for the oxygen and carbon atoms.25 The hydrogen atoms of the water molecules were assigned according to Bernal-Fowler rules,26 such that any water molecule always accepts and donates two hydrogen bonds

run P (ensemble) (atm)

T (K)

1 (NpT) 2 (NpT) 3 (NpT) 4 (NpT) 5 (NpT) 6 (NpT) 7 (NpT) 8 (NpT) 9 (NVT) 10 (NVT)

255 255 255 255 245 245 245 245 255 245

0 100 500 1000 0 100 500 1000 2000 1800

average crystal width CH4 [CH4]liquid grown of liquid occupancy (molecules/ (Å) region (Å) (%) nm3) 32.1 31.0 23.9 23.7 28.7 29.0 29.8 30.4 28.6 23.0

44.6 44.5 53.1 51.4 44.4 46.7 52.0 54.2 55.5 58.0

84 83 84 84 79 77 87 77 93 -

0.9 1.4 1.4 1.4 1.4 1.4 1.4 1.5 1.6 -

a It should be noted that since the average rate of growth is 1 Å/ns for all runs, the numerical values in column 4 also represent the run lengths in nanoseconds. We also point out that in the evaluation of the average liquid densities for methane (last column) only bulk solution regions well away from the influence of active transport were considered.

but otherwise has random orientation. From 10 000 different random hydrogen-bonded networks initially constructed, the one with the smallest total dipole moment was selected; in this way, a system with a polarization of less than 1% (i.e., its dipole moment was less than 1% of its maximum possible value) was generated. This initial crystal was relaxed for about 1 ns at the desired working temperature and pressure, then about 30% of the methane in the system was removed from a specific slab, and this region melted. Thus, two-phase systems, formed from a methane-water solution and a methane hydrate crystal, were obtained. The long-range interactions in the present simulations were handled with smooth particle mesh Ewald27 using interpolating splines of sixth order, a Fourier space grid of about 1 point/Å, and a real space cutoff of 10 Å with an alpha-parameter of 0.3 Å-1. The short-range interactions were evaluated within a cutoff of 10 Å. The equations of motion were integrated by a fourthorder predictor-corrector method28 using a time step of 1 fs. Energy conservation was tested, and the drift was found to be less than 10-5 kJ/mol per time step. A Berdensen scheme29 (with a relaxation time of 0.4 ps) was employed to control the pressure and separate Nose´-Hoover chain thermostats30 (of order 4) were attached to the translations and rotations of each molecule (with relaxation times of 0.3 ps) to generate the desired average temperature profile in the system. Further details of the MD methodology employed can be found in ref 24. The simulated systems consisted of 4352 molecules of water and 490 molecules of methane confined into a tetragonal box with initial dimensions of 34.2 Å × 34.2 Å × 136.5 Å. The water molecules were modeled using the TIP4P31 potential and the methane molecules were represented by single LJ sites,32 where this potential has been previously used in a variety of successful simulation studies of methane hydrates13b,14b,16b,18a,b. The cross terms for the LJ parameters were determined using the Lorentz-Berthelot mixing rules.28 Experience and testing with various models and systems has indicated that the results presented here will be essentially invariant to reasonable choices of potentials. The liquid-like and solid-like characters of molecules were determined on the basis of their translational mean-square displacements (with respect to averaged positions) during 50 ps trajectory segments.33 All of the configurations presented below will utilize averaged particle coordinates.

Crystal Growth of Methane Hydrate

J. Phys. Chem. B, Vol. 112, No. 8, 2008 2401

Figure 2. Initial and final configurations for MD trajectories 5 and 6. The oxygen atoms belonging to water molecules behaving as liquidlike and solid-like particles are represented by red and light blue spheres, respectively. The lines joining oxygen atoms represent nearest-neighbor molecules within a 3 Å cutoff. The methane molecules characterized as liquid-like are colored blue, while solid-like methane molecules appear as indigo spheres. (a) The initial configuration for the two trajectories; (b,c) The final configurations for runs 5 and 6, respectively. The two vertical yellow lines delimit the regions of new crystal grown.

3. Results and Discussions The simulations reported here were performed at temperatures of 245 and 255 K under both constant volume and constant pressure conditions, where pressures ranged from 0 to 2000 atm. The specifics of these simulation runs are summarized in Table 1. Between 23-32 Å of new crystal were formed within each individual 23-32 ns trajectory; Figure 2 illustrates representative amounts of crystal grown from aqueous methane solutions for these systems. We might expect the kinetics of crystallization to depend upon the physical conditions of the system, that is, its pressure, temperature and supersaturation level of methane in solution. In this study, (average) rates of crystal growth were consistently maintained at 1 Å/ns for all trajectories. High supersaturations of the aqueous methane solution, roughly 3-4 times the expected amounts of dissolved methane under 1000 atm pressure (see Table 1), provided the large driving force required to crystallize the methane hydrate within the time scale accessible to molecular simulations (tens of nanoseconds). The ratio between the average concentration of methane in the crystalline and liquid phases was typically between 3 and 3.5. We did not observe a systematic dependence of the kinetics of growth with pressure or temperature because of the relatively large uncertainties inherent in our results; details of rate dependences would require simulations 10-100 times longer. Consistent with experimental results,8 the crystals formed were not fully occupied with methane; that is, the newly formed hydrate contained empty cages. Methane cage occupancies between 73 and 95% were obtained (see Table 1). Profiles of densities and potential energies across the interface (i.e., along the direction of heterogeneity in the system) were measured

Figure 3. Profiles of (a) densities, (b) densities derivatives, (c) potential energy, and (d) potential energy derivatives across the growing interface. The black and blue lines are data from trajectory 7 for water (left axis) and methane (right axis), respectively. The units of energy are kJ/mol and molecules/Å3 for densities, while distances are in units of the box length. The green vertical line indicates roughly the position of the growing interface.

with respect to a moving frame, whose displacement matches the average rate of growth (i.e., the movement of the temperature pulse). The crude profile data from a MD simulation was smoothed using a 3-point Savitzky-Golay method34 such that numerical derivatives of interfacial profiles could be extracted. As evident from Figure 3, the density and potential energy derivative profiles of methane and water appear somewhat asymmetric and unlike the apparent Gaussian forms observed with one-component systems.23,24 In the case of the energy profile for methane, the decrease in potential energy is slower through the liquid side of the interface and apparently sharper on the crystal side of the interface. We also find that the methane function is shifted further toward the liquid side of the interface,

2402 J. Phys. Chem. B, Vol. 112, No. 8, 2008

Figure 4. Molecular structure at the solid/liquid interface of a hydrate crystal. The atoms are labeled as in Figure 2, where the lines joining oxygen atoms represent molecules with separations of less than 3 Å. The yellow lines emphasize the microfacets exposed toward the liquid. This averaged configuration is taken from trajectory 8.

by 2-4 Å, than is the result for water. This observation suggests “ordering” of methane molecules proceeds that of water, and this behavior can be seen as being associated with the strong tendency of methane to move from the liquid phase into the empty “pockets” formed at the hydrate surface (see ref 12). As in prior work,12,23,24 we use the same definitions here for interfacial width and position: the distance across the interface that encompasses 95% of the peak in the function defines the interfacial width, while the position of the maximum in the derivative defines the position of the interface. The density profiles of methane (and their derivatives) are similar in shape to those of water, and both predict similar widths and positions of the interface. We did not identify a systematic dependence of interfacial widths on pressure or temperature. The average interfacial width from density derivative profiles of both methane and water is in the range of 11-12 Å. The interfacial widths predicted by the derivatives of the potential energy profiles are larger for methane than for water, with values of 14-17 Å and 12-13 Å, respectively. The profile data obtained here for sII methane hydrates, while superior to that reported previously12,35 for sI hydrate (because of the reduction of finite-size effects resulting from larger systems and better sampling), indicates that the two structures have very similar interfacial characteristics. Several comments with respect to measures of interfacial width are in order. The solid/liquid interface of ice Ih is known to be molecularly rough, with this roughness extending over 4-6 molecular diameters.36 This value is consistent with previous simulations of a variety of systems, including methane hydrates,12 atomic systems,23,24 and ice.37 It has also been noted in these simulation studies that the specific value of the interface width depend on the measure (e.g., energy, density, or structure) employed in its determination. In addition, in very recent work with atomic systems,38 we have shown that the solid/liquid interface has a tendency to exhibit collective structural fluctuations that favor exposure of particular microfacets. Hence, we would expect any measurement of the interfacial width in such systems also to depend on the time and length scales over which averages are performed (i.e., the structural fluctuations that the interface is allowed to sample). We remark that in this picture of crystal growth it is not appropriate to imagine the deposition of individual water molecules from the liquid phase into the ordered hydrogen-bonded network of the hydrate crystal. In this

Vatamanu and Kusalik

Figure 5. Image looking down onto a [001] face of MH-sII during heterogeneous crystal growth, omitting all liquid-like molecules from view. The formation of both 51264 and 51263 cages at the interface is apparent. The yellow shaded area denotes a 51264 cage formed at the interface, where the large yellow sphere is the trapped methane molecule. The methane molecules trapped in surface 51263 cages (shaded light-blue) are represented by smaller green spheres. The remaining molecules are represented as in Figure 2, where the lines joining oxygen atoms represent water molecules with separations of less than 3 Å. This averaged configuration is taken from trajectory 9.

work, we have chosen simply to record the measured property along the direction of growth (the z-direction) and have made no attempt to try to factor out the complex interfacial topology that systems might exhibit at any instant in time. The interfacial widths measured here thus might be expected to correspond to results from experiments that probe nanometer length scales on nanosecond time scales. The structure of the crystal/liquid interface consists of a complex and continually changing pattern generated by the three-dimensional (3-D) arrangement of the partially formed water cages within the interface. Instead of maintaining a planar character, where each layer through the system exhibits a consistent degree of structure, interfacial layers display a varied blend of ordered and disordered local regions. In previous work with atomic systems, we have shown that the liquid/crystal interface favors (on average) a topology that exposes specific microfacets toward the liquid.38 In this work, we identified a similar character within the crystal/solution interface of MHsII, where the [001] interface tends to exhibit [110] microfacets (see Figure 4). We emphasize that this is an average tendency; thus, at any particular instant in time, one will find a varying detailed topology as a result of the inherent fluctuating nature of the interface. The much larger cross-sectional area employed in this study relative to that of our previous work with sI hydrate12 has greatly enhanced the extent of microfaceting observed, consistent with the behavior found in atomic systems.35 On several occasions across several of the trajectories, cages of the wrong type were observed to form at the interface. It has been previously reported that during growth of the [001] face of MH-sI, intermediate 51263 cages can form, where these cages do not belong to either MH-sI or MH-sII. In this work, we observed that growth of the [001] face of MH-sII can also lead to 51263 cages (see Figure 5). If these 51263 cages survive the ordering-disordering fluctuations at the interface, then the growth of a different crystal structure will be propagated (either

Crystal Growth of Methane Hydrate

Figure 6. Image of part of the crystal grown in trajectory 10 capturing formation of sI during the growth of hydrate sII. The solid within the shaded area corresponds to sI hydrate. The orientation of the system has been chosen to emphasize the sI pattern. The molecules are represented as in Figure 2, where again only solid-like particles are shown.

Figure 7. Examples of observed crystal defects: (a) two methane molecules (blue spheres) trapped in a large cage from run 7, and (b) three molecules of water (red spheres) trapped in a large cage from run 10. These images are averaged configurations (over 50 ps) and capture the realistic positions only of those atoms making the cage but not for the molecules trapped within the cages. The light blue spheres are the oxygen atoms of the water molecules forming the cages, where the lines joining the spheres represent water molecules with separations of less than 3 Å.

sI or sK12); alternatively, these cages will need to melt and be replaced by the large 51264 cages specific to a sII hydrate. The fluctuations occurring at an interface can result in crystalline structures that are different from that of the crystalline seed provided.12,35 When such structures can be sterically accommodated on the initial crystal template, a new kind of crystal (or polymorph) can be grown. It has been previously shown12 that a [001] face of sI hydrate can be successfully transformed into a [110] face of sII, where the two structures are made compatible via a transition layer of intermediate 51263 cages. Therefore, one would expect that the reverse should be a possibility, that is, the growth of MH-sI on a [001] MH-sII face. Indeed, we were able to identify such a behavior in one of our trajectories; this situation is shown in Figure 6 where a relatively large chunk of sI hydrate has formed during the growth of sII. Since the cross-sectional area used in this work is considerably larger than that in our earlier work,35 we suggest that such polymorph transitions (i.e., sI-sK-sII) during the heterogeneous crystal growth of methane hydrates are not dependent on finite-size effects. Finally, as in our previous growth studies of sI hydrate, we also observe significant numbers of defects in the crystals we have grown, apparently independent of the imposed conditions. As noted above, cage occupancy levels were typically about

J. Phys. Chem. B, Vol. 112, No. 8, 2008 2403 80%. Yet, some of our simulations resulted in two methane molecules (rather than one) being trapped in the large 51264 cages of MH-sII (see Figure 7a. Trajectories 6, 7, 9, and 10 produced 2, 1, 4, and 4 double methane occupied cages, respectively. While the majority of the double occupancies were found in runs performed at higher pressures, this data was not conclusive. Intriguingly, in one instance (trajectory 10), we identified a large 51264 cage containing three trapped water molecules (see Figure 7b). Molecular exchanges between the water molecules within the cage and those making up the cage were also observed on three separate occasions over a period of 9 ns. In each case, the swapping events involved different individual molecules and occurred over timescales of less than 100 ps. In trajectory 6, we identified one small 512 cage occupied with a single water molecule. Clearly, such defect structures, even at modest concentrations, could have a significant impact on the properties of hydrate crystals. 4. Conclusions In this paper, results of several molecular simulations of the heterogeneous crystal growth of methane hydrate sII are presented. Methane hydrate crystals were grown under steadystate conditions both at constant pressure and at constant volume at rates of 1 Å/ns. The MH-sII/aqueous methane solution interface was found to be between 11 and 12 Å wide from the point of view of particle density, and 12-17 Å wide from the potential energy viewpoint. The measured profile functions for potential energies and densities across the interface were not symmetric, and they did not exhibit a temperature or pressure dependence. Consistent with previous work, the topology of the sII hydrate interface was found to be complex with the [001] face exhibiting [110] microfacets. Intermediate 51263 cages and a polycrystalline transition from MH-sII to MH-sI via such cages were identified. Crystals grown were found to contain a number of defects. Typically about 20% of the hydrate cages were unoccupied, apparently independent of the applied temperature and pressure, while defects consisting of cages containing two methane molecules trapped in a large 51264 cage were also seen, primarily in simulations at higher pressures. Additionally, cages containing one or three water molecules (instead of methane) were observed. It is our hope that future careful experimental studies, perhaps utilizing X-ray or neutron scattering, may be able to verify the presence of such defects in methane hydrate crystals. Although the populations of these defects can be expected to be low, they may have reasonable impacts upon properties such as thermal conductivities. This work further emphasizes that the crystal growth of systems such as methane hydrates is rich with molecular level detail and demonstrates that molecular simulation is an excellent means by which to probe this behavior. Certainly, an understanding of observed macroscopic properties will be built on such molecular level insights. Our future work will explore CO2 and tetrahydrofurane hydrates as well as investigate further the growth of methane hydrates. Acknowledgment. We are grateful for the financial support of the Natural Sciences and Engineering Research Council of Canada and the University of Calgary. We also acknowledge computational resources provided by WestGrid. References and Notes (1) Atwood, J. L.; Davies, J. E. D.; MacNicol, D. D. Structural Aspects of Inclusion Compounds Formed by Inorganic and Organic Host Lattices; Inclusion Compounds, Vol. 1; Oxford University Press: Oxford, UK, 1991. (2) Sloan, E. D., Jr. Nature 2003, 426, 353.

2404 J. Phys. Chem. B, Vol. 112, No. 8, 2008 (3) (a) Claypool, G. E.; Kaplan, I. R. Natural Gases in Marine Sediments; Plenum: New York, 1974; Vol. 3, pp 99-139. (b) Barnes, R. O.; Goldberg, E. D. Geology 1976, 4, 297. (c) Boetius, A.; Suess, E. Chem. Geol. 2004, 205, 291. (d) Kevnvolden, K. A. Chem. Geol. 1988, 71, 41. (e) Chapman, P. K.; Haynes, W. Acta Astron. 2005, 57, 372. (f) Chatti, I.; Delahaye, A.; Fournaison, L.; Petitet, J. P. Energy ConVers. Manage. 2005, 46, 1333. (g) Pooladi-Darvish, M. J. Pet. Technol. 2004, 56, 65. (h) Laherrere, J. Energy Explor. Exploit. 2000, 18, 349. (4) (a) Wignall, P. B.; Newton, R. J.; Little, C. T. S. Am. J. of Science 2005, 305, 1014. (b) Chazelas, B.; Leger, A.; Ollivier, M. Sci. Total EnViron. 2006, 354, 292. (c) McElwain, J. C.; Wade-Murphy, J.; Hesselbo, S. P. Nature 2005, 435, 479. (d) Buffett, B.; Archer, D. Earth Planet. Sci. Lett. 2004, 227, 185. (e) Beauchamp, B. Comptes Rendus Geoscience 2004, 336, 751. (f) Svensen, H.; Planke, S.; Malthe-Sorenssen, A. et al. Nature 2004, 429, 542. (g) Renssen, H.; Beets, C. J.; Fichefet, T. et al. Paleoceanography 2004, 19, Art. No. PA2010. (h) Dickens, G. R. Earth Planet. Sci. Lett. 2003, 21, 169. (i) Benton, M. J.; Twitchett, R. J. Trends in Ecology & EVolution 2003, 18, 358. (j) Judd, A. G.; Hovland, M.; Dimitrov, L. I. et al. Geofluids 2002, 2, 109. (k) Wignall, P. B. Earth-Sci. ReV. 2001, 53, 1. (l) Brewer, P. G. Ann. N. Y. Acad. Sci. 2000, 912, 195. (m) Paull, C. K.; Ussler, W.; Dillon, W. P. Geophys. Res. Lett. 1991, 18, 432. (5) (a) Kleinberg, R. L.; Flaum, C.; Griffin, D. D. et al. J. Geophys. Res. 2003, 108 (B10), Art. No. 2508. (b) Fluteau, F. Comptes Rendus Geoscience 2003, 335, 157. (c) Rohl, U.; Bralower, T. J.; Norris, R. D., et al. Geology 2000, 28, 927. (d) Gornitz, V.; Fung, I. Global Biogeochem. Cycles 1994, 8, 335. (e) Kvenvolden, K. A. ReV. Geophys. 1993, 31, 173. (6) (a) Lunine, J. I.; Stevenson, D. J. Icarus 1987, 70, 61. (b) Lunine, J. I.; Stevenson, D. J. Astrophys. J. Suppl. Ser. 1985, 58, 493. (c) Grasset, O.; Sotin, C.; Deschamps, F. Planet. Space Sci. 2000, 48, 617. (7) (a) Kelleher, B. P.; Simpson, A. J.; Rogers, R. E., et al. Marine Chem. 2007, 103, 237. (b) Farabegoli, E.; Perri, M. C.; Posenato, R. Global and Planetary Change 2007, 55, 109. (c) Chuang, P. C.; Yang, T. F.; Lin, S. et al. Terrestrial Atmospheric and Oceanic Sciences 2006, 17, 903. (d) Matsushima, J. J. Geophys. Res. 2006, 111 (B10), Art. No. B10101. (e) Hill, T. M.; Kennett, J. P.; Valentine, D. L. et al. PNAS 2006, 103, 13570. (f) Leifer, I.; Luyendyk, B. P.; Boles, J. et al. Global Biochemical Cycles 2006, 20, Art. No. GB3008. (g) Inagaki, F.; Nunoura, T.; Nakagawa, S. et al. PNAS 2006, 103, 2815. (h) Priest, J. A.; Best, A. I.; Clayton, C. R. I. Geophys. J. Int. 2006, 164, 149. (8) Sloan, E. D., Jr. Clarthrate Hydrates of Natural Gases, 2nd ed.; Marcel Dekker: New York, 1998. (9) Mak, T. C. W.; McMullan, R. J. Chem. Phys. 1964, 42, 2732. (10) Loveday, J. S.; Nelmes, R. J.; Guthrie, M.; Belmonte, S. A.; Allan, D. R.; Klug, D. D.; Tse, J. S.; Handa, Y. P. Nature 2001, 410, 661. (11) (a) Shimizu, H.; Kumazaki, T.; Kume, T.; Sasaki, S. J. Phys. Chem. B 2002, 106, 30. (b) Chou, I. M.; Anurag Sharma, A.; Burruss, R. C.; Shu, J.; Mao, H. K.; Russell, J.; Hemley, R. J.; Goncharov, A. F.; Stern L. A.; Kirby S. H. PNAS. 2000, 97, 13484. (12) Vatamanu, J.; Kusalik, P. G. J. Phys. Chem. B 2006, 110, 15896. (13) (a) Rodger, P. M. Ann. N. Y. Acad. Sci. 2000, 912, 474. (b) Nada H. J. Phys. Chem. B 2006, 110, 16526. (14) (a) Ripmeester, J. A.; Ratcliffe, C. I.; Klug, D. D.; Tse, J. S. Ann. N. Y. Acad. Sci. 1994, 715, 161. (b) Forrisdahl, O. K.; Kvamme, B.; Haymet, A. D. J. Mol. Phys. 1996, 89, 819. (c) Rodger, P. M.; Forester, T. R.; Smith, W. Fluid Phase Equilib. 1996, 116, 326. (d) English, N. J.; Johnson, J. K.; Tylor, C. E. J. Chem. Phys. 2005, 123, 244503. (15) (a) Alavi, S.; Ripmeester, J. A.; Klug, D. D. J. Chem. Phys. 2006, 124, 014704. (b) Tse, J. S.; Klein, M.; McDonald, I. R. J. Chem. Phys. 1984, 81, 6146. (c) Sizov V. V.; Piotroskaya, E. M. J. Phys. Chem. B 2007,

Vatamanu and Kusalik 111, 2886. (d) Miyoshi, T.; Ohmura, R.; Yasuoka, K. J. Phys. Chem. C 2007, 111, 3799. (e) English, N. J.; MacElroy, J. M. D. J. Comput. Chem. 2003, 24, 1569. (f) Jiang, H.; Jordan, K. D.; Taylor, C. E. J. Phys. Chem. B 2007, 111, 6486. (16) (a) Tse, J. S.; Klug, D. D. J. Supramol. Chem. 2002, 2, 467. (b) Alavi, S.; Ripmeester, J. A.; Klug, D. D. J. Chem. Phys. 2007, 126, 124708. (17) (a) Carver, T. J.; Drew, M. G. B.; Rodger, P. M. J. Chem. Soc., Faraday Trans. 1995, 91, 3449. (b) Kvamme, B.; Huseby, G.; Forrisdahl, O. K. Mol. Phys. 1997, 90, 979. (c) Wallqvist, A. J. Chem. Phys. 1992, 96, 5377. (d) Storr, M. T.; Rodger, P. M. Ann. N. Y. Acad. Sci. 2000, 912, 669. (e) Anderson, B. J.; Tester, J. W.; Borghi, G. P.; Trout, B. L. J. Am. Chem. Soc. 2005, 127, 17852. (18) (a) Moon, C.; Taylor, P. C.; Rodger, P. M. J. Am. Chem. Soc. 2003, 125, 4706. (b) Guo, G. J.; Zhang, Y. G.; Zhao, Y. J. J. Chem. Phys. 2004, 121, 1542. (c) Guo, G. J.; Zhang, Y. G.; Liu, H. J. Phys. Chem. C 2007, 111, 2595. (19) Granasy, L.; Pusztai, T.; Warren, J. A. J. Phys.: Condens. Matter 2004, 16, R1205, and references therein. (20) Svandal, A.; Kvamme, B.; Granasy, L.; Pusztai, T.; Buanes, T.; Hove, J. J. Cryst. Growth 2006, 287, 486. (21) (a) Moon, C.; Hawtin, R. W.; Rodger, P. M. Faraday Discuss. 2007, 136, 367. (b) Lee, H.; Seo, Y.; Se, Y.-T.; Moudrakovski, I. L.; Ripmeester, J. A. Angew. Chem., Int. Ed. 2003, 42, 5048. (c) Park, Y.; Kim, D.-Y.; Lee J. W.; Huh, D.-G.; Park, K.-P.; Lee, J.; Lee, H. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 12690. (d) Tegze, G.; Granasy, L.; Kvamme, B. Phys. Chem. Chem. Phys. 2007, 9, 3104. (22) Tegze, G.; Pusztai, T.; Toth, G.; Granasy, L.; Svendal, A.; Buanes, T.; Kuznetsova, T.; Kvamme, B. J. Chem. Phys. 2006, 124, 234710. (23) Gulam Razul, M. S.; Tam, E. V.; Lam, M. E.; Linden, P.; Kusalik, P. G. Mol. Phys. 2005, 103, 1929. (24) Vatamanu, J.; Kusalik, P. G. J. Chem. Phys. 2007, 126, 124703. (25) McMullan, R. K.; Jeffrey, G. A. J. Chem. Phys. 1965, 42, 2725. (26) (a) Bernal, J. D.; Fowler, R. H. J. Chem. Phys. 1993, 1, 515. (b) Hayward, J. A.; Reimers, J. R. J. Chem. Phys. 1997, 106, 1518. (27) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. J. Chem. Phys. 1995, 103, 8577. (28) Allen, M. P.; Tildesley, D. J. Computer Simulations of Liquids, Clarendon Press: Oxford, 1993; pp 340-342. (29) Berendsen, H. C. J.; Van Gunsteren, W. F. Molecular Dynamics simulations: techniques and approaches. In Molecular Liquids, dynamics and interactions; NATO ASI Series C135; Reidel: New York, 1984. (30) (a) Nose, S. J. Chem. Phys. 1984, 81, 511. (b) Hoover, W. G. Phys. ReV. A 1985, 31, 1695. (c) Martyna, G. J.; Klein, M. L.; Tuckerman, M. J. Chem. Phys. 1992, 97, 2635. (31) Jorgensen, W. L.; Chahdrasekhar, J.; Madura, J. D.; Impery, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (32) Jorgensen, W. L.; Madura, J. D.; Swenson, C. J. J. Am. Chem. Soc. 1984, 106, 6638. (33) Gulam Razul, M. S.; Hendry, J. G.; Kusalik, P. G. J. Chem. Phys. 2005, 123, 204722. (34) Press, W. H.; Flannery, B. P.; Teukolsky, S. A.; Vetterling, W. T. Numerical Recipes in Fortran: The Art of Scientific Computing; Cambridge University Press: New York, 1992. (35) Vatamanu, J.; Kusalik, P. G. J. Am. Chem. Soc. 2006, 128, 15588. (36) Beaglehole, D.; Wilson, P. J. Phys. Chem. B 1993, 973, 11053. (37) Gulam Razul, M. S. Ph.D. Thesis, Dalhouise University, Halifax, 2005. (38) Vatamanu, J.; Kusalik, P. G. Phys. ReV. B 2007, 76, 035431.