Molecular Dynamics Study of Carbon Dioxide Hydrate Dissociation

ARTICLE pubs.acs.org/JPCA

Molecular Dynamics Study of Carbon Dioxide Hydrate Dissociation Sapna Sarupria and Pablo G. Debenedetti* Department of Chemical and Biological Engineering, Princeton University, Princeton, New Jersey 08544, United States ABSTRACT: We present results from a molecular dynamics study of the dissociation behavior of carbon dioxide (CO2) hydrates. We explore the effects of hydrate occupancy and temperature on the rate of hydrate dissociation. We quantify the rate of dissociation by tracking CO2 release into the liquid water phase as well as the velocity of the hydrateliquid water interface. Our results show that the rate of dissociation is dependent on the fractional occupancy of each cage type and cannot be described simply in terms of overall hydrate occupancy. Specifically, we find that hydrates with similar overall occupancy differ in their dissociation behavior depending on whether the small or large cages are empty. In addition, individual cages behave differently depending on their surrounding environment. For the same overall occupancy, filled small and large cages dissociate faster in the presence of empty large cages than when empty small cages are present. Therefore, hydrate dissociation is a collective phenomenon that cannot be described by focusing solely on individual cage behavior.

I. INTRODUCTION Clathrate hydrates are crystalline solids in which guest molecules reside in polyhedral water cages stabilized by hydrogen bonds. The guest molecules range from gases such as methane, ethane, propane, isobutane, carbon dioxide, nitrogen, and hydrogen to heavier compounds like cyclopentane and cyclobutanone. Hydrates are important in hydrocarbon processing1,2 and could play a major role as an energy source.35 The use of hydrates for hydrogen storage,68 water desalination,911 gas storage and transportation,1216 and separation of mixtures1719 is also being investigated.20 Hydrates also feature in the strategies being considered for carbon dioxide (CO2) sequestration. It has been suggested that CO2 can be stored as a solid in the form of clathrate hydrates or as a pool of liquid below a cap of its hydrate.21,22 It has also been suggested that the economic viability of CO2 storage as hydrates can be improved by combining it with the recovery of methane from its hydrate reserves, such that methane is replaced with CO2 in the hydrate.14 Engineering any technology along these lines requires a fundamental understanding of both the thermodynamics and the kinetics of CO2 hydrate formation and dissociation. The study presented here focuses on the dissociation behavior of carbon dioxide hydrates. This work is part of our ongoing broad investigation of the thermodynamics and kinetics of CO2 hydrate formation and dissociation. Molecular-based computer simulations are a powerful tool to study the microscopic structure of clathrate hydrates and its correlation with macroscopic properties. Since the initial investigation of methane hydrates by Tester and co-workers23 and Tse et al.,24,25 several simulation studies have been reported. While the majority of these investigations are on methane hydrates, other systems have also been studied. The various aspects explored through simulations include the free energy of replacing r 2011 American Chemical Society

methane with other guest molecules,26,27 the effect of helper gas molecules on the stability of hydrates,2831 nucleation,3237 formation3844 and dissociation4552 of gas hydrates, and the effect of inhibitors on hydrate formation and stability.5358 These studies have provided valuable insights into the thermodynamic stability of hydrates, and their formation and dissociation behavior in bulk and at interfaces, and have also contributed to the development of better predictive models for hydrate stability.59,60 Motivated by this body of work, we use molecular dynamics (MD) simulations to investigate various aspects of CO2 hydrate dissociation. The kinetics of hydrate dissociation is quantified in terms of the rate of CO2 release into bulk water and the velocity of the hydratewater interface. An important conclusion from our study is that the rate of hydrate dissociation is dependent not only on the overall occupancy of the hydrate, but also on the cage-specific occupancy. This subtlety has not been appreciated previously. In addition, this dependence cannot be explained by the dissociation behavior of individual cages because the latter is affected by the surrounding environment. These results suggest that hydrate dissociation is a collective phenomenon. The paper is structured as follows: The setup of the simulation system, details of the MD simulations, and the method used to classify water molecules as hydrate- or liquid-like are described in section II. Results from our study are discussed in section III. Conclusions and suggestions for further studies are presented in section IV. Special Issue: Victoria Buch Memorial Received: November 14, 2010 Revised: March 6, 2011 Published: March 23, 2011 6102

dx.doi.org/10.1021/jp110868t | J. Phys. Chem. A 2011, 115, 6102–6111

The Journal of Physical Chemistry A

ARTICLE

Figure 1. Snapshot of the hydratewater system. (a) Oxygen and hydrogen atoms are shown in red and white, respectively. (b) Classification of water molecules as hydrate-, ice-, or liquid-like for the snapshot shown in (a). Hydrate-like water molecules are colored in red and liquid-like molecules are blue. No ice-like water molecules were identified in this snapshot. The z-dimension corresponds to the horizontal direction, normal to the hydratewater interface.

II. METHODS A. Setup of the Simulation System. Simulations of periodically replicated crystals in the absence of a solidliquid interface require superheating to melt the ordered solid, and this hinders the accurate estimation of melting points. To overcome this, Fernandez et al.61 suggested simulating systems in which solid and liquid phases are in contact. They demonstrated that the melting point of ice can be reproduced with good accuracy using this method. Subsequently, Myshakin et al.48 used a similar methodology to study the melting of the methane hydrate system. In our work, we have used this approach to study the melting behavior of CO2 hydrates. The initial hydrate structure was generated from the coordinates provided in Hollander and Jeffrey62 (also reported in ref 2). The oxygen coordinates were replicated four times in each dimension to give a 4.8 4.8 4.8 nm3 cubic box of hydrate. Hydrogens were assigned random orientations while maintaining the geometry of the water molecules. The energy of the resulting configuration was minimized keeping the oxygens frozen, using the Steepest Descent algorithm in Gromacs 4.0.5.63 This was followed by a 20 ps MD simulation of the hydrate structure, keeping the oxygen atoms frozen, but allowing the water molecules to rotate. The water molecules in the resulting configurations were oriented such that the hydrogen atoms of any given water molecule were directed toward the

