Computer Simulation of Methane Hydrate Cage Occupancy

2886

J. Phys. Chem. B 2007, 111, 2886-2890

Computer Simulation of Methane Hydrate Cage Occupancy Vladimir V. Sizov* and Elena M. Piotrovskaya Department of Chemistry, St. Petersburg State UniVersity, UniVersitetskii prospekt, 26, Peterhof, 198504 St. Petersburg, Russia ReceiVed: September 9, 2006; In Final Form: NoVember 11, 2006

Grand canonical Monte Carlo simulation results are presented for bulk sI methane hydrate. The description of hydrogen-bonded clathrate network allows the water molecules to move and rotate. A more idealized rigid structural model based on the van der Waals-Platteeuw (vdWP) theory is used for comparison. Occupancy isotherms and pressure versus temperature occupancy diagrams are computed for temperatures below 260 K and pressures up to 400 bar. It is found that the results obtained with the vdWP-like model are in qualitative agreement with experiment, though this model fails to account for structural transformations of water network in the vicinity of the melting point.

Introduction Gas hydrates are crystalline inclusion compounds consisting of gas molecules encaged in cavities of a hydrogen-bonded network of water molecules. These compounds can exist at elevated pressure for temperatures close to or below the melting point of ice. Depending on the properties of the guest gas molecules and the details of hydrate formation procedure, different gas hydrate structures can be obtained, the most important structures being sI, sII, and sH.1 Typically, smaller gas molecules (such as methane, ethane, and carbon dioxide) tend to form sI hydrates, while larger molecules preferentially form sII (propane, iso-butane) and sH (cyclohexane, cycloheptane) hydrates. These structures differ in the size of cavities in the clathrate network of water molecules, as well as in the number of cavities of different types in the unit cell. The smallest cavity found in gas hydrates is the pentagonal dodecahedral cage (512) comprising twelve pentagons. Larger cavities include tetracaidecahedral (51262) and hexacaidecahedral (51264) cage, which can be found in sI and sII hydrates, respectively. The industrial importance of gas hydrates (see refs 1 and 2 and references cited therein) required quantitative theoretical tools for the prediction of hydrate phase equilibria. The most successful of such tools is the theory of van der Waals and Platteeuw (vdWP),3 developed as early as 1958. This theory assumes that interactions between the guest molecule and the aqueous network are relatively weak and are limited only to the nearest neighbors of the gas molecule. The interactions between guest molecules in nearby cages are neglected, while the case of multiple cage occupancy is not considered at all. Thus, the behavior of every guest molecule does not depend upon the presence of other guest molecules. This classical model was later extended by Parrish and Prausnitz4 for multicomponent gas hydrates. More recent modifications of vdWP theory (see, for example, refs 5-7) accounted for the flexibility of hydrogenbonded network, multiple cage occupancy, and a number of other phenomena neglected by the original formulation. However, the original vdWP approach still remains the most important and widely used theoretical tool for the investigation of gas hydrates. * Corresponding author. Telephone: +7(812)4284066. Fax: +7(812)4286939. E-mail: [email protected].

