J. Phys. Chem. B 2010, 114, 9563–9571
9563
Crystal Growth Simulations of H2S Hydrate Shuai Liang and Peter G. Kusalik* Department of Chemistry, UniVersity of Calgary, 2500 UniVersity DriVe Northwest, Calgary, Alberta, Canada ReceiVed: March 22, 2010; ReVised Manuscript ReceiVed: May 20, 2010
In this paper, we report a molecular simulation study exploring the crystal growth behavior of H2S hydrates within two-phase (H2S hydrate crystal and H2S aqueous solution) and three-phase (H2S hydrate crystal, liquid H2S, and H2S aqueous solution) systems. The microscopic mechanisms of growth, as well as the interfacial properties during the heterogeneous crystal growth process, are probed. We find that the H2S hydrate can be grown at a higher rate than methane hydrates under comparable conditions (Vatamanu, J.; Kusalik, P. G. J. Phys. Chem. B 2006, 110, 15896). The three-phase simulations, which also allow us to identify the simulation conditions on the experimental phase diagram, demonstrate that the present models reasonably reproduce the phase behavior of this system. We find that the crystal interface has a strong affinity for water molecules. We observed a relatively low level of defects in the newly formed H2S hydrate crystal. 1. Introduction Gas hydrates are crystalline solids with a host-guest structure in which water forms cages around the gas molecules.1 Many types of gaseous molecules, such as methane, ethane, propane, isobutene, nitrogen, hydrogen sulfide, and carbon dioxide, can be guests in hydrates. Gas hydrates have attracted great interest in large part because of their relevance to the natural gas industry,1-3 as well as their potential impacts in energy storage,4-6 global warming,7-10 geology,11 and gas separations.12,13 Of the common components in natural gas, H2S is known to form hydrates at the lowest pressure and persists to the highest temperature.13-24 In addition, H2S has a significant effect on hydrate formation of gas mixtures. Previous studies have shown that H2S forms hydrate more readily than do other components of natural gas.14,15,18 Thus, the presence of H2S in natural gas may help initiate the formation of hydrates, which when formed are difficult to remove and may cause problems in the production, transportation, and processing of natural gas. To help address these practical problems, a better understanding of the molecular mechanisms of gas hydrate crystal formation is of particular importance. At present, computer simulation approaches such as molecular dynamic (MD) simulations are among the only methods that can provide direct molecular details of the crystal formation/melting process.25,26 Several MD simulation studies addressing the formation and dissociation kinetics of CH4 hydrates have been carried out.27-31 This group has successfully performed MD simulation studies of the crystal growth of both structure I (sI) and structure II (sII) CH4 hydrates32-34 under steady-state conditions. H2S and CH4 have similar molecular diameters (4.58 and 4.36 Å, respectively), and they are observed to form sI hydrates, containing two dodecahedra (512) and six tetrakaidecahedra (51262) cages per unit cell, where the guests occupy both the 512 and 51262 cages.1 Beyond these similarities, H2S molecules can interact with the H2O molecules more strongly since H2S is a polar molecule, and the solubility of H2S in water is about 72 times higher than that of CH4 molecules.1 Therefore, the * To whom correspondence should be addressed. E-mail: pkusalik@ ucalgary.ca.
growth behavior of the H2S hydrates is likely to be somewhat different from that of CH4 hydrates. In this paper we report our MD simulation results to elucidate the microscopic crystal growth behavior of H2S hydrates from two-phase (H2S hydrate crystal and H2S aqueous solution) and three-phase (H2S hydrate crystal, liquid H2S, and H2S aqueous solution) systems. These simulations indicate that the growth of H2S hydrates can be much faster than that of CH4 hydrate, partially due to the high solubility and mobility of H2S molecules in water. In addition, compared to CH4 hydrates, a much lower level of crystal defects was observed with H2S hydrates, which is consistent with the high stability observed experimentally1,14,15,18 for this crystal. This paper is organized as follows: Section 2 briefly describes the potential models for H2O and H2S molecules used in this paper; it also provides some details of the hydrate systems and our MD methodology. In section 3, we first compare two available potential models for H2S molecules and then discuss the results obtained from crystal growth simulations of H2S hydrate systems. Finally, our concluding remarks are given in section 4. 2. Methodology Throughout this work, the intermolecular potential for water is modeled using the TIP4P/2005 potential,35 which is a rigid four-site model consisting of one Lennard-Jones center and three fixed point charges. This potential has been shown to have significantly improved water properties with respect to the original TIP4P model35 and has been widely used in computer simulations of water and gas hydrate systems.36-39 For H2S molecules, two rigid four-site potential models were chosen and compared in this work. The first (which will be referred to as model I in the following) was proposed by Forester et al.40 and was demonstrated to provide a good description of a range of properties of the bulk liquid and solid H2S phases. Kristo´f and Liszi41 reparametrized this four-site model to reproduce the vapor-liquid coexistence properties of H2S (we refer to this as model II). As part of this study, we have compared the lattice parameter, potential energy, hydrate melting point, and crystal growth behavior using these two potential models and found no measurable difference between them (as
10.1021/jp102584d 2010 American Chemical Society Published on Web 07/02/2010
9564
J. Phys. Chem. B, Vol. 114, No. 29, 2010
Liang and Kusalik
Figure 1. (a) Schematic representation of the MD setup used for the crystal growth simulations. The movement of the temperature pulse is in the z-direction. Periodic boundary conditions are applied in all three Cartesian directions. (b, c) Configurations from two-phase and three-phase systems, respectively. The H2O molecules have been labeled as solidlike (green spheres) and liquidlike (red spheres) as discussed in the text. The yellow spheres indicate H2S molecules.
will be detailed in section 3.1). Since under our simulation conditions no vapor phase of H2S is expected, we have chosen to utilize model I in our production runs of H2S hydrate crystal growth. The interactions between H2S and H2O molecules were calculated as the usual sums of Coulomb and Lennard-Jones potentials. The cross terms for the Lennard-Jones parameters were determined using Lorentz-Berthelot mixing rules.26 The MD simulations were performed following our methodology described previously.32,33,42 The equations of motion were integrated with a fourth-order predictor-corrector algorithm and a time step of 1 fs. A spherical cutoff of 11 Å was utilized for the short-ranged forces, while the electrostatic interactions were evaluated with the smooth particle mesh Ewald method.43 Periodic boundary conditions were imposed in all three Cartesian directions. Further details of the MD methodology we have employed can be found elsewhere.32,33,42 An initial crystal structure of sI H2S hydrate was constructed with lattice parameter a ) 12 Å.44-46 In this structure, H2S molecules occupied all the small (512) and large (51262) cages of the crystal. The original hydrate system was composed of 1840 H2O molecules and 320 H2S molecules. The system cross section was 2 × 2 unit cells (24 Å × 24 Å), and the simulation cell was roughly 120 Å long. A schematic setup of the present system is shown in Figure 1a. A typical simulation run involved the following steps: (1) The crystal structure was first equilibrated for about 5 ns at a temperature of 285 K and pressure of 100 atm. (2) A high but narrow temperature pulse (500 K) was then applied briefly to initiate localized melting of the crystal. (3) Subsequently, about 50% of the crystal was melted using a step temperature profile (with temperatures of 275 and 300 K). (4) During the crystal growth simulation that followed, a narrow temperature pulse of 335 K was moved through the system (along the z-direction; see Figure 1) to melt the crystal at one interface while a constant undercooling of 275 K was maintained in the remainder of the system to induce crystallization at the other interface. If the amount of crystal melted at one interface is equivalent to the amount of crystal formed at the second interface, then steadystate conditions can be achieved. All the results reported here are in the NPT ensemble with a pressure of 100 atm.
In this study, two systems, one two-phase system and one three-phase system, were prepared as the starting points for steady-state crystal growth simulations. For the two-phase system, 25% of the H2S molecules were randomly removed from the original fully occupied hydrate crystal. Then by melting about 50% of the crystal, a system consisting of H2S hydrate crystal and H2S aqueous solution was prepared as shown in Figure 1b. Although the resulting H2S solution is highly supersaturated, bubble formation was not observed during our simulations (on a scale of 20-30 ns). For the three-phase system, about 50% of the fully occupied crystal was first melted, and then a simulation of about 5 ns was performed at 285 K and 100 atm to allow formation of a separate H2S liquid phase. Since this temperature and pressure are close to the melting point of the hydrate crystal (as discussed in section 3.1), no significant crystal growth or melting was observed during this equilibration simulation. The resulting three-phase system, consisting of hydrate crystal, liquid H2S, and H2S aqueous solution, is shown in Figure 1c. While various rates for the movement of the temperature pulse were explored, the results reported here are with constant moving rates of 2.5 Å/ns for the two-phase system and 0.3 Å/ns for the three-phase system, which reasonably correspond to the observed growth rates (as reported below). The simulation times of two-phase and threephase crystal growth runs were about 25 and 70 ns, respectively. As part of our previous work,32,33,42 we have tested the influence of the position of the temperature pulse on the growth behavior of the crystal. We found no evidence of correlations as long as the temperature pulse was restricted to a narrow region well away (i.e., more than 25 Å) from the growing interface. The labeling of molecules as liquidlike/solidlike was performed by utilizing a translational mean-squared displacement criterion as described previously.42,47 On the basis of our experiences,32-34,42,48,49 solidlike and liquidlike character can be determined from a molecular history over roughly 10-100 ps. The segment times used here are 10 ps in two-phase simulations and 50 ps for three-phase simulations where this choice acknowledges that the crystal grows at a much lower rate within the three-phase system. All of the configurations presented below utilize averaged molecular coordinates48 over these time windows.
Crystal Growth Simulations of H2S Hydrate
J. Phys. Chem. B, Vol. 114, No. 29, 2010 9565
Figure 2. Melting temperature determination for model I (a) and model II (b) for H2S hydrates, as described in the text. Each temperature point includes results from five statistically independent system runs. ∆E is defined as the energy difference between the initial and final total energies of the system (averaged over 5 ps). The error bars represent the standard deviation of the energy changes observed.
3. Results and Discussion 3.1. Comparison of Two H2S Potential Models. Our previous results suggest that the maximum growth rate of hydrate is obtained with a temperature 10-20 K below the melting temperature.32,33 Thus, it is useful to find the melting temperature (Tm) of our model H2S hydrate systems. Since melting temperatures can be rather sensitive to the interaction potential, we will also use these melting temperatures to help validate the two H2S potential models (I and II). For these purposes, we track the total energy of the system to determine an estimate for Tm. Specifically, after a 5 ns equilibration period, we record the energy change of the system at a constant desired temperature, where no heat pulse has been applied. If this temperature is higher than Tm, the crystal is expected to melt and the total energy of the system will increase. Otherwise, if the temperature is lower than Tm, the H2S aqueous solution is
expected to crystallize and thus the total energy of the system will decrease. If we run the simulation at Tm, the energy of the system should maintain roughly a constant value. By sampling the energy changes for a series of temperatures, Tm can be estimated. To avoid the uncertainty introduced by the H2S concentration change or phase separation during the measuring interval, we have utilized a three-phase system in this measurement. The energy changes averaged over five statistically independent 1 ns long runs at seven temperatures from 277 to 295 K were measured with both model I and model II, as shown in parts a and b, respectively, of Figure 2. A linear fit of the ∆E dependence on T gave estimates of Tm ) 290.9 ( 5.2 K for model I and Tm ) 288.7 ( 7.6 K for model II. No measurable difference between the values from the two models was observed, which are also reasonably consistent with the reported
9566
J. Phys. Chem. B, Vol. 114, No. 29, 2010
Liang and Kusalik
Figure 3. Molecular configurations (a) at the starting point and (b) after 25 ns of steady-state crystal growth from the two-phase system. The molecules are colored as in Figure 1.
experimental melting temperatures for H2S hydrate (about 30 K above the ice point for the water model).15,20,35 Although the dependence of ∆E on T will not be linear in general, we believe that the linear fitting used here should be a reasonable approximation over this 18 K temperature window. The large fluctuations of ∆E arise from the relatively short time scale of these simulation runs, as we will discuss below. The reasonable agreement of the measured melting temperature with experiment confirms the quality of the present potential models. We also compared the lattice parameters of cubic sI hydrate crystals using the two H2S potential models. These simulations were performed at a temperature of 125 K and pressure of 10 atm with a fully occupied crystal. The resulting lattice parameters for both models were in the range of 11.82 ( 0.01 Å, which is close to the experimental value of 11.96 Å under similar conditions.45 These lattice parameters are found to be invariant to the chosen H2S potential model. The potential energies for the two systems using either of the two different H2S potential models were both about 54.9 ( 0.1 kJ/mol for these conditions. Further comparison of the two models was performed by investigating the crystal growth behavior of the hydrates (data not shown). We observed no measurable difference in growth rates, crystal defects, or interfacial character within either two-phase or three-phase simulations. As mentioned above, since under our simulation conditions only liquid H2S and no vapor was expected, we chose to utilize model I in our production runs of hydrate crystal growth. 3.2. Crystal Growth within a Two-Phase System. Starting from the two-phase H2S/H2O system shown in Figure 3a, crystal growth was investigated at a temperature of 275 K. After 25 ns, about 65 Å of new crystal was found to have formed, and the final configuration is presented in Figure 3b. As can be seen in Figure 3, the original hydrate crystal melted completely and the original liquid crystallized completely. Figure 4a shows the position of the growing crystal interface as a function of the simulation time. Here the position of the interface was determined by using the z-profile of mean-square displacements (MSDs) of H2O molecules, following a recently developed interface detection approach.47 We note that the position of the interface is defined as a simple 1-D quantity (averaged over x-y-planes) instead of a more realistic 3-D quantity, but it should provide a reasonable estimate over the long simulation times reported here. The numerical derivative of the time dependence of the interfacial position gives instantaneous growth rates, as shown in Figure 4b. The mean growth rate is 2.62 Å/ns, which is significantly higher than that determined for methane hydrates under comparable conditions.33 We note that the level of supersaturation of the gas solution is much lower here. The relatively high growth rate observed for
this H2S hydrate indicates the importance of kinetic factors, including mass transfer, in hydrate growth. The large fluctuations of the instantaneous growth rates, while being somewhat sensitive to the method of detection, are indicative of the intrinsically stochastic (ordering/disordering) nature of the crystal growth process. The distribution of the instantaneous growth rates is presented in Figure 4c. A well-defined Gaussianshaped distribution is apparent, which further points to the stochastic behavior of the crystal growth. While the results presented in Figure 4 utilize MSD profiles of H2O molecules, we have also performed the same measurements with the H2S molecules and have obtained entirely consistent results (see the Supporting Information). A notable character of the growth kinetics shown in Figure 4a is that the crystal growth exhibits a mixed pattern of steplike and continuous growth. For example, from about 5.5 to 8.5 ns, the profile exhibits a relatively flat period, which suggests that the crystal growth was stagnant for nearly 3.0 ns. Immediately afterward, an almost vertical increase of about 5 Å was observed (see Figure 4a, around 8.5-9.0 ns), suggesting a rapid advance of the growing hydrate interface. The corresponding molecular configurations are shown in Figure 5a-c. We can see the crystal/ liquid interface maintained a roughly planar structure, indicative of a layer-by-layer growth mechanism. Such a mechanism, which is promoted in smaller systems by finite size effects, has been observed in the basal interface of small ice/water systems.50-53 With further inspection of molecular configurations during periods of more continuous growth, we find that the newly grown crystal exhibits clear microfaceted structuring at its interface. Two examples of such configurations, at 15 and 22 ns, respectively, are given in parts a and b of Figure 6. We note the orientation of these microfacets is apparently not uniquely defined; i.e., at any particular instant time, one can find a varying detailed topology as a result of the inherent fluctuating nature of the interface. These observations of microfacets are consistent with the behavior found in atomic systems49 and sII hydrate systems.32 The crystal growth in this relatively small (2 × 2) system exhibits both layer-by-layer and microfaceted mechanisms. We remark that, on the basis of experience, when using systems with larger cross-sectional areas, one can expect microfaceted growth to be the dominant process. We observed that the H2S hydrate crystal growth has far fewer vacancy defects than exhibited by the growth of methane hydrates.32-34 The cage occupancy of the newly formed crystal is about 92%, leaving roughly 8% of the water cages unfilled. The H2S molecules appear to occupy equally the two types of cages in the sI hydrate crystal, consistent with experimental observations.54 Usually, empty cages are surrounded by filled
Crystal Growth Simulations of H2S Hydrate
J. Phys. Chem. B, Vol. 114, No. 29, 2010 9567
Figure 4. Crystal growth kinetics from the two-phase system: (a) position of the growing interface, determined using the MSD profiles of H2O molecules, as a function of the simulation time, (b) instantaneous growth rate as given by the numerical derivation of the interface position vs time curve shown in (a), (c) distribution of the instantaneous growth rate, where the red line is the Gaussian fit with a center of 2.62 Å/ns.
9568
J. Phys. Chem. B, Vol. 114, No. 29, 2010
Figure 5. Three cropped configurations from the MD trajectory of the two-phase crystal growth with molecules colored as in Figure 1. The vertical yellow line indicates the position of the interface at 5.6 ns. We can see that the crystal growth was stagnant for nearly 3.0 ns between 5.6 ns (a) and 8.3 ns (b), followed by rapid growth between 8.3 ns (b) and 8.7 ns (c).
Figure 6. Cropped system configurations (a) at 15 ns and (b) at 22 ns, showing molecular structures at the solid/liquid interface from twophase crystal growth. The molecules are colored as in Figure 1. The yellow lines emphasize the microfacets exposed to the liquid.
Figure 7. Two examples of empty cages observed from a two-phase crystal growth simulation. Two neighboring empty 51262 cages are shown in (a), and one empty 512 cage is shown in (b).
cages, which seems to be necessary for their stabilization. Interestingly, two neighboring empty 51262 cages were found in this simulation, as can be seen in Figure 7a. Figure 7b gives an example of an empty small 512 cage. No other defects, such as cages containing more than one gas or water molecules observed in methane hydrates,32-34 were observed here. This low level of defects is consistent with the relatively high stability of H2S hydrates measured experimentally.1,14,15,18 It is interesting to point out that the structural transition involving 51263 cages observed in methane hydrates33 appears also to occur but with a much lower frequency. This behavior will be explored in greater detail in future work. 3.3. Crystal Growth within a Three-Phase System. The crystal growth kinetics from the three-phase system (see Figure 8) was also investigated. We find that in such as system the H2S bubble can become an effective barrier to the transport of
Liang and Kusalik
Figure 8. Four cropped molecular configurations from the three-phase crystal growth simulation: (a) initial configuration, (b) after 50 ns, (c) after 55 ns, (d) after 70 ns. The molecules are colored as in Figure 1. Note that, in (c), a water bridge has formed across the liquid H2S bubble.
H2O molecules. As can be seen in Figure 8b, the bubble remained at basically the same position during the first 50 ns of this crystal growth simulation, which implies almost no net water molecules have crossed the bubble during this interval. Interestingly, once the hydrate crystal has grown close to this liquid H2S bubble, a water bridge forms, temporarily allowing rapid water transport from the right (solution) side of the bubble to its left (growing interface) side. We have found that the liquid H2S/H2O interface is generally not planar, and molecular chains composed of several H2O molecules are occasionally observed to form, providing a “bridge” for H2O molecules to cross the liquid H2S bubble. When both sides of the H2S bubble are aqueous solution, this bridge soon breaks. With the hydrate crystal close to the H2S bubble, the formation of such a bridge appears to allow a large number of H2O molecules to cross the H2S bubble (Figure 8c), providing a means of “fast” water transport. As can be seen in Figure 8d, about 12 Å of liquid H2O accumulates between the crystal and H2S bubble as a result of this fast transport. An animation of this process can be found in the Supporting Information. This fast transport process suggests that the contact of hydrate-H2S-H2O is less favorable than the hydrate-H2O-H2S contact. Figure 9 shows the growing interface position in the threephase system as a function of time. Similar to the crystal growth from the two-phase system, this growth kinetics curve also exhibits a mixture of steplike and continuous behavior. However, in Figure 9 the periods of apparent growth stagnation are not as clear as was the case in the two-phase system (at the same temperature and pressure; see Figure 4), and most of the time the crystal grows in a reasonably continuous manner. This is perhaps because of the lower concentration of H2S in solution. Since empty “pockets” (half cages) at the crystal/liquid interface have a limited chance of capturing a H2S molecule, only those sites with H2S molecules are likely to form complete cages. The interface consequently exhibits a more complex topology (some examples can be found in the Supporting Information). With the lower concentration of H2S molecules in the solution, we get a much slower growth rate of 0.39 Å/ns, which is almost an order of magnitude lower than was achieved in the twophase system. We remark that the H2S solution composition is roughly 3 times higher in the two-phase system.
Crystal Growth Simulations of H2S Hydrate
J. Phys. Chem. B, Vol. 114, No. 29, 2010 9569
Figure 9. Position of the growing interface in the three-phase crystal growth simulation as a function of the simulation time, where the position was determined from the MSD profiles of H2O molecules as discussed in the text.
Figure 10. Cropped molecular configurations from three-phase crystal growth, with molecules colored as in Figure 1. The two vertical yellow lines indicate that the newly formed crystal at 47 ns (a) was remelted at 52 ns (b).
On several occasions, instead of remaining constant or advancing, the interfacial position apparently recedes, as can be seen clearly in Figure 9 (around 30 and 50 ns). This suggested that, at these times, the crystal had melted back under these growing conditions. Figure 10 shows the corresponding molecular configurations at 47 and 52 ns, from which we can see that about 3-4 Å of the newly formed crystal apparent at 47 ns was melted 5 ns later. We note that the distribution of crystal growth rates from this system, where the H2S solution is not highly supersaturated, is rather similar to that pictured in Figure 4c. Similar crystal melting phenomena were not observed in the case of the two-phase H2S hydrate system (see section 3.2), where the high level of supersaturation of H2S provides a very strong driving force for crystallization and hence the average growth rate is considerably faster and the likelihood of observing a melt-back fluctuation is significantly lower. As in the two-phase system, a very low level of vacancy and other crystal defects was observed in the crystal grown in this three-phase system. Interestingly, we did observe that H2O molecules that make up a cage can temporarily enter empty cages and become trapped in a distorted cage, but such structures do not survive. A more detailed analysis of this phenomenon is currently under way. Finally, we emphasize an advantage of the three-phase system simulation. We note that the (mechanical) pressure in the MD simulation and the (gas) pressure in a real system are not the same quantity in the two-phase system simulation, since the
third (vapor) phase is not present. In the absence of this vapor phase, the connection between pressure and the composition in model simulation results can be problematic. Three-phase simulations of hydrate nucleation or growth have been reported in a few studies27,29,55,56 for methane hydrate systems. The methane/water interface, with a relatively high interfacial tension,57 effectively impedes gas molecules from penetrating into the aqueous phase. Even under a pressure as high as 50 MPa, the equilibration time can still be on microsecond scales.55 Thus, direct simulations of the crystal growth and/or nucleation processes of methane hydrates in three-phase systems are computationally expensive to realize. The relatively high solubility and mobility of H2S molecules in water give them an advantage when three-phase simulations are performed. An evaluation of the H2S concentration during the simulation gives a composition yH2S ) 0.031(5) for the solution near the crystal growing interface and yH2S ) 0.038(4) for the solution near the melting interface, where the values in parentheses indicate standard errors. Experimental data20 predict that at the quadruple (water, hydrate, vapor, liquid H2S) point (see Figure 11) the H2S composition in the aqueous solution is 0.0323. Assuming that the solubility of H2S in water is essentially constant with increased pressure above this point (up to 20 MPa),58,59 we see that the present model system roughly reproduces the solubility of the real system. We can identify our (nonequilibrium) crystal growth simulation conditions roughly on the experimental phase diagram, as shown in Figure 11. The applied temperature (275 K) has been shifted according to the measured melting temperature (about 290 K; see section 3.1) of the model system. The appearance of the point in the hydrate-liquid H2S stable region of the experimental phase diagram represents an important link between our MD simulations and the real systems. 4. Conclusions In this paper, we reported the results from MD simulations of the crystal growth of H2S hydrates within two-phase and three-phase systems in the NPT ensemble. In the case of the two-phase system, a growth rate of about 2.62 Å/ns was obtained at 275 K and 100 atm, which is significantly higher than that determined for methane hydrates grown from a solution with a
9570
J. Phys. Chem. B, Vol. 114, No. 29, 2010
Liang and Kusalik will explore the role of the gas molecules, the size effects, and the influence of the degree of supersaturation and of the temperature on crystal growth behavior. Acknowledgment. We are grateful for the financial support of the Natural Sciences and Engineering Research Council of Canada. We also acknowledge computational resources made available via WestGrid (www.westgrid.ca) and the University of Calgary.
Figure 11. Illustration of the nonequilibrium simulation conditions as transcribed onto the experimental pressure-temperature phase diagram of a H2S hydrate system. The blue star roughly indicates the conditions of our three-phase crystal growth. Note the temperature has been shifted by 15 K according to the corresponding estimate of the melting temperature of the model system (see the text). In the diagrams, the symbols Lw, H, V, and LH2S represent liquid water, hydrate, vapor, and liquid H2S, respectively. Note that only the two quadruple points of Q1 and Q2 are accurate experimental values taken from ref 1.
much higher level of supersaturation, indicating the importance of kinetic factors such as simple mass transfer in hydrate growth. The crystal growth in this relatively small (2 × 2) system exhibits both layer-by-layer (which can be seen as a consequence of finite size effects) and microfaceted mechanisms. In the case of the three-phase system, the crystal was grown from a solution with much lower H2S concentrations at 275 K and 100 atm. Under such conditions, the crystal growth rate is about 0.39 Å/ns, which is almost 1 order of magnitude less than that from the two-phase system. Occasionally, and correlated with the local compositional changes of the solution near the growing interface, the newly formed crystal was observed to melt back. The behavior of the three-phase simulation confirms that the hydrate crystal prefers to be “wet” by the aqueous phase. In addition, the newly formed crystals from both two-phase and three-phase systems were found to contain a lower level of vacancy and other crystal defects than exhibited by the growth of methane hydrates.32-34 This low level of defects is relatable to the relatively high stability of H2S hydrates observed experimentally. The three-phase simulation allows us to roughly identify the simulation conditions on the experimental phase diagram. An interesting experimental observation reported recently60 is that H2S can rapidly replace dimethyl ether (DME) from sII hydrates, which requires rapid transport of gas molecules through the hydrate. We found no specific evidence of H2S transport through the crystal in our simulations. This may be because the time scale of the current simulations is simply too short. Nonetheless, we did observe water transport in the crystal facilitated by vacant cages. This perhaps points to a possible limitation of the present model system. H2S molecules cannot be further polarized when they interact strongly with water molecules, and thus, it is harder for H2S to compete for water H-bonds in our simulations (relative to the real system). Careful examination and analyses for possible molecular transports through the solid aimed at providing molecular mechanisms explaining the experimental observations is ongoing. We believe that the microscopic insights obtained here are of fundamental interest and will help advance the understanding of the crystal growth behavior of gas hydrates. Future work
Supporting Information Available: Interface position data determined from H2S molecules within the two-phase simulation, some examples of the interface topologies within the threephase simulation, and an animation of the fast water transport across the liquid H2S bubble within the three-phase simulation. This material is available free of charge via the Internet at http:// pubs.acs.org. References and Notes (1) Sloan, E. D.; Koh, C. A. Clathrate Hydrates of Natural Gases; CRC Press: Boca Raton, FL, 2007. (2) Sloan, E. D. Nature 2003, 426, 353. (3) Boswell, R. Science 2009, 325, 957. (4) Strobel, T. A.; Kim, Y.; Andrews, G. S.; Ferrell, J. R.; Koh, C. A.; Herring, A. M.; Sloan, E. D. J. Am. Chem. Soc. 2008, 130, 14975. (5) Lee, H.; Lee, J. W.; Kim, D. Y.; Park, J.; Seo, Y. T.; Zeng, H.; Moudrakovski, I. L.; Ratcliffe, C. I.; Ripmeester, J. A. Nature 2005, 434, 743. (6) Florusse, L. J.; Peters, C. J.; Schoonman, J.; Hester, K. C.; Koh, C. A.; Dec, S. F.; Marsh, K. N.; Sloan, E. D. Science 2004, 306, 469. (7) Archer, D. Biogeosciences 2007, 4, 521. (8) Bohannon, J. Science 2008, 319, 1753. (9) McElwain, J. C.; Wade-Murphy, J.; Hesselbo, S. P. Nature 2005, 435, 479. (10) Svensen, H.; Planke, S.; Malthe-Sorenssen, A.; Jamtveit, B.; Myklebust, R.; Eidem, T. R.; Rey, S. S. Nature 2004, 429, 542. (11) Kleinberg, R. L.; Flaum, C.; Griffin, D. D.; Brewer, P. G.; Malby, G. E.; Peltzer, E. T.; Yesinowski, J. P. J. Geophys. Res., [Solid Earth] 2003, 108, 2508. (12) Kang, S. P.; Lee, H. EnViron. Sci. Technol. 2000, 34, 4397. (13) Yamamoto, Y.; Komai, T.; Yoon, J. H.; Kang, S. P.; Okita, S.; Kawamura, T. J. Chem. Eng. Jpn. 2003, 36, 971. (14) Carroll, J. J. Natural Gas Hydrates: A Guide for Engineers; Gulf Professional Publishing: Amsterdam, 2003. (15) Carroll, J. J.; Mather, A. E. Can. J. Chem. Eng. 1991, 69, 1206. (16) Wohler, F. Ann. Chem. Pharm. 1840, 33, 125. (17) Schicks, J. M.; Lu, H.; Ripmeester, J. A.; Ziemann, M. The Characteristics of Gas Hydrates Formed From H2S and CH4 Under Various Conditions. Proceedings of the 6th International Conference on Gas Hydrates (ICGH 2008), Vancouver, British Columbia, Canada, July 6-10, 2008. (18) Carroll, J. J. Presented at the GPA Europe Spring Meeting, Dublin, 2004. (19) Davidson, D. W.; Garg, S. K.; Gough, S. R.; Hawkins, R. E.; Ripmeester, J. A. Can. J. Chem. 1977, 55, 3641. (20) Selleck, F. T.; Carmichael, L. T.; Sage, B. H. Ind. Eng. Chem. 1952, 44, 2219. (21) Wright, R. H.; Maass, O. Can. J. Res. 1932, 6, 94–101. (22) Ng, H. J.; Robinson, D. B. Fluid Phase Equilib. 1985, 21, 145. (23) Cady, G. H. J. Phys. Chem. 1981, 85, 3225. (24) Ng, H. J.; Robinson, D. B. Fluid Phase Equilib. 1985, 19, 21. (25) Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Applications; Academic Press: San Diego, 2002. (26) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: Oxford, U.K., 1989. (27) Hawtin, R. W.; Quigley, D.; Rodger, P. M. Phys. Chem. Chem. Phys. 2008, 10, 4853. (28) Moon, C.; Taylor, P. C.; Rodger, P. M. J. Am. Chem. Soc. 2003, 125, 4706. (29) Nada, H. J. Phys. Chem. B 2006, 110, 16526. (30) English, N. J.; Johnson, J. K.; Taylor, C. E. J. Chem. Phys. 2005, 123, 244503. (31) Baez, L. A.; Clancy, P. Ann. N. Y. Acad. Sci. 1994, 715, 177. (32) Vatamanu, J.; Kusalik, P. G. J. Phys. Chem. B 2008, 112, 2399. (33) Vatamanu, J.; Kusalik, P. G. J. Phys. Chem. B 2006, 110, 15896. (34) Vatamanu, J.; Kusalik, P. G. J. Am. Chem. Soc. 2006, 128, 15588. (35) Abascal, J. L. F.; Vega, C. J. Chem. Phys. 2005, 123, 234505.
Crystal Growth Simulations of H2S Hydrate (36) Aragones, J. L.; Conde, M. M.; Noya, E. G.; Vega, C. Phys. Chem. Chem. Phys. 2009, 11, 543. (37) Dyer, P. J.; Docherty, H.; Cummings, P. T. J. Chem. Phys. 2008, 129, 024508. (38) Ghoufi, A.; Goujon, F.; Lachet, V.; Malfreyt, P. J. Chem. Phys. 2008, 128, 154716. (39) Tribello, G. A.; Slater, B. J. Chem. Phys. 2009, 131, 024703. (40) Forester, T. R.; McDonald, I. R.; Klein, M. L. Chem. Phys. 1989, 129, 225. (41) Kristof, T.; Liszi, J. J. Phys. Chem. B 1997, 101, 5480. (42) Vatamanu, J.; Kusalik, P. G. J. Chem. Phys. 2007, 126, 124703. (43) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. J. Chem. Phys. 1995, 103, 8577. (44) Jeffrey, G. A. J. Inclusion Phenom. 1984, 1, 211. (45) Davidson, D. W.; Garg, S. K.; Gough, S. R.; Handa, Y. P.; Ratcliffe, C. I.; Tse, J. S.; Ripmeester, J. A. J. Inclusion Phenom. 1984, 2, 231. (46) Gutt, C.; Asmussen, B.; Press, W.; Johnson, M. R.; Handa, Y. P.; Tse, J. S. J. Chem. Phys. 2000, 113, 4713. (47) Liang, S.; Kusalik, P. G. Chem. Phys. Lett., in press. DOI: 10.1016/ j.cplett.2010.05.088. (48) Kusalik, P. G.; Gillis, K.; Vatamanu, J. 21st International Symposium on High Performance Computing Systems and Applications
J. Phys. Chem. B, Vol. 114, No. 29, 2010 9571 (HPCS’07), Saskatoon, Canada, 2007; Springer-Verlag Berlin: Berlin, Germany; pp 51-58. (49) Vatamanu, J.; Kusalik, P. G. Phys. ReV. B 2007, 76, 035431. (50) Nada, H.; Furukawa, Y. J. Cryst. Growth 1996, 169, 587. (51) Nada, H.; Furukawa, Y. J. Phys. Chem. B 1997, 101, 6163. (52) Nada, H.; Furukawa, Y. Jpn. J. Appl. Phys., Part 1 1995, 34, 583. (53) Nada, H.; Furukawa, Y. Surf. Sci. 2000, 446, 1. (54) Cady, G. H. J. Phys. Chem. 1983, 87, 4437. (55) Walsh, M. R.; Koh, C. A.; Sloan, E. D.; Sum, A. K.; Wu, D. T. Science 2009, 326, 1095. (56) Zhang, J. F.; Hawtin, R. W.; Yang, Y.; Nakagava, E.; Rivero, M.; Choi, S. K.; Rodger, P. M. J. Phys. Chem. B 2008, 112, 10608. (57) Schmidt, K. A. G.; Folas, G. K.; Kvamme, B. Fluid Phase Equilib. 2007, 261, 230. (58) Koschel, D.; Coxam, J. Y.; Majer, V. Ind. Eng. Chem. Res. 2007, 46, 1421. (59) Gillespie, P. C.; Wilson, G. M. Vapor-Liquid and Liquid-Liquid Equilibria: Water-Methane, Water-Carbon Dioxide, Water-Hydrogen Sulfide, Water-n-Pentane; Reserach Report RR-48, Project 758-B77; Gas Processors Association: Tulsa, OK, 1982. (60) Buch, V.; Devline, J. P.; Monreal, I. A.; Jagoda-Cwiklik, B.; UrasAytemiz, N.; Cwiklik, L. Phys. Chem. Chem. Phys. 2009, 11, 10245.
JP102584D