neighboring oxygen atoms. The classical (nondissociative) nature of TIP4P/2005 potential and the short MD relaxation run, therefore, resulted in configurations that were compliant with the BernalFowler rules.64 The net dipole moment of the hydrate simulation cell containing 2944 water molecules was 44% lower than the corresponding quantity for a simulation cell containing the same number of water molecules in the liquid state, at 298 K (due to finite size effects the latter is nonzero: the average over 20 snapshots is 4% of the maximum possible value). This provided the initial hydrate configuration. The cages were determined using the algorithm developed by Jacobson et al.65 This determination took the periodic boundary conditions into consideration. Carbon dioxide molecules were added to the center of the cages. The configuration of water molecules in contact with the hydrate phase was obtained from a 1 ns equilibration run of pure water at 298 K and 1 bar. The box of water was combined with the hydrate structure to obtain the starting configuration of our simulations, and was energy-minimized using the Steepest Descent algorithm in Gromacs 4.0.563 and relaxed by a 40 ps MD simulation at 273 K. Note that the hydrate and water boxes were combined along the z-dimension such that the hydrate water interfacial area was in the xy plane. The total number of water molecules in the system was 6585 and the number of CO2 molecules depended on the occupancy, as described below. A snapshot of the hydratewater system during a simulation is shown in Figure 1a. 6103

dx.doi.org/10.1021/jp110868t |J. Phys. Chem. A 2011, 115, 6102–6111


ARTICLE

Table 1. Hydrate Systems Studied (θL: Occupancy of Large Cages; θS: Occupancy of Small Cages) occupancy (%)

θL

θS

case

notation

1

100%

100

1

1

2 3

87%S 87%L

87.5 87.2

1 0.83

0.5 1

CO2 forms sI hydrate, comprising six 51262 (large) cages and two 512 (small) cages per unit cell. We use the common notation XnYm to denote a cage consisting of n X-sided and m Y-sided polygons. Occupancy is calculated as the number of cages filled per unit cell (i.e., per eight cages for sI structure). We studied three different occupancy states of the CO2 hydrate (Table 1). The two cases with overall occupancy of ca. 87% differed in the distribution of their empty cages. In one case (87%S), half of the small cages were empty, and in the other case (87%L), 17% of the large cages were empty. In the latter case, the empty large cages were selected randomly. Both the 87% overall occupancy cases had approximately one empty cage per unit cell in the hydrate structure. These two cases allowed us to explore the effect of the distribution of guest molecules in the different cage types on the kinetics and mechanisms of hydrate dissociation, while maintaining the same overall CO2 occupancy. B. Details of the Molecular Dynamics Simulations. The TIP4P/2005 model66 was used to represent water. TIP4P/2005 has been shown to reproduce a wide range of water properties with comparable accuracy to that of other commonly used classical force fields, while at the same time capturing its complex solid-state phase behavior with distinctively superior accuracy.67 Carbon dioxide was described as a rigid linear molecule with five interacting sites as per the TraPPE force field.68 Recent studies of pure CO2 in different phases by Anderson et al. suggest that CO2 flexibility plays only a limited role on its properties.69 All molecular dynamics (MD) simulations were run in the isothermalisobaric (NpT) ensemble using Gromacs 4.0.5.63 Constant temperature and pressure were maintained using a Nose-Hoover thermostat70 and a Parrinello-Rahman barostat,71 respectively. Semi-isotropic pressure coupling was used such that the box length along the z-dimension was scaled independently of the x- and y-dimensions. Simulations were run at 30.5 bar and different temperatures, from 240 to 290 K. Electrostatic interactions were calculated using the particle mesh Ewald (PME) method.72 LorentzBerthelot mixing rules were used for Lennard-Jones cross interactions between different atom types. Nonbonded interactions were cutoff at 1 nm and long-range dispersion corrections were applied to energy and pressure. SHAKE73 was used to maintain the geometry of water molecules with a relative tolerance of 104 nm. Each simulation was run for 50 ns or until the hydrate melted completely, with a time step of 2 fs. Coordinates were saved every picosecond for further analysis. C. Structural Analysis of Water Molecules. We followed the method described in previous studies47,48 to classify water molecules as solid- or liquid-like. An order parameter, q, was calculated for each water molecule. q measures the deviation in the arrangement of a given water molecule and its nearest neighbors from tetrahedrality and is defined as qi ¼

ni 1

ni

∑ ∑ ðjcos θjik jcos θjik þ 1=9Þ2 j¼1 k¼j þ 1

ð1Þ

where ni denotes the nearest neighbors of water molecule i and θjik is the angle between the oxygen atoms of water molecules j, i, and k.

We considered water molecules within 0.35 nm of the central molecule to be its nearest neighbors,47 this distance corresponding approximately to the first minimum in the water oxygenoxygen pair correlation function in the hydrate phase at 240 K and 30.5 bar. In a perfect tetrahedral arrangement, the term in parentheses of eq 1 will be zero resulting in zero q-value. Baez and Clancy47 showed that a cutoff value of 0.4 can be used to distinguish between solid- (qi < 0.4) and liquid-like (qi > 0.4) water molecules. To further distinguish between ice-like and hydrate-like water molecules, the presence of cyclic pentagons was investigated.47 Cyclic pentagons were defined as consisting five water molecules such that the edges were shorter than the nearest-neighbor distance. Water molecules in the hydrate structure are part of 4,5 or 6 cyclic pentamers. Therefore, water molecules with qi < 0.4 and part of 4, 5, or 6 cyclic pentamers were classified as hydrate-like water molecules. Ice-like water molecules were those with qi < 0.4 but not a part of any cyclic pentagon. Other water molecules were classified as liquid-like. This was followed by updating the state of the water molecule based on that of its nearest neighbors.47,48 That is, if a water molecule had three or more nearest-neighbors in a different state, then the state of that water molecule was updated to that of its nearest neighbors. In Figure 1b, water molecules are colored based on their classification as hydrate-, ice-, or liquid-like. The corresponding configuration is shown in Figure 1a, and it can be seen that the hydrate-like water molecules are correctly identified. The CO2 molecules were classified as being in the hydrate (small, large, or partial cage) or the bulk liquid phase based on the number of hydrate-like water molecules in the first solvation shell of the CO2 molecule, ns: a CO2 was considered to be in a large hydrate cage if ns > 23, in a small cage if 18 < ns < 21, and in a partial cage if 10 < ns < 18.47,48 Partial cages comprise incomplete small and large cages. Any CO2 molecule with less than 10 hydrate-like water molecules in the first hydration shell was considered to be in liquid water.47,48 A water molecule was considered to be in the hydration shell of a CO2 molecule if the distance between the CO2 carbon and the water oxygen was less than 0.6 nm. This distance corresponds to the location of the first minimum in the CO2 carbonwater oxygen radial distribution function in the hydrate.