Methane hydrate, which is the system of interest to us in this study, has been receiving a lot of attention in theoretical studies due to the abundance of methane in terrestrial hydrates. The results of theoretical analysis were verified by 13C NMR measurements8 and Raman spectroscopy.9 Despite the simplifications discussed above, good agreement with the predictions of vdWP theory was obtained. The validity of the basic approximations of the vdWP theory was studied in ref 10 from the point of view of computer simulations. The vdWP theory and its modifications were found deficient in the treatment of guest-guest interactions and the influence of cage occupancy on the thermodynamic properties of gas hydrates. However, the most controversial assumption of the vdWP model is the neglect of the relaxation of the hydrogen-bonded network forming the clathrate cages. The host structure of the hydrate in vdWP treatment is considered as a static field in which guest molecules evolve. The validity of this assumption can be examined by straightforward application of computer simulation methods. In particular, an obvious way to test the realism of rigid clathrate approach is to compare the simulated behavior of gas hydrates suggested by two structural models. The first model demands for a rigid clathrate structure to ensure vdWP-like behavior. Since the original theory is largely based on idealized models used in statistical mechanics of adsorption, one can state that the gas molecules are thus “adsorbed” in the external field of the clathrate network. This adsorption is localized since it is limited to the clathrate cages rather than to the whole space occupied by the gas hydrate. The second model goes beyond the rigid clathrate treatment by introducing translational and rotational degrees of freedom for water molecules. The goal of this study is to investigate the dependence of clathrate cage occupancy upon temperature and pressure by means of computer simulations. Simulations for the vdWP-like rigid clathrates will be carried out to allow for a comparison with the predictions of thermodynamic theories. Going beyond vdWP approximations with the flexible clathrate model we can obtain molecular-level information on the possible effects of clathrate framework non-rigidity with particular emphasis on the gas hydrate composition and structure. Simulation Details. Computer simulations were carried out for a system containing eight (2 × 2 × 2) sI methane hydrate

10.1021/jp0658905 CCC: $37.00 © 2007 American Chemical Society Published on Web 02/28/2007

Methane Hydrate Cage Occupancy

J. Phys. Chem. B, Vol. 111, No. 11, 2007 2887

TABLE 1: Water and Methane Model Parameters species water

model

site

SPC/E O H methane UA CH4 AA C H

/kB, K σ, nm 78.2 0 148.1 33.2 15.1

q, e

geometry

3.166 -0.8476 O-H ) 1.0 Å 0 0.4238 H-O-H ) 109.47° 3.73 0 3.5 -0.24 C-Η ) 1.09 Å 2.5 0.06 H-Χ-Η ) 109.47°

unit cells with 3D periodic boundary conditions. The number of water molecules (368 water molecules in eight unit cells) in the simulation cell was kept constant while the number of methane molecules varied with the chemical potential settings. In a manner similar to simulations of adsorption in a static external field, the chemical potential of the host molecules forming the clathrate cages was not considered explicitly. The resulting simulation conditions required the application of the grand canonical ensemble Monte Carlo method (GCMC),11 which is widely used for the simulation of adsorption. This method was also used in ref 5 to calculate the average occupancy in propane hydrate at 273.15 K for several gas pressures and to estimate the free energy of cage occupation. It should be noticed here that most simulation studies of gas hydrates (see refs 10 and 12-16 for some recent studies) were carried out using the molecular dynamics method17 in the canonical (NVT) and isothermal-isobaric (NPT) ensembles. While these simulations are highly valuable for observing gas hydrate behavior at molecular level, and especially for obtaining time-dependent properties, the composition of the system cannot be varied in the simulation run in response to changes in pressure or temperature. Instead, multiple simulations are needed to study the effects of gas content on structural and thermodynamic properties of hydrates. Monte Carlo simulations were carried out in several isothermal series for different values of methane chemical potential. In all series the chemical potential was being gradually decreased. The final configuration obtained in a Monte Carlo run was used as the starting point for the next run with a lower value of the chemical potential. This technique resembles the calculation of adsorption isotherms, which in this paper will be referred to as occupancy isotherms. A typical GCMC simulation run consisted of 10 million Monte Carlo steps (up to 50 million for high occupancy). This is considered sufficient for the evaluation of cage occupancy even for the flexible clathrate model, especially when the change of the chemical potential between two runs is not large. Typically, every occupancy isotherm comprised twelve chemical potential values ranging from -1.2 kJ/mol to -14.8 kJ/ mol. Eight occupancy isotherms were calculated both for rigid and flexible clathrate models, the temperatures considered being 200, 210, 220, 230, 235, 240, 250, and 260 K. For convenience the occupancy was plotted versus the pressure of methane in the bulk phase (P) rather than the chemical potential of methane. The gas pressure for a given chemical potential was obtained as a statistical average via virial calculations in Monte Carlo simulations of bulk methane with three-dimensional periodic boundary conditions. Therefore, the values of P should not be directly associated with overall pressure. The SPC/E model18 was chosen for water along with the united-atom (UA) Lennard-Jones model for methane.19 The validity of the UA representation of methane was verified by a series of test calculations with an all-atom (AA) five-site methane model.20 The Lennard-Jones potential parameters, partial charges, and molecular structure information for all molecular models mentioned above are given in Table 1. Longrange interactions were cut off at 20 Å, thus requiring the reconstruction of the periodic images of the simulation cell.