III. RESULTS AND DISCUSSION A. Determining the Melting Point of CO2 Hydrates. The melting point is sensitive to intermolecular interactions, and therefore, we begin by determining the melting point of CO2 hydrate for the chosen force-fields. For the purpose of this work, the relevant region of the hydrate phase diagram comprises twophase regions separated by lines representing coexistence of three phases. For example, CO2 hydrate (H) is in equilibrium with liquid water (LW) in the temperature range of 273.1283.1 K and pressures of 12.546.5 bar.2 According to the phase rule, coexistence of the two phases (LW, H) in a three-component system (CO2, water, and hydrate) entails three degrees of freedom (e.g., T, P, CO2/H2O ratio). For a given pressure and composition, if the temperature is increased starting from inside this LWH region, the system approaches the three-phase (hydrateliquid watervapor) coexistence line. Beyond this three-phase coexistence temperature, the hydrate melts and the system is in the liquid watervapor (LWV) region. Therefore, one can determine the dissociation/melting temperature of the hydrate by starting in the LWH region and heating the hydrate gradually. The temperature at which the hydrate begins to melt is 6104


The Journal of Physical Chemistry A its dissociation or melting temperature. We followed this strategy to determine the melting point of CO2 hydrate in our study. Hydrates are nonstoichiometric in that not every cage is occupied by a guest molecule. Consequently, different experimental studies have reported different hydration numbers and occupancy of hydrates. This variable-composition region defines a narrow parabolic domain in the T-x plane, where the solid hydrate region, instead of being represented by a vertical line (fixed x), spans a narrow area encompassing the experimentally observed range of possible hydrate compositions.2,60,74 However, for most hydrates the composition range of the solid hydrate phase is not known accurately. In the case of CO2 hydrates, most experimental estimates of the hydration number are between 5.75 and ∼7.20,7577 which correspond to an occupancy range of 100 to 80%. We, therefore, investigated three different CO2 hydrate compositions: 100% occupancy, and two cases of 87% hydrate occupancy which differed in the distribution of CO2 between small and large cages. Further details are given in the Methods section and Table 1. The preponderance of experimental evidence suggests that CO2 occupies more than 95% of the large cages, and vacancies are mainly observed in small cages for the sI hydrate structure.2,75,76 It is, however, important to note that estimating the cage-specific occupancy from experimental results, which measure relative, not absolute, occupancies, involves approximations and assumptions.75 In addition, hydrate composition exhibits appreciable sensitivity to experimental factors such as conditions of sample preparation75 (e.g., growth near a vaporliquid interface or in bulk74). Understanding and rationally manipulating such experimental variables could allow the possibility of producing hydrate phases with tunable and controllable composition. Recent results from molecular dynamics simulations of hydrate growth identified empty large cages (51262 and 51264) as defects in the growing hydrate lattice.40,43 The surrounding occupied cages stabilized these unoccupied large cages. These observations indicate that it is possible to form hydrate structures with non-negligible fraction of empty large cages. Motivated by this, we investigate two cases of 87% overall occupancy in this study, including one where 17% of the large cages are empty (see Table 1), allowing us to explore the spectrum of possible dissociation behaviors. To determine the melting point of CO2 hydrate, we simulated the hydratewater system at different temperatures and observed the change in the water density profile within the hydrate region as a function of time. The melting point was taken to be the lowest temperature at which changes in the water density profile within the hydrate region were observed during the 50 ns of the simulation. Initially, the hydrate phase spanned the 0.6 e z e 4.0 nm region of the computational cell, and excluded the partial cages. This enabled us to identify dissociation due to perturbation in complete cages, and not due to diffusion of water molecules from partial cages present in the starting configuration of the simulations. In contrast to previous studies on ice,61 we found that for our hydratewater systems the melting point could not be determined by exclusively following the energy of the system. Diffusion of CO2 and water molecules from the partial cages at the hydratewater interface into the liquid phase led to an increase in the energy of the hydratewater system at temperatures where negligible changes in hydrate structure were observed. Figure 2 shows the density profile of water in the hydrate structure for 87%S and 87%L (see Methods section and Table 1 for notation). The region enclosed within the two dashed lines

ARTICLE

Figure 2. Number density of water molecules in the hydratewater simulation system along the z-dimension. The density profile after 50 ns simulation at 240 K is shown in red in all panels. The density profiles in blue are those obtained after 50 ns simulation at (a) 254 K for 87%S, (b) 256 K for 87%S, (c) 252 K for 87%L, and (d) 254 K for 87%L. The region within the dashed lines represents the initial hydrate phase, within which there are no partial cages. The density profiles are averaged over 100 ps.

represents the initial hydrate region, just outside of which are the partial cages and hydrate-water interface. For 87%S, although the peak heights near the interfacial region decrease at 254 K (Figure 2a), the inside layers are only perturbed at 256 K (Figure 2b). In the case of 87%L, the density profile within the hydrate region begins to change at 254 K (Figure 2c,d). At the higher temperatures (256 K for 87%S; 254 for 87%L), a clear decrease in peak heights within the hydrate region, and disappearance of peaks close to the interfacial regions (but inside the hydrate region) was observed. Based on this measure, the estimates of melting points are approximately 258, 256, and 254 K for 100%, 87%S, and 87%L, respectively. Since we used increments of 2 K in our simulations, the error on the dissociation temperature of the hydrates is also 2 K. Thus, there is only a modest dependence of dissociation temperature on hydrate occupancy. This is consistent with the narrow parabolic region representing the hydrate phase in the T-x phase diagram,2,60,74 which suggests that hydrates with different composition can have similar melting points. The experimentally determined melting point of CO2 hydrate at 30.5 bar is 280.3 K.78 While full quantitative agreement between the experimental and simulation-based melting point cannot be expected due to the limitations of the water and CO2 force-fields, qualitative trends can be compared. Specifically, the melting point of CO2 hydrate at 30.5 bar estimated in our simulations is higher than the melting point of TIP4P/2005 ice (=250 K61) for all cases of hydrate occupancy studied here, consistent with the experimental results.78 The difference in the melting points (=Tm(hydrate)Tm(ice)) of 48 K obtained from simulations is in good agreement with the experimental value of 7.15 K. B. Rate and Mechanism of Hydrate Dissociation. 1. Rate of Hydrate Dissociation. When the hydrate melts during the simulation, the thermostat provides energy to the system, and a steady 6105



ARTICLE

Table 2. Time Taken (in ns) by the Hydrate Slab to Melt: The Time of Melting Is Taken to Be the Approximate Time at Which the Energy of the HydrateWater System Plateaus

Figure 3. Evolution of the energy in the hydrate-water system, E relative to the final bulk system energy, ÆEeæ at 270, 280, and 290 K for the three occupancy cases studied. ÆEeæ was averaged over 2425, 78, and 34 ns for 270, 280, and 290 K, respectively. The increase in energy indicates dissociation of the CO2 hydrate. Because both E and ÆEeæ are negative, an increase in ÆEeæ/E corresponds to an increase in E (i.e., E becomes less negative).

increase in energy is observed. As seen in Figure 3, at high temperatures the energy increases consistently to reach a final plateau indicating that the hydrate slab melts completely. The steep increase in the energy just prior to the final plateau corresponds to the melting of the last layer of the hydrate. This steep increase has also been observed in the melting of ice.61 The rate at which the energy increases is correlated to the rate of melting. The 87%L CO2 hydrate dissociates faster than its 87%S and 100% counterparts. Table 2 shows the time taken for the entire hydrate slab in the simulation to melt at high temperatures. For comparison, the time taken for the empty hydrate to melt is also given. Clearly, the sI hydrate structure is unstable in the absence of CO2 molecules and collapses within few picoseconds of the simulation. The small difference between the times taken by the 100% and 87%S hydrates to melt suggests that the rate of melting does not change significantly when the occupancy of the hydrate is decreased from 100 to 87% by removing CO2 molecules from the small cages. On the other hand, the considerably smaller time taken for 87%L hydrate to melt compared to the 100% case indicates that removing guest molecules from large cages increases the rate of dissociation. Thus, we see that the kinetics of hydrate dissociation is sensitive to cage-specific occupancy, not just to overall occupancy. This is an important observation, because previous simulation and modeling studies48,49,51,52 of hydrate dissociation have focused only on the overall occupancy. This may also underlie the discrepancy between the results of Mysakhin et al.48 and English et al.49 on methane hydrate dissociation. Mysakhin et al.48 observed the rate of hydrate dissociation changing with hydrate composition for a slab geometry,

temperature (K)

100%

87%S

87%L

empty

270

19.0

19.3

9.2

0.26

280

5.1

5.2

2.7

0.16

290

1.9

1.6

0.9

Figure 4. Temperature dependence of the rate of CO2 release (rCO2). (a) The rate of CO2 release calculated from linear fits to the number of CO2 molecules released into bulk water as a function of time. (b) Arrhenius plot. Error bars were calculated from five simulation runs at the given temperature.

while English et al. did not see any significant composition dependence on the dissociation rate of spherical nanocrystals of methane hydrate.49 English et al.49 do not report the cagespecific occupancy. We further quantified the rate of dissociation in the following sections by calculating the rate of CO2 release into bulk water, and the hydratewater interface velocity. 2. Rate of CO2 Release from the Hydrate. The rate of CO2 molecules released into bulk water is correlated with the rate of hydrate dissociation. A CO2 molecule is considered to be released into bulk water if its classification changes from being in the hydrate to being in bulk liquid phase during the simulation (classification criteria are defined in the Methods section). The number of CO2 molecules released into bulk water in the initial part of the simulation (e.g., 025 ns for T < 260 K) varied approximately linearly with time and the slope of the line gave the initial rate of CO2 release (rCO2), shown in Figure 4a. Clearly, rCO2 is similar in the case of 100% and 87%S, while that for 87%L is higher. Consistent with results from the previous section, this indicates that, for the same overall occupancy, the presence of empty small or large cages has different effects on the kinetics of hydrate dissociation. The temperature dependence of the rate of CO2 release displays Arrhenius behavior (Figure 4b). The activation energies are 62.30, 66.72, and 66.17 kJ/mol for 100%, 87%L, and 87%S, respectively. Interestingly, the magnitude of the activation energy 6106



ARTICLE

Figure 5. (a) Arrhenius plot for rate of CO2 release, (rCO2) for 87%S hydrate. (b) Fraction of solid water molecules along the z-dimension of the simulation box for 87%S hydrate at 240 K (average over 49.549.6 ns time interval). Black points and curve: Data when water molecules were classified as solid- or liquid-like considering both q-values and cyclic pentamers. Magenta points and curve: Data when classification of water molecules was based on q-value alone.