The average number of methane molecules (N) obtained from GCMC was used to calculate the total occupancy θ ) N/64. The cage occupancy is a quantity important both to experimental and theoretical studies of gas hydrate, since it is a measure of the gas content in the hydrate phase. Furthermore, the hydration number (n) can be readily obtained from the total occupancy as n ) 5.75/θ. For the rigid clathrate model, we use the traditional GCMC technique with four types of trial moves: molecule movement, rotation, creation and destruction. The probabilities of trial steps were chosen as 40%, 40%, 10%, and 10%, respectively. For the united atom methane model, particle rotations lead to identical configurations, so only three types of steps were performed: movement (75%), insertion (12.5%), and destruction (12.5%). In principle, even higher probabilities of insertions/ destructions can be justified due to low acceptance ratio of these steps. The application of the GCMC technique to flexible clathrate model is significantly more complicated. First of all, the water molecules, which now possessed rotational and translational degrees of freedom, were not assigned a chemical potential value and were not inserted/deleted by appropriate GCMC steps. Therefore, the GCMC procedure had to be optimized to obtain good sampling for both the number of guest molecules and the hydrate structure corresponding to every value of cage occupancy. The reasonable balance between guest molecules insertion/ deletion and the relaxation of aqueous cages can be achieved by varying the GCMC step probabilities and/or the step acceptance ratio. In several series of preliminary simulations, different combinations of step probabilities were compared. Insertions and destructions of molecules were attempted with step probabilities from 10% to 47.5%. Depending on the cage occupancy (i.e., temperature and chemical potential of guest molecules), the acceptance ratio of these steps ranged from 0.2% to 41.5%, the greater acceptance ratios usually being observed for lower cage occupancy values. Due to the dramatic variations in the acceptance of insertion/destruction steps the GCMC production runs were carried out with fixed step probabilities rather than with fixed relaxation versus insertion/destruction acceptance probabilities ratio, since the latter choice implied the computationally expensive tuning of step probabilities, which, however, could not guarantee the adequate sampling for both water and gas molecules. In preliminary simulations the Monte Carlo steps involving displacement of molecules were carried out with step probabilities from 3% to 40%. For 0.1 Å maximum displacement, the acceptance ratio was typically between 60% and 63%, while higher values of maximum displacement lead to significantly less favorable acceptance statistics: 35% for 0.2 Å and 26% for 0.25 Å. The rotation of molecules was attempted with probability from 2% to 40%; the average acceptance ratio of 68% was observed for 5.7° maximum rotation, 43% for 11.5°, and 37% for 14.3°. Along the occupancy isotherms the acceptance ratios of displacements and rotations decreased with the chemical potential of methane. The final set of step probabilities chosen for GCMC simulations with the flexible clathrate model was identical to the one used for the rigid model. Results and Discussion Typical occupancy isotherms obtained at 200 K for both UA and AA methane models using rigid and flexible clathrate representations are shown in Figure 1. The comparison of

2888 J. Phys. Chem. B, Vol. 111, No. 11, 2007

Figure 1. Typical occupancy isotherms for rigid and flexible clathrate models calculated using UA and AA methane potentials (circles flexible model, UA methane; triangles - flexible model, AA methane; squares - rigid model, UA methane).

Figure 2. log P vs T cage occupancy diagram for the rigid model of clathrate network.