in all three cases is similar within the accuracy of the fitting, indicating that the faster rate of dissociation of 87%L originates primarily in the pre-exponential factor. The value of the activation energies reported above corresponds approximately to the breaking of three hydrogen bonds. We also calculated the rate of CO2 released into bulk in the case where water molecules were classified as solid- or liquid-like based only on their q-value. No distinction between ice- and hydrate-like molecules was made in this approach. While the rates obtained by either methods are similar at high temperatures, the latter method underestimates the rate of CO2 release at low temperatures, see Figure 5a. The classification of water molecules is sensitive to the method used primarily in the interfacial region, Figure 5b. The approach that does not consider pentagons (i.e., based only on the local tetrahedrality) results in higher fraction of solid-like water molecules in this region and, therefore, fewer CO2 molecules in bulk water. This leads to lower estimates of CO2 release rate. The considerable difference in the CO2 release rates obtained from the two methods at low temperature indicates that at these temperatures released CO2 molecules are concentrated near the interfacial area and do not diffuse away into bulk water. Previous methane hydrate dissociation studies have also reported that at lower temperatures hydrate dissociation is characterized by diffusion of methane molecules out of partial cages at the hydrate-water interface.48,51 3. Velocity of the HydrateWater Interface. Figure 5b shows the fraction of solid-like water molecules (fsolid) along the length of the simulation box. We found that almost no ice-like water molecules were identified in most cases. Therefore, solid-like water molecules mostly comprise hydrate-like water molecules and fsolid ∼ fhydrate. In the hydrate region, fsolid ∼1 and in the liquid water phase, fsolid ∼ 0. There is a smooth transition from the hydrate to liquid water and the region over which this transition occurs is the hydratewater interface. The width of the hydratewater interface was measured by fitting the curves to the functional form:

fsolid

z z0 ¼ 0:5ðfH þ fL Þ þ 0:5ðfH fL Þ tanh d

ð2Þ

where fH, fL, z0, and d are fitting parameters that represent the fraction of solid-like water molecules in the hydrate and bulk liquid phase, the z-value at which fsolid = 0.5 (fH þ fL), and half of interface width, respectively.

Figure 6. Fraction of solid water molecules along the z-dimension of the simulation box. The fraction of solid water molecules are calculated in slabs of finite thickness and averaged over 100 ps. z-Dimension represents the direction normal to the hydratewater interface. Fraction of solid water molecules (a) at 240 K for 87%S after 0.1 (black), 25.1 (red), and 49.9 (blue) ns of simulation; and (b) at 270 K for 87%S after 0.2 (black), 2.2 (red), 5.2 (blue), and 14.2 (gray) ns of simulation. Arrows point in the direction of increasing time.

Figure 5b compares the fsolid profile computed using only the q-value to distinguish between solid- and liquid-like water molecules to the corresponding profile obtained using both qvalue and participation in cyclic pentamers to distinguish between solid- and liquid-like water molecules. The latter method gave rise to a narrower interface. This indicates that the water molecules near the hydrate display higher tetrahedrality than in bulk water. The width of the hydratewater interface obtained from the fitting varied between 0.3 and 0.6 nm. We did not observe any systematic increase or decrease in the interfacial width (not to be confused with the thickness of the entire hydrate slab) as the hydrate dissociated. Further, no clear temperature dependence of the interface width was observed. The change in the thickness of the hydrate (high fsolid) region with time provides a measure of the rate of hydrate dissociation. We measured the thickness of the hydrate region as the distance between the location at which fsolid = 0.8 at the two hydrate water interfaces. The rate of change of this quantity gave the hydratewater interface velocity. At temperatures below 250 K, negligible changes in the hydrate thickness were observed. At higher temperatures, a clear decrease in the hydrate thickness was observed during the simulation (Figure 6). Therefore, estimates of the hydratewater interface velocity were obtained only for T g 270 K, where a linear trend was observed, and are shown in Figure 7. At temperatures between 250 and 270 K, although the hydrate thickness decreased during the simulation, the decrease was not linear in time. In some cases we observed that intervals of decrease in hydrate thickness were separated by periods of negligible change in hydrate thickness. This suggests a different mechanism of hydrate dissociation at lower temperatures in which the rate-limiting step is probably the nucleation of a surface defect. We are investigating this in ongoing work. Figure 7 shows the temperature dependence of the hydratewater interface velocity. The velocity increases from ∼0.1 m/s at 270 K to ∼1.0 m/s at 290 K for 87%S and 100% hydrate. Similar velocities have been calculated for the icewater interface by Nicolson et al.79 Consistent with the discussion in the previous section, the hydratewater interface velocity results also indicate that empty large cages decrease the stability of sI structure of CO2 hydrate more than empty small cages. The interface velocity displays Arrhenius temperature dependence, as seen in Figure 7. The activation energies obtained from the fits are 68.4, 69.8, and 69.3 kJ/mol for 87%S, 87%L, and 100% 6107



ARTICLE

Figure 9. Dissociation of empty and filled large (51262) cages for 87%L occupancy hydrate, at 260 and 270 K. Calculations correspond to cages located initially along an xy plane such that the z-coordinate of their center of mass was approximately at z = 3.05 nm (see Figure 6 for location of coordinate origin).

Figure 7. Temperature dependence of the hydratewater interface velocity. (a) Hydratewater interface velocity at different temperatures (T). (b) Arrhenius plot. Error bars were calculated from five simulation runs at the given temperature.

Figure 8. Dissociation of empty and filled small (512) cages for 87%S occupancy hydrate, at 260 and 270 K. Different colors represent different cages along the same initial z-coordinate. Calculations correspond to cages located initially along an xy plane such that the z-coordinate of their center of mass was approximately at z = 0.7 nm for the empty small cages and at z = 1.3 nm for the filled small cages (see Figure 6 for location of coordinate origin).

occupancy, respectively. Interestingly, the activation energies for the three cases are similar within numerical accuracy. This result is consistent with the activation energies calculated from the temperature dependence of the rate of CO2 release into bulk water. Both estimates suggest that the activation energy does not display any significant occupancy dependence, and corresponds to the breaking of ∼35 hydrogen bonds, confirming that the faster dissociation kinetics of 87%L are reflected entirely in the pre-exponential factor.

C. Cage Dissociation. We explored the differences in the dissociation of individual cages in order to investigate the higher rate of dissociation of 87%L compared to that of 87%S. The dissociation of a cage was followed by monitoring the number of intact edges in that cage. Edges of a given cage were considered to be intact as long as the distance between the involved water molecules was within 0.35 nm. This distance is the same as the cutoff used to determine the nearest-neighbors of a water molecule in the calculation of the order parameter q,47 and in the identification of cages. As the cage begins to melt/dissociate the number of intact edges in the cage decreases. Representative trajectories of cage dissociation are shown in Figures 810. The trajectories are shown for cages in the same z-plane in the initial configuration of the simulation to enable comparison between cages with similar dissociation onset times. Consistent with a previous study,80 we found that filled 512 cages in 87%S persist longer than empty cages (Figure 8). It can be seen that the longer lifetime of filled 512 cages compared to empty ones persists even at higher temperatures. The time taken for filled cages to dissociate completely is approximately twice that of the empty cages at both 260 and 270 K. Note the similarities in the trajectories of a given cage type in the same z-plane (e.g., trajectories of empty small cages at 260 K shown in Figure 8), suggesting relatively little variation in dissociation behavior across the hydratewater interface for a given type of cage. This has also been observed in previous studies of methane hydrate decomposition.47,4952 This distinct stabilization of 512 cages by guest molecules probably also contributes to their considerably earlier appearance compared to larger cages in homogeneous nucleation of hydrates as observed by Rodger and co-workers36,42 and Walsh et al.34 In interesting contrast, the dissociation of large (51262) cages in 87%L does not show appreciable differences between filled and empty cages (Figure 9) within the statistics sampled in our simulations. Both filled and empty large cages dissociate completely in ∼3040 ns at 260 K and ∼48 ns at 270 K. It thus appears that in the case of 51262 cages, the presence of CO2 does not affect the lifetime of the cage. In Figure 10, trajectories of filled small and large cages are shown for 87%L and 87%S. It can be seen that the trajectories 6108



Figure 10. Dissociation of filled small (512) and large (51262) cages for 87%L and 87%S occupancy hydrates at 270 K. Calculations correspond to cages located initially along an xy plane such that the z-coordinate of their center of mass was approximately at z = 1.85 nm (see Figure 6 for location of coordinate origin).

corresponding to the same cage differ between the two hydrates having the same overall fractional occupancy. In the 87%L hydrate, both 512 and 51262 cages dissociate twice as fast as in case of 87%S. Thus, the time taken for similar filled cages to dissociate is different in the two hydrates with the same overall occupancy but different cage-specific occupancy. Therefore, individual cage behavior is dependent on the surrounding hydrate structure. In a different simulation study, it was reported that the hydrate structure is influenced by presence of empty cages.81 It was found that the lattice constant of xenon hydrate crystal with an empty small cage was lower than that the corresponding fully-occupied crystal. However, negligible changes were observed in the presence of an empty large cage.81 Our results along with these observations suggest that hydrate dissociation is a collective phenomenon that cannot be described based on the simple distinction between small and large or filled and empty cages. The environment surrounding a cage plays a significant role in determining its dissociation rate.

IV. CONCLUSIONS We used molecular dynamics simulations to study the dissociation of CO2 hydrates in contact with liquid water. We explored the effect of cage occupancy on the rate of hydrate dissociation, and investigated the effect of the presence of empty small and empty large cages while maintaining the same overall hydrate occupancy. We quantified the dissociation process through the rate of CO2 release from the hydrate into bulk water and through the velocity of the hydratewater interface. Our study reveals two novel aspects of hydrate dissociation: In the first place, the rate of hydrate dissociation is dependent not only on the overall occupancy, but also on the cage-specific occupancy. Therefore, it is important to consider the relative occupancies of the small and large cages, not just overall hydrate occupancy, when considering stability and kinetics of hydrate dissociation. Second, hydrate dissociation is cooperative and cannot be explained based on the simple distinction between small and large or filled and empty cages. Furthermore, the dissociation of cages is dependent on their surrounding hydrate structure. It will indeed be interesting to explore this behavior in experiments. Experimental studies74,75 indicate that hydrate

ARTICLE

composition is sensitive to factors such as conditions of sample preparation (e.g., growth near vaporliquid interface or in bulk). Recent simulation studies43 show evidence that faster rate of hydrate growth results in lower hydrate occupancy. These observations suggest that it may be possible to tune hydrate occupancy, enabling systematic studies focused on the dependence of dissociation behavior on hydrate composition. The present study suggests several additional directions for future investigations. Understanding the microscopic basis of the difference in pre-exponential factors governing the dissociation rate, and the insensitivity of the activation energy to details of cage occupancy is important. We plan to investigate these questions in future work. To further develop our understanding of hydrate dissociation, similar studies for sII and sH hydrate structures should be performed. Computational evidence suggests that different crystallographic faces display different melting behavior in case of ice.82 Whether the same is true for hydrate structures should be investigated. It has been observed that hydrates of large molecules can be stabilized by using small helper gas molecules (e.g., isobutane þ methane2). It would be interesting to explore the role of cooperativity in the formation and dissociation of mixed hydrates. Additionally, effects of confinement and of the presence of salt on the kinetics of hydrate dissociation, and on the cooperative behavior reported herein merit investigation. Both these conditions are relevant to storage of CO2 as hydrates in sediments below the ocean floor. Our work so far has focused on hydrate dissociation. The formation kinetics is of course of paramount importance.83 Its investigations using state-of-the-art path sampling techniques,8488 including also studies of cooperative behavior in hydrate formation,34 is being pursued in ongoing work. Ultimately, the goal is to develop a general framework for the computational study of hydrate structure and stability, including both formation and dissociation.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected].