isotherms obtained for UA and AA methane models suggests that the simple Lennard-Jones model can provide reasonable quality of cage occupancy calculations at a significantly lower computational cost than for all-atom representation. The occupancy isotherms presented in Figure 1 reveal the major difference between two clathrate structure models used in this study. For the rigid model, high occupancy (over 90%) is obtained already at relatively low pressure (ca. 10 bar at 200 K). Isotherms for the flexible model indicate significantly lower cage occupancies, which is the immediate result of introducing additional degrees of freedom (movement and rotation of water molecules) to the hydrate model. Though the shape of isotherms in Figure 1 is not unexpected, neither of the two models is able to produce an ideal Langmuir-like occupancy isotherm. More information about the dependence of cage occupancy on pressure and temperature can be obtained from pressure versus temperature occupancy diagram presented in Figure 2 for the rigid clathrate network. The occupancy diagrams were obtained by projecting the equal occupancy contour lines from a three-dimensional occupancy-pressure-temperature plot on the pressure-temperature plane. The contours are plotted with a 0.1 step, the numbers on the contour lines indicating the occupancy for every contour line. For reasons discussed below the contours corresponding to 0.4 and lower occupancies are plotted as dashed rather than solid lines. The vdWP-like description of gas hydrate structure leads to a highly idealized behavior of the system. The occupancy contours seem almost linear in log P versus T coordinates. The spacing between contour lines increases with increasing occupancy, which is an obvious consequence of the shape of

Sizov and Piotrovskaya

Figure 3. log P vs T cage occupancy diagram for the flexible model of clathrate network.

occupancy isotherms presented in Figure 1. As can be readily appreciated, the flexible model of gas hydrate structure shows major deviations from the vdWP-like description (Figure 3). First of all, the occupancy contour lines for the flexible clathrates are no longer straight. In fact, all contours display a sharp bend between 230 and 240 K, while before and after this temperature range the contours are approximately linear. Though the bends at 235 K cannot be directly connected to the melting of the gas hydrate, they still indicate a major change in the gas absorption conditions. The temperature of this transition is fairly close to the SPC/E ice Ih melting point of 225 K observed in computer simulations.21,22 Some recent studies suggest even lower melting points for SPC/E ice Ihs21524-26 or 214 K.27 These melting points were located through rigorous free energy calculations, while in less sophisticated simulations metastable SPC/E ice Ih can be observed up to 300 K.28 This issue severely complicates the comparison of structural changes observed for gas hydrates with the melting of ice. Thus, the nature of the transition under study cannot be determined from superficial analysis of occupancy curves and deserves a more detailed discussion which will be presented later. For any given pressure and temperature, the occupancy of rigid clathrate is higher than for the flexible one. This effect is especially pronounced for higher cage occupancies. The agreement of the flexible model with the vdWP-like model should be significantly better for the temperatures below 200 K, which are not considered in this work. The reasons for the mentioned unidentified transition at 235 K (Figure 3) should be sought in the analysis of structural changes in the clathrate network accompanying the increase of temperature. Since the gas hydrates have lower density than ice Ih, the methane hydrate melting point under the conditions of our simulations should not be strongly dependent upon the spatial constraints of constant-volume simulations. Thus, a temperature of 230 K or even 235 K can be a reasonable estimate of methane hydrate melting point. In contrast with this assumption, our simulations do not show the partitioning of methane and water above this temperature, which prevents us from defining the transition in question as a phase transition. The melting of solid water phases is accompanied by a sharp increase of the mobility of water molecules. In the case of gas hydrates this leads to large-scale deviations from the cage structure corresponding to the global energy minimum. As a result, in the vicinity of the melting transition some cages may be temporarily or permanently “closed” to the guest molecules even though the clathrate structure on the whole is preserved, i.e., no melting occurs. The significant difference between cage

Methane Hydrate Cage Occupancy

J. Phys. Chem. B, Vol. 111, No. 11, 2007 2889