’ ACKNOWLEDGMENT We thank Prof. Valeria Molinero and Liam C. Jacobson for kindly sharing their code for determining cages in the hydrate structure with us. This work was supported by funding from BP Amoco as part of Carbon Mitigation Initiative at Princeton University. Computations were performed at the Terascale Infrastructure for Groundbreaking Research in Engineering and Science facility at Princeton University. ’ REFERENCES (1) Makogon, I. F.; Makogon, Y. F. Hydrates of Hydrocarbons, 1st ed.; Pennwell Books: Tulsa, OK, 1997. (2) Sloan, E. D.; Koh, C. A. Clathrate Hydrates of Natural Gases, 3rd ed.; CRC Press/Taylor-Francis: Boca Raton, FL, 2008. (3) Boswell, R.; Collett, T. S. Energy Environ. Sci. 2011, Advance article; doi: 10.1039/C0EE00203H. (4) Makogon, Y.; Holditch, S.; Makogon, T. J. Pet. Sci. Eng. 2007, 56, 14–31. (5) Milkov, A. V. Earth-Sci. Rev. 2004, 66, 183–197. (6) Lang, X.; Fan, S.; Wang, Y. J. Nat. Gas Chem. 2010, 19, 203–209. (7) Struzhkin, V. V.; Militzer, B.; Mao, W. L.; Mao, H.-K.; Hemley, R. J. Chem. Rev. 2007, 107, 4133–4151. 6109


The Journal of Physical Chemistry A (8) Lee, H.; Lee, J.; Kim, D.; Park, J.; Seo, Y.; Zeng, H.; Moudrakovski, I.; Ratcliffe, C.; Ripmeester, J. Nature 2005, 434, 743–746. (9) Javanmardi, J.; Moshfeghian, M. Appl. Ther. Eng. 2003, 23, 845–857. (10) J., B. A. Chem. Eng. Prog. 1975, 71, 80–87. (11) Knox, W. G.; M., H.; Jones, G. E.; Smith, H. B. Chem. Eng. Prog. 1961, 57, 66–71. (12) Hao, W.; Wang, J.; Fan, S.; Hao, W. Energy Convers. Manage. 2008, 49, 2546–2553. (13) Bortnowska, M. Pol. Maritime Res. 2009, 16, 70–78. (14) 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–12694. (15) Javanmardi, J.; Nasrifar, K.; Najibi, S.; Moshfeghian, M. Appl. Ther. Eng. 2005, 25, 1708–1723. (16) Mao, W. L.; Mao, H. K. Proc. Natl. Acad. Sci. U.S.A. 2004, 101, 708–710. (17) Cha, I.; Lee, S.; Lee, J. D.; Lee, G. W.; Seo, Y. Environ. Sci. Technol. 2010, 44, 6117–6122. (18) Peng, X. M.; Hu, Y. F.; Liu, Y. S.; Jin, C. W.; Lin, H. J. J. Nat. Gas Chem. 2010, 19, 81–85. (19) Seo, Y.; Moudrakovski, I. L.; Ripmeester, J. A.; Lee, J.; Lee, H. Environ. Sci. Technol. 2005, 39, 2315–2319. (20) Sloan, E. D. Nature 2003, 426, 353–363. (21) Rochelle, C. A.; Camps, A. P.; Long, D.; Milodowski, A.; Bateman, K.; Gunn, D.; Jackson, P.; Lovell, M. A.; Rees, J. Geol. Soc. Spec. Publ. 2009, 319, 171–183. (22) Tohidi, B.; Yang, J.; Salehabadi, M.; Anderson, R.; Chapoy, A. Environ. Sci. Technol. 2010, 44, 1509–1514. (23) Tester, J. W.; Bivins, R. L.; Herrick, C. C. AIChE J. 1972, 18, 1220. (24) Tse, J. S.; Klein, M. L.; McDonald, I. R. J. Chem. Phys. 1983, 78, 2096–2097. (25) Tse, J. S.; Klein, M. L.; McDonald, I. R. J. Phys. Chem. 1983, 87, 4198–4203. (26) Geng, C. Y.; Wen, H.; Zhou, H. J. Phys. Chem. A 2010, 113, 5463–5469. (27) Yezdimer, E. M.; Cummings, P. T.; Chialvo, A. A. J. Phys. Chem. A 2002, 106, 7982–7987. (28) Papadimitriou, N. I.; Tsimpanogiannis, I. N.; Stubos, A. K.; Martin, A.; Rovetto, L. J.; Peters, C. J. J. Phys. Chem. Lett. 2010, 1, 1014–1017. (29) Susilo, R.; Alavi, S.; Ripmeester, J.; Englezos, P. Fluid Phase Equilib. 2008, 263, 6–17. (30) Susilo, R.; Alavi, S.; Ripmeester, J. A.; Englezos, P. J. Chem. Phys. 2008, 128, 194505. (31) Alavi, S.; Ripmeester, J.; Klug, D. J. Chem. Phys. 2006, 124, 204707. (32) Jacobson, L. C.; Hujo, W.; Molinero, V. J. Am. Chem. Soc. 2010, 132, 11806–11811. (33) Jacobson, L. C.; Hujo, W.; Molinero, V. J. Phys. Chem. B 2010, 114, 13796–13807. (34) Walsh, M. R.; Koh, C. A.; Sloan, E. D.; Sum, A.; Wu, D. T. Science 2009, 326, 1095–1098. (35) Radhakrishnan, R.; Trout, B. L. J. Chem. Phys. 2002, 117, 1786–1796. (36) Moon, C.; Taylor, P. C.; Rodger, P. M. J. Am. Chem. Soc. 2003, 125, 4706–4707. (37) Moon, C.; Hawtin, R. W.; Rodger, P. M. Faraday Discuss. 2007, 136, 367–382. (38) Tung, Y.-T.; Chen, L.-J.; Chen, Y.-P.; Lin, S.-T. J. Phys. Chem. B 2010, 114, 10804–10813. (39) Nada, H. J. Phys. Chem. B 2010, 113, 4790–4798. (40) Liang, S. A.; Kusalik, P. G. J. Phys. Chem. B 2010, 114, 9563–9571. (41) Liang, S. A.; Kusalik, P. G. Chem. Phys. Lett. 2010, 494, 123–133. (42) Hawtin, R. W.; Quigley, D.; Rodger, P. M. Phys. Chem. Chem. Phys. 2008, 10, 4853–4864.

ARTICLE