Figure 4. Snapshots of methane hydrate structures with various cage occupancies obtained for UA methane model: (a) θ ) 1.19, (b) θ ) 0.95, (c) θ ) 0.69, (d) θ ) 0.55, (e) θ ) 0.36, and (f) θ ) 0.27. Dashed lines show the hydrogen-bonded water network.

occupancy values for rigid and flexible clathrate models could be a consequence of this effect, which can only occur in simulation when the latter model is used. However, we found no direct evidence of “cage closures” due to excessive relaxation of the clathrate network for the particular simulation technique used in this work. Therefore, the relatively low cage occupancy values obtained for the flexible clathrate model are likely to be caused by the peculiarities of SPC/E water behavior in gas hydrates. Obviously, the SPC/E water model is less successful in reproducing the structure of solid water phases (including gas hydrates) than it is for liquid water. This assumption can be verified by simulations with more accurate water models, which can also provide a better estimate of the melting point of water. Though our simulations do not directly characterize the stability of methane hydrate (free energy, dissociation pressure, etc.), one can attempt to compare the predictions of both clathrate models with experimental data using indirect indicators. The results of Raman and NMR measurements, as well as theoretical calculations, provide information on the hydration number n, which can be compared to simulation results. Since most experimental and theoretical values of n fall between 5.9 and 6.31,8,9,29-33 or at most 6.8,34 in Figures 2 and 3 we should compare the location of regions with θ > 0.8 (n < 7.2) or θ > 0.9 (n < 6.4) on the P-T diagrams (Figures 2 and 3) to the experimentally observed conditions of hydrate stability. This approach suggests the equivalence between the coexistence line on the phase diagram and the contour of 90% total cage occupancy, corresponding to the stable methane hydrate. Though the hydration number may vary along the coexistence line,35 the range of these variations allows one to use a constantoccupancy line for qualitative comparison with experimental phase diagrams.

Choosing n < 6.4 as an indication of methane hydrate being near the coexistence line, one can conclude that the flexible clathrate model predicts that methane hydrate is likely to be stable at relatively high pressures (ca. 160 bar at 200 K). The rigid model produces significantly lower pressure estimates for the existence of stable hydrates (5 bar at 200 K or 55 bar at 260 K), which are in somewhat better agreement with experimental data on hydrate/water+gas equilibrium than the results obtained from the flexible model. For the latter model, the gas pressures observed for total occupancy of 0.9 can be correlated to experimental methane hydrate phase diagrams1,36-38 only if one assumes that the wide (up to 90 K28) metastability region above the SPC/E ice Ih true melting point can be found for sI hydrate as well, i.e., that the simulation temperature of 215 K corresponds to experimental temperature of 300-305 K. Theoretical calculations (see, for example, refs 33 and 39) also do not support the results obtained using the flexible clathrate model. The inability of the SPC/E water model to provide the basis for accurate correlation of simulation temperatures with the experimental ones precludes a quantitative comparison of experimental phase diagrams with simulation data. Thus, at the moment we do not have enough data to interpret the apparent shortcomings of the flexible clathrate model as either the consequence of the intrinsic drawbacks of the model itself, or inadequacy of the simulation technique, or poor quality of the SPC/E water potential. Despite the uncertainties in the computed P-T properties of methane hydrate, the analysis of simulated structures can provide information on the relationship between total cage occupancy and the degree of sI hydrate crystalline structure preservation. Figure 4 presents snapshots obtained from simulations in the framework of flexible clathrate model for six different average occupancies. Some fragments of clathrate structure can be