(43) Vatamanu, J.; Kusalik, P. G. J. Phys. Chem. B 2006, 110, 15896–15904. (44) Park, S.; Sposito, G. J. Phys. Chem. B 2003, 107, 2281–2290. (45) Iwai, Y.; Nakamura, H.; Arai, Y.; Shimoyama, Y. Mol. Simul. 2010, 36, 246–253. (46) Rodger, P. Ann. N.Y. Acad. Sci. 2000, 912, 474–482. (47) Baez, L.; Clancy, P. Ann. N.Y. Acad. Sci. 1994, 715, 177–186. (48) Myshakin, E. M.; Jiang, H.; Warzinski, R. P.; Jordan, K. D. J. Phys. Chem. A 2009, 113, 1913–1921. (49) English, N. J.; Johnson, J. K.; Taylor, C. E. J. Chem. Phys. 2005, 123, 244503. (50) Wan, L. H.; Yan, K. F.; Li, X. S.; Huang, N. S.; Tang, L. G. Acta Chim. Sin. 2009, 67, 2149–2154. (51) English, N. J.; Phelan, G. M. J. Chem. Phys. 2009, 131, 074704. (52) Alavi, S.; Ripmeester, J. A. J. Chem. Phys. 2010, 132, 144703. (53) Wathen, B. Kwan, P. Jia, Z. Walker, V. High Performance Computing Systems and Applications; Springer: Berlin/Heidelberg, 2010; Vol. 5976; pp 117133. (54) Kvamme, B.; Kuznetsova, T.; Aasoldsen, K. Mol. Simul. 2005, 31, 1083–1094. (55) Hawtin, R. W.; Rodger, P. M. J. Mater. Chem. 2006, 16, 1934–1942. (56) Moon, C.; Taylor, P. C.; Rodger, P. M. Can. J. Phys. 2003, 81, 451–457. (57) Kvamme, B.; Huseby, G.; Forrisdahl, O. K. Mol. Phys. 1997, 90, 979–991. (58) Wallqvist, A. J. Chem. Phys. 1992, 96, 5377–5382. (59) Chialvo, A.; Houssa, M.; Cummings, P. J. Phys. Chem. B 2002, 106, 442–451. (60) Wierzcgowski, S.; Monson, P. Ind. Eng. Chem. Res. 2006, 45, 424–431. (61) Fernandez, R.; Abascal, J.; Vega, C. J. Chem. Phys. 2006, 124, 144506. (62) Hollander, F.; Jeffrey, G. A. J. Chem. Phys. 1977, 66, 4699–4705. (63) Hess, B.; Kutzner, C.; Van Der Spoel, D.; Lindahl, E. J. Chem. Theory Comput. 2008, 4, 435. (64) Bernal, J. D.; Fowler, R. H. J. Chem. Phys. 1933, 1, 515–548. (65) Jacobson, L. C.; Hujo, W.; Molinero, V. J. Phys. Chem. B 2009, 113, 10298–10307. (66) Abascal, J. L. F.; Vega, C. J. Chem. Phys. 2005, 123, 234505. (67) Vega, C.; Abascal, J. L. F.; Conde, M. M.; Aragones, J. L. Faraday Discuss. 2009, 141, 251–276. (68) Potoff, J. J.; Siepmann, I. J. AIChE J. 2001, 47, 1676–1682. (69) Anderson, K. E.; Mielke, S. L.; Siepmann, J. I.; Truhlar, D. G. J. Phys. Chem. A 2009, 113, 2053–2059. (70) Nose, S. Mol. Phys. 1984, 52, 255–268. (71) Parrinello, M.; Rahman, A. J. Appl. Phys. 1981, 52, 7182–7190. (72) Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. 1993, 98, 10089–10092. (73) Ryckaert, J.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327–341. (74) Huo, Z.; Hester, K.; Sloan, E. D., Jr,; Miller, K. T. AIChE J. 2003, 49, 1300–1306. (75) Uchida, T. Waste Manage. 1997, 17, 343. (76) Sun, R.; Duan, Z. Geochim. Cosmochim. Acta 2005, 69, 4411–4424. (77) Circone, S.; Stern, L. A.; Kirby, S. H.; Durham, W. B.; Chakoumakos, B. C.; Rawn, C. J.; Rondinone, A. J.; Ishii, Y. J. Phys. Chem. B 2003, 107, 5529–5539. (78) Wendland, M.; Hasse, H.; Maurer, G. J. Chem. Eng. Data 1999, 44, 901–906. (79) Nicholson, B. F.; Clancy, P.; Rick, S. W. J. Cryst. Growth 2006, 293, 78–85. (80) Guo, G.; Zhang, Y.; Zhao, Y.; Refson, K.; Shan, G. J. Chem. Phys. 2004, 121, 1542–1547. (81) Ikeda-Fukazawa, T.; Yamaguchi, Y.; Nagashima, K.; Kawamura, K. J. Chem. Phys. 2008, 129, 224506. 6110



ARTICLE

(82) Nada, H.; Furukawa, Y. Surf. Sci. 2000, 446, 1–16. (83) Debenedetti, P. G.; Sarupria, S. Science 2009, 326, 1070– 1071. (84) Allen, R. J.; Valeriani, C.; ten Wolde, P. R. J. Phys.: Condens. Matter 2009, 21, 463102. (85) Escobedo, F. A.; Borrero, E. E.; Araque, J. C. J. Phys.: Condens. Matter 2009, 21, 333101. (86) Dellago, C.; Bolhuis, P. G. Advanced Computer Simulation Approaches for Soft Matter Sciences III; Advances in Polymer Science; Springer-Verlag: Berlin, 2009; Vol. 221; pp 167233. (87) van Erp, T. S.; Bolhuis, P. G. J. Comput. Phys. 2005, 205, 157–181. (88) Bolhuis, P. G.; Chandler, D.; Dellago, C.; Geissler, P. L. Annu. Rev. Phys. Chem. 2002, 53, 291–318.

6111


Molecular Dynamics Study of Carbon Dioxide Hydrate Dissociation

Recommend Documents