2890 J. Phys. Chem. B, Vol. 111, No. 11, 2007 observed even for the lowest occupancy value (27%), though for occupancies smaller than ca. 40% the hydrate structure is strongly degraded. Since the true stability of such low-occupancy structures cannot be directly determined from our simulations, in Figure 3 the occupancy contours corresponding to cage occupancy of 40% and less are represented by dotted rather than solid lines. For the system shown in Figure 4a the number of methane molecules exceeds the number of cages available (76 molecules in 64 cages). Double occupancy of some cages leads to noticeable local distortions of the hydrate structure as compared to lower occupancy (see Figure 4b for comparison). Such a phenomenon as double occupancy in methane hydrate is not a frequently observed one. Furthermore, the 12 excessive methane molecules can be found in 512 cages as well as in the 51262 ones. As a result, the θS/θL ratio is equal to unity for this particular configuration. Since the 512 cages are definitely too small to hold two methane molecules, the clathrate network around the doubly occupied cages is distorted so seriously that these cages can no longer be treated as 512. In experimental studies this structural degradation would result in noticeable changes in the measured unit cell parameters of the gas hydrate. More detailed analysis of simulation results was hindered by the lack of algorithm for reliable identification of large (51262) and small (512) cavities for strongly distorted hydrate structures. Independent monitoring of occupancies θS and θL for small and large cages, respectively, was therefore highly problematic. Evaluation of partial occupancies for simulations in the framework of the rigid clathrate model was carried out by calculating the distances between methane molecules and centers of cavities and assigning each methane molecule to the nearest cage. The θS/θL occupancy ratio did tend to increase with increasing gas fugacity. For medium-pressure values, the occupancies of large and small cages were approximately equal (θS/θL ∼ 1); for lower pressure values, the large cages showed higher occupancies than the small ones, while at high pressure a slight trend toward preferential filling of small cages could be found for some temperatures. The maximum and minimum calculated θS/θL ratios were 1.135 and 0.870. The results of the partial occupancy analysis for the rigid clathrate model are in general agreement with theoretical predictions and experimental data.1,8,9,29-33 Some studies insist on higher occupancy values for the large cages than those obtained from GCMC simulations or indicate significantly lower filling of small cages; for example, the θS/ θL ratio can be as low as 0.73.40 Conclusions As one compares the results obtained for rigid and flexible clathrate models, the advantages and disadvantages of both approaches to the description of the hydrate structure become obvious. Under the conditions of our GCMC simulation the rigid vdWP-like model provides qualitative agreement with experimental data and theoretical predictions. However, at higher temperature it fails to account for the structural changes accompanying the melting of the hydrate. Therefore, a single model cannot provide satisfactory description for both the hydrate-to-water phase transition and the hydrate composition below the melting point. Possibly, a quasi-rigid model with limited mobility of water molecules can be developed for this purpose. The impact of Monte Carlo simulations results is not limited to methane hydrates. On the contrary, hydrates of other common gases should be considered to study the influence of temperature and pressure on the molecular-level properties of gas hydrates

Sizov and Piotrovskaya and the relationship between the gas content and hydrate structure stability. We plan to expand this work by considering sI carbon dioxide hydrate and sII hydrates. Acknowledgment. Funding for this work came from INTAS Grant No. 03-51-5537. The authors wish to thank Dr. Martin B. Sweatman and Prof. Vladimir R. Belosludov for advice and fruitful discussion. References and Notes (1) Sloan, E. D. Clathrate Hydrates of Natural Gases, 2nd ed.; Marcel Dekker: New York, 1998. (2) Koh, C. A. Chem. Soc. ReV. 2002, 31, 157. (3) van der Waals, J. H.; Platteeuw, J. C. AdV. Chem. Phys. 1959, 2, 1. (4) Parrish, W. R.; Prausnitz, J. M. Ind. Eng. Chem., Process Des. DeV. 1972, 11, 26. (5) Tanaka, H. J. Chem. Phys. 1994, 101, 10833. (6) Westacott, R. E.; Rodger, M. P. Chem. Phys. Lett. 1996, 262, 47. (7) Tanaka, H.; Nakatsuka, T.; Koga, K. J. Chem. Phys. 2004, 121, 5488. (8) Ripmeester, J. A.; Ratcliffe, C. I. J. Phys. Chem. 1988, 92, 337. (9) Sum, A. K.; Burruss, R. C.; Sloan, E. D. J. Phys. Chem. B 1997, 101, 7371. (10) Chialvo, A. A.; Houssa, M.; Cummings, P. T. J. Phys. Chem. B 2002, 106, 442. (11) Norman, G. E.; Filinov, V. S. High Temp. 1969, 7, 216. (12) van Klaveren, E. P.; Michels, J. P. J.; Schouten, J. A.; Klug, D. D.; Tse, J. S. J. Chem. Phys. 2002, 117, 6637. (13) Yezdimer, E. M.; Cummings, P. T.; Chialvo, A. A. J. Phys. Chem. A 2002, 106, 7982. (14) Schober, H.; Itoh, H.; Klapproth, A.; Chihaia, V.; Kuhs, W. F. Eur. Phys. J. E 2003, 12, 41. (15) Alavi, S.; Ripmeester, J. A.; Klug, D. D. J. Chem. Phys. 2005, 123, 024507. (16) Okano, Y.; Yasuoka, K. J. Chem. Phys. 2006, 124, 024510. (17) Alder, B. J.; Wainwright, T. E. J. Chem. Phys. 1959, 31, 459. (18) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. (19) Jorgensen, W. L.; Madura, J. D.; Swenson, C. J. J. Amer. Chem. Soc. 1984, 106, 6638. (20) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. J. Amer. Chem. Soc. 1996, 118, 11225. (21) Bryk, T.; Haymet, A. D. J. J. Chem. Phys. 2002, 117, 10258. (22) Bryk, T.; Haymet, A. D. J. Mol. Simul. 2004, 30, 131. (23) Galloway, T. J.; Ruska, W.; Chappelear, P. S.; Kobayashi, R. Ind. Eng. Chem. Fundam. 1970, 9, 237. (24) Sanz, E.; Vega, C.; Abascal, J. L. F.; MacDowell, L. G. Phys. ReV. Lett. 2004, 92, 255701. (25) Vega, C.; Sanz, E.; Abascal, J. L. F. J. Chem. Phys. 2005, 122, 114507. (26) Vega, C.; Abascal, J. L. F.; Sanz, E.; MacDowell, L. G.; McBride, C. J. Phys.: Cond. Matter 2005, 17, 3283. (27) Fernandez, R. G.; Abascal, J. L. F.; Vega, C. J. Chem. Phys. 2006, 124, 144506. (28) McBride, C.; Vega, C.; Sanz, E.; MacDowell, L. G.; Abascal, J. L. F. Mol. Phys. 2005, 103, 1. (29) Galloway, T. J.; Ruska, W.; Chappelear, P. S.; Kobayashi, R. Ind. Eng. Chem. Fundam. 1970, 9, 237. (30) de Roo, J. L.; Peters, C. J.; Lichtenthaler, R. N.; Diepen, G. A. M. AIChE J. 1983, 29, 651. (31) Handa, Y. P. J. Chem. Thermodyn. 1986, 18, 915. (32) Klauda, J. B.; Sandler, S. I. J. Phys. Chem. B 2002, 106, 5722. (33) Sun, R.; Duan, Z. Geochim. Cosmochim. Acta 2005, 69, 4411. (34) Seo, Y.-T.; Lee, H. Korean J. Chem. Eng. 2003, 20, 1085. (35) Circone, S.; Kirby, S. H.; Stern, L. A. J. Phys. Chem. B 2005, 109, 9468. (36) Makogon, T. Y.; Sloan, E. D. J. Chem. Eng. Data 1994, 39, 351. (37) Nakano, S.; Moritoki, M.; Ohgaki, K. J. Chem. Eng. Data 1999, 44, 254. (38) Jager, M. D.; Sloan, E. D. Fluid Phase Equil. 2001, 185, 89. (39) Cao, Z.; Tester, J. W.; Sparks, K. A.; Trout, B. L. J. Phys. Chem. B 2001, 105, 10950. (40) Klapproth, A.; Goreshnik, E.; Staykova, D.; Klein, H.; Kuhs, W. F. Can. J. Phys. 2003, 81, 503.