Kinetics of CH4 and CO2 Hydrate Dissociation and Gas Bubble

Article pubs.acs.org/JPCA

Kinetics of CH4 and CO2 Hydrate Dissociation and Gas Bubble Evolution via MD Simulation M. Uddin* Alberta Innovates - Technology Futures, Edmonton, Alberta T6N 1E4, Canada

D. Coombe Computer Modeling Group Ltd., Calgary, Alberta T2L 2A6, Canada ABSTRACT: Molecular dynamics simulations of gas hydrate dissociation comparing the behavior of CH4 and CO2 hydrates are presented. These simulations were based on a structurally correct theoretical gas hydrate crystal, coexisting with water. The MD system was first initialized and stabilized via a thorough energy minimization, constant volume−temperature ensemble and constant volume− energy ensemble simulations before proceeding to constant pressure−temperature simulations for targeted dissociation pressure and temperature responses. Gas bubble evolution mechanisms are demonstrated as well as key investigative properties such as system volume, density, energy, mean square displacements of the guest molecules, radial distribution functions, H2O order parameter, and statistics of hydrogen bonds. These simulations have established the essential similarities between CH4 and CO2 hydrate dissociation. The limiting behaviors at lower temperature (no dissociation) and higher temperature (complete melting and formation of a gas bubble) have been illustrated for both hydrates. Due to the shift in the known hydrate stability curves between guest molecules caused by the choice of water model as noted by other authors, the intermediate behavior (e.g., 260 K) showed distinct differences however. Also, because of the more hydrogen-bonding capability of CO2 in water, as reflected in its molecular parameters, higher solubility of dissociated CO2 in water was observed with a consequence of a smaller size of gas bubble formation. Additionally, a novel method for analyzing hydrate dissociation based on H-bond breakage has been proposed and used to quantify the dissociation behaviors of both CH4 and CO2 hydrates. Activation energies Ea values from our MD studies were obtained and evaluated against several other published laboratory and MD values. Intrinsic rate constants were estimated and upscaled. A kinetic reaction model consistent with macroscale fitted kinetic models has been proposed to indicate the macroscopic consequences of this analysis.

1.1. Gas Hydrate Dissociation and Gas Bubble Formation. Observations from Kim et al.7 from laboratory tests indicate that the CH4 hydrate decomposition rate is first order in fugacity difference, the temperature effect follows the Arrhenius relation, and the dissociation rate depends on stirring rate. Additionally, it appeared that gas bubble accumulation at the dissociation interface and the localized temperature gradient at the dissociation surface (effect of heat transfer) can significantly reduce the dissociation rate. Uddin et al.5 refitted the test data from Kim et al.7 to explore the effects of mass transfer on dissociation. We therefore interpret the effect of stirring rate as seen by Kim et al.7 as an efficient removal of newly form gas bubbles to allow intrinsic hydrate dissociation. As indicated by Ripmeester et al.8 in their molecular dynamics

1.0. INTRODUCTION The exploitation of natural gas deposits worldwide may provide the ultimate solution to diminishing hydrocarbon reserves due to the enormous quantity of hydrocarbons stored in them. The major determinates will be to establish economically viable and environmentally friendly production methods to achieve this. Exploratory field tests in the high Arctic on CH4 hydrate dissociation1−3 and CO2 hydrate replacement (Alaska gas hydrate production trial using CO2/CH4 exchange, www.netl.doe. gov) have been recently undertaken, supported by extensive laboratory research, extensive field monitoring, and field simulations of the process. We have been involved in much of this latter work4−6 and have developed useful macroscopic kinetic representations of both processes. However, such models require fitted parameters, and a more fundamental understanding of the nature of these parameters seems advisable. Here we will focus on molecular dynamics simulation as a method to establish such parameters. © 2014 American Chemical Society

Received: November 1, 2013 Revised: February 22, 2014 Published: February 26, 2014 1971

dx.doi.org/10.1021/jp410789j | J. Phys. Chem. A 2014, 118, 1971−1988

The Journal of Physical Chemistry A

Article

calculation of methane hydrate dissociation, the dissociation rate is expected to decline over time as gas bubble formation inhibits further hydrate decay. This they attribute to a mass transfer resistance cause by the presence of newly form gas bubbles. The above independent studies suggest a significant uncertainty in any assumed constant intrinsic rate of hydrate dissociation. Follow-up studies by Clarke and Bishnoi9−11 indicate analogous phenomena for CO2 hydrate, C2H6 hydrate, and hydrate mixture dissociation as well. More concretely, the driving force for hydrate dissociation under various temperature and pressure conditions is expressed as the difference (or distance) from the hydrate stability curve. Figure 1 summarizes the known bulk scale hydrate stability Figure 2. Schematic showing the gas evolution process in hydrate dissociation: (I) hydrate phase containing highly concentrated gas molecules trapped within the hydrogen bonded cages; (II) water phase containing dispersed gas bubbles and hydrate containing dissociated gas molecules, rapidly partitioned into dissolved gas and dispersed bubble components (as the concentration is much greater than the supersaturation limit); (III) water phase containing free gas, where the gas phase is evolved by connecting or the coalescence of the dispersed bubbles.

and the possibility of hydrate crystal decomposition and reformation.

2.0. STRUCTURALLY CORRECT SI CH4 AND CO2 HYDRATES 2.1. Setup SI Hydrate Crystal. The hydrate crystal structures that were discovered in natural environment are two cubic structures, known as structure I and structure II, and one hexagonal structure, known as structure H; see Lu et al.15 This paper focuses on properties and dissociation of structure I CH4 and CO2 hydrates, only. Structure I is a cubic structure with the space group pm3n̅ and a lattice constant of about 12 Å. A unit cell of structure I hydrate contains 6 tetrakaidecahedrons and 2 pentagonal dodecahedrons cavities or cages. The 14sided cavity is a tetrakaidecahedron consisting of 12 pentagonal and 2 hexagonal faces (51262). The 12-sided cavity has twelve pentagonal faces (512) and is a pentagonal dodecahedron. The pentagonal dodecahedron or small cage is a building unit to both structure I and structure II hydrates. The cavity is almost spherical with an average cavity radius of about 3.95 Å. The small cages are situated on the corners and in the middle of the unit cell, and their centers have cubic site symmetry (m3). The small cavities are connected through their vertices, leading to somewhat larger, oblate spaces between them. These spaces are large ellipsoidal (51262) cages, which are centered on points with tetragonal 42̅m site symmetry. The average cavity radius of the large cages is 4.33 Å. It can thus be occupied by molecules with a diameter of up to 6 Å. The hydrate structure I is therefore formed by guest molecules with an average van der Waals diameter of 4.2 to 6 Å (e.g., CH4, H2S, CO2, or C2H6). Figure 3 shows the side (x or y-axis) and top (z-axis) views of the minimal hydrate crystal structure, respectively. In particular, the x and y faces have two large hexagonal cavity tunnels (the remaining off-axis portion of the tunnel is composed of 12 pentagonal faces), such that the x tunnels and y tunnels are slightly offset from each other because these tunnels cannot cross each other. In contrast, the z face shows 9 hexagonal cavity tunnels. Here, the hexagonal faces of the cavities form bridges between successive tunnels. Six water molecules in the

Figure 1. CH4 and CO2 hydrate stability curves based on various laboratory studies reported in the literature. Quadruple point Q1 involves CH4 hydrate equilibrium between CH4 hydrate, ice, liquid water, and CH4 vapor (hCH4-I- Lw-VCH4); quadruple point Q2 involves CO2 hydrate equilibrium between CO2 hydrate, ice, liquid water, and CO2 vapor (hCO2-I-Lw-VCO2); quadruple point Q3 involves CO2 hydrate equilibrium between CO2 hydrate, liquid water, liquid CO2 and CO2 vapor (hCO2-Lw- LCO2-VCO2).

curves for both CH4 and CO2 hydrates.12−14 Here it is seen that the CH4 hydrate slope changes at the water freezing point, but dissociated CH4 remains in the vapor phase under all conditions. In contrast, CO2 appears as hydrate, vapor, and liquid phases within a relatively small P−T envelope and slopes are also affected by water freezing. Also the CO2 stability curve crosses the CH4 stability curve at slightly above 283.15 K. These observations further clearly demonstrate that a CH4− CO2 replacement would be a very complex process which is needed for an accurate prediction. This is one motivation to conduct systematic MD simulations for CH4 and CO2 hydrates coexisting with water, but at the nano scale. 1.2. Goals and Objectives. The main objective of this paper is to explore and compare CH4 and CO2 hydrate dissociation and bubble evolution via molecular dynamics simulation. We also wish to investigate if Arrhenius temperature dependence can be demonstrated from these simulations and the role of bubble evolution on this process. With this, we hope to eventually establish a sound theoretical basis for the mechanisms and analysis of CO2 replacement of CH4 hydrate as an attractive enhanced methane gas recovery process. These ideas are illustrated schematically in Figure 2, where we will focus our analysis on the dynamic meta zone leading to intermediate phase (II). In this region, highly concentrated trapped gas and hydrogen bonds with short life expectancy coexist leading to fluctuating localized temperature gradients 1972



Article

for CH4 hydrate and 90% for CO2 hydrate observed in practice but this still provides a useful base case. English et al.18 have considered the role of partial occupancy on hydrate dissociation via MD simulation and found only a weak dependence. On the basis of this model development, all physical descriptions of the methane hydrate reported by Sloan19 were honored. More particularly, they discussed the hexagonal chain formed the backbone bridging of SI hydrate. On the basis of our physical model, we also distinguish the hexagonally bonded water molecules types (two H2O molecules connected hexagons which surrounded by four larges cages and other four H2O surrounded by one small and three large cages). These structural differences may affect dissociation behavior at the scale of molecular dynamics. 2.2. Setup MD Simulation Boxes. The CH4 hydrate case for MD simulation was generated by further adding a fluid phase to the z-face of the crystal structure. The length of each fluid phase is equal to the length of the hydrate crystal. The number of water molecules was based on a chosen water density. The CO2 hydrate case for numerical study was generated analogously by again adding a fluid phase of equal length to the hydrate crystal. These cases provide the initial conditions for the simulations which follow are shown in Figure 4. We have explored the role of adding a number of salt (Na and Cl) molecules separately, but this is not reported here. Macroscopically, this would represent a class II hydrate reservoir, where hydrate crystal coexists with saline water. 2.3. MD Simulation Methodology. These simulation models were implemented in the MD simulator, Groningen Machine for Chemical Simulations, GROMACS (www. gromacs.org). Here the three molecule types (H2O, CH4, CO2) were treated as rigid with specified interatomic lengths and angles. Several polarization models for H2O (SPC versus T1P4P) and CO2 (simple three-point versus five-point virtual site) were explored, but for the simulations reported here we settled on the T1P4P and simple three-point representations, respectively. We recognize that this choice will shift the predicted water-ice and hydrate stability curves of Figure 1 significantly to the left, i.e., to lower temperatures (see Conde and Vega20 plus numerous discussions in the GROMACS Web site). The intermolecular van der Waals forces were specified with Lennard-Jones 6-12 potentials (Table 1). Lorentz− Berthelot mixing rules were used for the LJ cross interactions between different atom types. Additionally, the long-range electrostatic forces utilized smoothed particle mesh Ewald summation. Periodic boundary conditions were chosen with a cell length of 3.414, 3.414, and 13.980 nm in x, y, and z directions. Note with assumption of periodicity in the z direction, hydrate melting at the z = 0 crystal face can cause dissociated molecules to appear near the z = L face. This will be observed in the simulations which follow. The stability of the initial hydrate crystal was first demonstrated by energy minimization followed by constant molecule number, volume, and temperature (NVT) ensemble and constant molecule number, volume, and energy (NVE) ensemble equilibration simulations, both run over 500 ps time periods. The dissociation process itself was run utilizing a constant molecule number, pressure, and temperature (NPT) ensemble. Here the temperature thermostat chosen is Nose− Hoover whereas the pressure barostat used is Parrinello− Rahman. Reflecting the known hydrate stability curves for CH4

Figure 3. Pictures showing a physical model of the SI hydrate crystal: (a) top (x−y) view and (b) site (x−z) view (overall dimension, x = 34.14 Å, y = 34.14 Å, z = 19.13 Å; 600 bonded H2O molecules and 55 empty cages). The model highlights a key structural stability element: the hexagonal chain works as a backbone in the growth of SI hydrate; the (1)H2O water molecule, connected by two hexagons oriented 90° to each other, is surrounded by four large cages, and the (2)H2O water molecule is surrounded by the three large and one small cages.

unit cell form bridges between these dodecahedra in such a way that two molecules are connected to the large cavity whereas the remainder are connected to the smaller cavities. These structural characteristics have been determined by X-ray scattering measurements by Palin and Powell,16 as further illustrated by Figure 1 of van der Waals and Platteeuw.17 These hydrate crystal characteristics are now converted to a structurally correct hydrate model for molecular dynamics simulation by generating appropriate coordinate locations for each atom in the crystal. This provides the appropriate framework for locating the guest molecules at the center of each cavity. 100% cavity occupancy is assumed throughout this model development. We recognize that this assumption is significantly different from the fractional occupancies of 95% 1973



Article

Figure 4. MD simulation systems for simulating dissociations of the CH4 and CO2 hydrates coexisting with water. The CH4 and CO2 hydrate crystals were solvated with an equal volumewater in the normal z-direction (initial box dimension, x = 3.414 nm, y = 3.414 nm, z = 13.980 nm). Here, all the cages (100% occupancy) are filled with guest CH4 and CO2 molecules.

Table 1. Partial Charges of CH4, CO2, and H2O Molecules and Atomic Site−Site Lennard-Jones Potential Parameters for CH4 Hydrate−H2O and CO2 Hydrate−H2O Systems

a

site

mass

charge q (e)

σ (A)

ε (kJ/mol)

C (CH4) H (CH4) C(CO2) O(CO2) HW(H2O)a OW(H2O) MW(H2O)

12.011 1.008 12.011 15.994 1.008 15.994 0.0

−0.24 0.06 0.70 −0.35 0.0 0.52 −1.04

3.500 2.500 3.750 2.960 0.0 3.154 0.0

0.2761 0.1255 0.4393 0.8786 0.0 0.6485 0.0

TIP4P water model.

hydrate and CO2 hydrate, the hydrate dissociation process was explored at several selected temperatures (T = 240, 260, 280, 300 K) and pressures (P = 2, 10 MPa for CH4 hydrate and P = 1, 2, 10 MPa for CO2 hydrate). Utilizing a time step size of 1 fs, simulations were typically run to 4 ns. This was felt sufficient to analyze most aspects of the dissociation process as the gas phase evolution is completed at higher temperatures, and bulk phase properties such as density, volume, and energy status are not significantly changed by the simulation end time.

Figure 5. Potential energy stabilizations for the CH4 and CO2 hydrates coexisting with water over the length of the energy minimization simulations.

crystal size effect). That there are a few more free CO2 molecules released, as compared with the CH4 hydrate case, may reflect the increased solubility of CO2 in water relative to CH4 solubility. 3.2. Analysis NVT and NVE Simulations. NVT and NVE ensemble simulations were then conducted to further establish the correctness of the initial hydrate structure. Here we show only the final structures after the NVE simulations in Figure 7. For CH4 hydrate, an additional few free CH4 molecules are generated again due to the periodic boundary conditions and only approximate stability on a bounding edge, which is a manifestation of the small crystal size necessarily utilized in these simulations, as mentioned previously. For CO2 hydrate, there are also a few free CO2 molecules released, due to the periodic boundary conditions and only approximate stability on a bounding edge. Figure 8 shows the guest molecules density profiles along the z-axis and corresponding water molecules orientation with respect normal z-axis at various temperatures, indicating even at 320 K, guest molecules remains in the hydrate and there is no breakage of the crystals. Finally, the corresponding temperature profiles are shown in Figure 9. Here the limited surface breakage of the hydrate at the water−hydrate contact is responsible for the small temperature

3.0. RESULTS OF HYDRATE STABILITY SIMULATIONS 3.1. Analysis Energy Minimization. Energy minimization is first used to establish the correctness of the initial hydrate structure and eliminate unphysical atomic motion. Figure 5 illustrates the results of the calculation and the essential stability of the proposed structures. CO2 molecules settle relatively more quickly as energy minimization achieved earlier. The magnitude of potential energy is higher in the CO2 hydrate. We feel it is imperative to establish the stability of the initial hydrate structure before dissociation simulations are conducted; otherwise, additional artificial dynamic processes can occur in the simulations. The stability of these initial structures is further demonstrated in Figure 6. For CH4 hydrate, the few free CH4 molecules generated are due to the periodic boundary conditions employed and the approximate stability obtained on a bounding edge, which is actually a manifestation of the small crystal size necessarily utilized in these simulations. For CO2 hydrate, the few free CO2 molecules generated are again due to the periodic boundary conditions employed and the approximate stability obtained on a bounding edge (i.e., small 1974



Article

Figure 6. 3D view showing the atomic configurations for the CH4 and CO2 hydrates coexisting with water at the end of the energy minimization simulations. For clarity images on the left show only the guest molecules.

Figure 7. Atomic configurations of the CH4 and CO2 hydrates at the end of 500 ps of NVE simulation (initial reference temperature and pressure, T0 = 300 K and P0 = 10 MPa).

decline. This should be compared with Figures 6 and 7 of Ripmeester et al.,8 which shows a more extensive temperature

decline and which indicates their structures are unstable and break completely. 1975



Article

Figure 8. Variations in density of guest molecules (CH4 and CO2) along the z-axis and the orientation of water molecules with respect to normal zaxis. NVE simulations of the CH4 and CO2 hydrates at the initial reference temperatures T0 = 260, 280, 300, 320 K and pressure P0 = 2, 10 MPa.

Figure 9. Variations of overall system temperature. NVE simulations of the CH4 and CO2 hydrates at the initial reference temperatures T0 = 260, 280, 300, 320 K and pressure P0 = 2, 10 MPa.

follow. We note that the number of guest molecules enclosed in the initial hydrate structure greatly exceeds the solubility limits of either CH4 or CO2 in the adjacent water (e.g., 1 m3 hydrate

In summary, our simulations unequivalently demonstrate the essentially stability of the proposed hydrate crystal structure as a proper initial state for the decomposition simulations which 1976



Article

Figure 10. CH4 and CO2 hydrate dissociations at a temperature of 240 K. Atomic configurations showing the displacements of the CH4 and CO2 molecules at the end of 4 ns of NPT simulations.


1977



Article



1978



Article

Figure 14. Variations in potential energy over the length of 4 ns NPT simulations of the CH4 and CO2 hydrate dissociations at temperatures 240, 260, 280, and 300 K.

accentuating the finite crystal size behavior. The higher temperature (T = 280, 300 K) cases show distinct hydrate dissociation behavior with eventual formation of a gas phase bubble. Because of significant solubility of CO2 in water, it is expected that a higher number of CO2 molecules remain dissolved in water at the end of the bubble formation phase. As a consequence, the CO2 bubble formed is noticeably smaller that the corresponding methane bubble (generated from equal number of hydrate guest molecules in the simulations). For both hydrate types, the calculated limits of hydrate stability are shown to occur at conditions shifted from that of the experimental bulk stability curves of Figure 1. For example, at 10 MPa, the experimental CH4 hydrate stability is 286 K (the red or Sloan curve in Figure 1) whereas our simulated stability temperature is approximately 260 K. This is a known effect20 attributed to the choice of water model for the MD simulations. This choice further implies no consideration of ice formation in our models at these temperatures. An earlier study21 using the same empirical water model T1P4P performed MD simulation of ice nucleation and growth leading to water freezing. They performed MD trajectories for different temperatures, but only temperatures lower than 230 K underwent the ice crystallization using the T1P4P water model. Some further comments about the equilibrium solubility limits of CH4 in water and CO2 in water should also be made, on the basis of these simulations. From the volume of gas bubbles formed at higher temperatures, it is clear that CH4 has lower solubility in water than CO2 as expected. More concretely, and as one example, with 315 CO2 and 5152 H2O molecules in initial MD hydrate model, we count 110 free CO2 and 205 CO2 in gas bubbles at p = 2 MPa, T = 300 K (Figure 13b). This translates to a CO2 solubility of 110/5152 = 0.019, which is consistent with experimental CO2 molar solubility. 4.2. Analysis Time Plots. The same CH4 and CO2 hydrate dissociation cases can be further explored by utilizing various time plots. Figure 14 shows the variation in potential energy over the hydrate dissociation process for both gases. These plots show the bond changing activity of potential energy evolution over time. More particularly, for CH4 hydrate the potential energy plot indicates there is no change (hydrate

dissociates into 173 m3 gas at STP), yet the initial hydrate configuration remains stable. Thus there is no “concentration” driving force to melt the hydrate under the initial conditions, and the small number of free guest molecules generated during the initialization is due strictly to a surface structure effect. The issue of gas solubility limits will be shown in the next section under P, T conditions whereby the hydrate actually dissociates.

4.0. RESULTS OF NPT HYDRATE DISSOCIATION SIMULATIONS 4.1. Analysis Spatial Plots. The CH4 hydrate dissociation process was explored at several selected temperatures (T = 240, 260, 280, 300 K) at a chosen pressure of 10 MPa, which represents a realistic pressure value to encounter CH4 hydrate. The CO2 hydrate dissociation process was explored at the same selected temperatures but at a chosen pressure of 2 MPa, which represents a realistic pressure value to encounter CO2 hydrate. These choices of conditions are consistent with the hydrate stability curves for the two gases as shown in Figure 1. Spatial distribution plots of CH4 hydrate dissociation after 4 ns of simulation at temperatures T = 240, 260, 280, 300 K are shown in Figures 10−13, respectively. Except for the finite crystal size effect mentioned above, the T = 240 K and T = 260 K structures can be seen to remain stable, with the higher temperature of temperature of 260 K possibly accentuating the finite crystal size behavior. The higher temperature (T = 280, 300 K) cases show distinct hydrate dissociation behavior with eventual formation of a gas phase bubble. We note that the hydrate breakup tends to occur from the outer crystal faces with free CH4 molecules eventually diffusing together to form a single bubble (most clearly illustrated in Figure 13). Also, because of limited solubility of CH4 in water, it is expected that only some CH4 molecules remain dissolved in water at the end of the bubble formation phase. Analogously, spatial distribution plots of CO2 hydrate dissociation after 4 ns of simulation at temperatures T = 240, 260, 280, 300 K are also shown in Figures 10−13, respectively. The T = 240 K and T = 260 K structures can be seen to remain stable (except for the finite crystal size effect mentioned above), with the higher temperature of temperature of 260 K possibly 1979



Article

formation of a gas phase filling part of the simulation volume. In particular at T = 300 K, a step change in system volume or density for CH4 hydrate at about 2500 ps is observed. This corresponds to the dramatic formation of a single large free gas bubble, as also noted in the spatial plot of Figure 13a. For CO2 hydrate we note, however, the higher initial overall density here relative to the CH4 hydrate cases, because the two hydrates differ significantly in density. Dynamically for CO2 hydrate, there is again no change (dissociation) at temperature of 240 K, whereas at all three higher temperatures, dissociation appears to have stabilized after 4 ns. The significantly lower densities (higher volumes) achieved at longer times at the three higher temperatures reflect the formation of a gas phase filling part of the simulation volume. In contrast to CH4 hydrate at T = 300 K, no dramatic step change in system volume or density for CO2 hydrate is observed. This results from the much lower volume of the single free gas bubble seen here, as also noted in the spatial plot of Figure 13b. These observations are further reflected in the average CH4 and CO2 density profiles along the z-axis at 4 ns as shown in Figure 18 for the four temperature cases. We recall that at initial time, CH4 and CO2 only exist in the hydrate in half the simulation box. Because of the hydrate crystal structure, this results in a somewhat oscillatory pattern in half the box and zero in the remaining portion. At 4 ns and T = 240 K, Figure 18 demonstrates that this pattern is essentially unchanged, whereas at T = 260 K, this structure is beginning to evolve as the hydrate starts to dissociate. At the two highest temperatures at 4 ns, a CH4 or CO2 gas bubble has formed in a significant region of the simulation box, and the profiles shown in Figure 18 reflect this state. More particularly, these curves reflect the spatial distributions of gas molecules shown in Figures 10−13 at various temperatures. As the temperature increases, the hydrate structure begins to break and the excess gas above their water solubility limits, forms gas bubbles of varying sizes, shapes, and locations as shown. There is larger CH4 gas bubble because of its lower solubility in water. Water molecule orientation profiles, as shown in Figure 19 at 4 ns similarly reflect the changing state of the system when the CH4 or CO2 hydrate dissociates. Again the initial state, with

decomposition) at temperature of 240 K, whereas complete decomposition appears to have occurred by 4 ns at the two highest temperatures of 280 and 300 K. Only at T = 260 K is the dissociation process still ongoing at 4 ns. For CO2 hydrate, the potential energy plot indicates there is no change (hydrate decomposition) at temperature of 240 K, whereas complete decomposition appears to have occurred by 4 ns at all three highest temperatures of 260, 280, and 300 K. Figure 15 shows the predicted system kinetic energy after 4 ns for each hydrate system as a function of temperature. The

Figure 15. Kinetic energy as a function temperature over the length of 4 ns NPT simulations of the CH4 and CO2 hydrate dissociation.

kinetic energy is seen to increase linearly with temperature, whereas the difference in absolute energy can be attributed to the difference in molecular masses of CH4 and CO2. Figures 16 and 17 show the variation in system density and system volume (i.e., inverse properties) over the CH4 and CO2 dissociation process to 4 ns for the four temperature cases. The initial volume/density states are very similar at the four temperatures because the water and hydrate phases are show small variation with temperature. Dynamically for CH4 hydrate, there is no change (dissociation) at a temperature of 240 K, whereas at 260 K dissociation appears to underway. The significantly lower densities (higher volumes) achieved at longer times at the two higher temperatures reflect the

Figure 16. Change in the system volume of the CH4 and CO2 hydrate dissociation over the length of 4 ns of NPT simulations at temperatures 240, 260, 280, and 300 K. 1980



Article

Figure 17. Change in the system density of the CH4 and CO2 hydrate dissociations over the length of 4 ns of NPT simulations at temperatures 240, 260, 280, and 300 K.

Figure 18. Variation in the density of CH4 and CO2 molecules along the z-axis at the end of 4 ns of NPT simulations at temperatures 240, 260, 280, and 300 K.

ordered water molecules in the hydrate crystal and disordered (random) water molecules in the adjacent water bath, is essentially retained after 4 ns for the T = 240 K case. As temperature increases and the hydrate begins to melt, the net water orientation tends to a random profile throughout the simulation box. The radial distribution function profiles for CH4−H2O and CH4−CH4 pairs also reflect the local structure existing during the simulation. Typically, these give useful information on the local structure in a globally homogeneous system. But to interpret these plots in our case, we must remember these are calculated averages over the whole simulation box, whereas the actual limiting global distributions tend to be inhomogeneous: hydrate in one-half of the box and aqueous fluid in the other, or a large gas bubble separating aqueous fluid. However, useful information can still be found. Figure 20a shows such plots at 4 ns for the four temperature cases on interest. The T = 240 K case, where negligible hydrate dissociation occurs even after 4 ns essentially illustrates the oscillating pattern of CH4

distribution within the hydrate. As the hydrate dissociates (e.g., at higher temperatures), this oscillating pattern disappears and a single large gas bubble is created. Then the CH4−H2O radial distribution function averages to very small values for most distances as the only CH4−H2O contact is at the boundary of a large gas bubble. There is limited contribution of CH4−H2O in the aqueous phase because of the low CH4 solubility. Conversely, the CH4−CH4 radial distribution function peaks at very high values, again reflecting one large pure CH4 gas bubble. The radial distribution function profiles for CO2−H2O and CO2−CO2 pairs also reflect the local structure existing during the simulation. Figure 20b shows such plots at 4 ns for the four temperature cases on interest. The T = 240 K case, where negligible hydrate dissociation occurs even after 4 ns essentially illustrates the oscillating pattern of CO2 distribution within the hydrate. This oscillating pattern disappears and a single large gas bubble is created as the hydrate dissociates at higher temperatures. Because of the increased solubility of CO2 in 1981



Article

Figure 19. Variation in the orientation of water molecules with respect to the normal (z-axis) of the simulation box at the end of 4 ns of NPT simulations of the CH4 and CO2 hydrate dissociations at temperatures 240, 260, 280, and 300 K.

Figure 20. Variation in the radial distribution functions for pairs of CH4−H2O and CH4−CH4 and of CO2−H2O and CO2−CO2 at the end of 4 ns of NPT simulations for the CH4 and CO2 hydrate dissociation at temperatures 240, 260, 280, and 300 K.

4.3. Effects of Pressure. Throughout the EM, NVE, and NVT ensemble simulations we assessed the role of pressure at

water, relative to CH4, these changes to the radial distribution functions are less dramatic than for the CH4 case. 1982



Article

Figure 21. Effect of pressure on the system volume over the length of 4 ns of NPT simulations of the CH4 and CO2 hydrate dissociations at temperature 260 K and pressures 1, 2, and 10 MPa.

Figure 22. Effect of pressure on system density over the length of 4 ns of NPT simulations of the CH4 and CO2 hydrate dissociations at temperature 260 K and pressures 1, 2, and 10 MPa.

coefficients which are reported in Table 2. For comparison, experimental diffusion coefficients for CH4 and CO2 in water over a temperature range of 273−308 K have been reported by Maharajh and Walkley22 and compared against several empirical relationships. Here both molecules have been found to have similar values and with a slight increase with temperature (1.0 × 10−5 cm2/s and larger). Their reported diffusion values over the ranged from 1.05 × 10−5 to 1.96 × 10−5 cm2/s for CH4−H2O and 1.00 × 10−5 to 2.41 × 10−5 cm2/ s for CO2−H2O. Zeebe23 has used mean square displacement analysis from molecular dynamics to calculate CO2 diffusion in water and its temperature dependence, again obtaining similar values. Our results show more extreme temperature behavior, but our system is quite different and this explains our observations. At low temperature without hydrate breakage, the effective diffusion coefficient represents an average of gas diffusion in water along with limited gas molecule motions in the hydrate

various temperatures while preserving the hydrate structure. The final pressure value was 10 MPa at this point. We have also conducted some limited analysis of the role of pressure by repeating the CH4 hydrate simulations at P = 2 MPa (base case 10 MPa), and the CO2 hydrate simulations at pressures of 1 and 10 MPa (base case 2 MPa). Although we will not show many of these results explicitly, we have found a limited effect of pressure, and visually all results shown above were repeated at different pressure levels. As an example, Figures 21 and 22 show system volume and density variations with pressures over time at T = 260 K. The equivalent variation of density along the z-axis at 4 ns as a function of pressure is shown in Figure 23.

5.0. ANALYSIS KINETIC BEHAVIOR 5.1. Diffusion Coefficient. Plots of mean square displacements for CH4 and CO2 show essentially linear behavior over 2 ns simulations at any temperature (actual plots not shown). These can be analyzed to generate effective diffusion 1983



Article

Figure 23. Effect of pressure on the CH4 and CO2 molecules density along the z-axis at the end of 4 ns of NPT simulations of the CH4 and CO2 hydrate dissociations at temperature 260 K and pressures 1, 2, and 10 MPa.

5.3. Hydrate Dissociation Kinetics. Further quantitative comparison of the dissociation process can be established by analyzing the breakage of hydrogen bonds during dissociation for both guest molecules. Essentially we will convert the number of hydrogen bonds broken per time, to the amount of guest molecules released per time. Several additional temperature simulations were run and analyzed to give a more complete picture of the temperature behavior. The rate of hydrogen bonds breaking increases with an increase in temperature, and the activation energy is positive. This process occurs only if a sufficient amount of energy is made available to permit the vibration between molecules to become so vigorous as to overcome the van der Waals forces and to break the hydrogen bonds. For the structure I CH4 hydrate, for example, this process can be defined by the simple reaction CH4·nH2O → CH4 + nH2O + ΔH. Here CH4 is the trapped guest molecule within the hydrogen bonded cages by the H2O molecules, and ΔH is the endothermic heat of reaction. The hydration number, n, is approximately n = 5.75, when all cages are filled with guest molecules. As each oxygen atom associated with 4 hydrogen bonds, one can assume approximately 23 (i.e., 4 × 5.75) hydrogen bonds needed to break for one CH4 molecule to escape. More concretely, prediction of gas hydrate dissociation rates is calculated as follows: Determine the H-bond breaking rates [H bond/ps = (initial number of hydrogen bonds − number of hydrogen bonds in completely dissociated hydrates)/dissociation time] of hydrogen bonds associated in the gas hydrate crystals. Assume 23 hydrogen bonds need to break for one guest molecule to escape from a completely filled gas hydrate. The gas molecule escaping rate is now expressed in molecules/ps. This is converted to the molar gas escaping rate, [(g mol)/ps] by use of Avogadro’s number [6.023 × 1023 molecules/(g mol)]. Recognizing the hydrate unit cell dimension is 1.2 nm, a surface dissociation rate [mol/(nm2 ps)] can be calculated by dividing by the unit cell surface area of 1.4 nm2. Finally, using a pressure difference of 1 MPa between the set pressure and hydrate equilibrium pressure for each gas hydrate (Figure 1), the upscaled dissociation rate in mol/(m2 MPa s) is obtained, where we have also upscaled time and distance units. The choice of 1 MPa pressure difference during the NPT ensemble

Table 2. Predicted Effective Diffusion Coefficients of the Guest Molecules, CH4 and CO2 temp (K) 240 250 260 265 270 280 290 300

effective CH4 diff coeff D × 10−5 effective CO2 diff coeff D × 10−5 (cm2/s) (cm2/s) 0.0294 0.1093 0.3081 0.5462 0.7270 1.7405 4.5013 4.7342

(±0.0010) (±0.0074) (±0.0423) (±0.0834) (±0.3876) (±0.6052) (±0.3117) (±0.6842)

0.3212 0.3918 0.7191 0.7269 0.9484 0.9622 1.1636 1.8476

(±0.0169) (±0.0803) (±0.0382) (±0.0045) (±0.0484) (±0.0118) (±0.0323) (±0.0065)

crystal. At high temperatures with hydrate breakage, the effective diffusion coefficient represents an average of gas diffusion in water along with extensive gas molecule motions in a free gas bubble. Thus our effective diffusion coefficients can bracket the expected basic diffusion behavior of gas molecules in water. 5.2. Hydrogen Bond Breakage Behavior. Analysis of hydrogen bond behavior will prove useful in our analysis of dissociation kinetics. Hydrogen bonds are determined from a GROMACS feature based on cutoffs for the angle hydrogendonor−acceptor and the distance donor−acceptor. Here OH is viewed as a donor, whereas O is an acceptor, and the cutoff values employed are cutoff radius 0.35 nm and cutoff angle 30°. Figure 24 shows a spatial plot of the number of active hydrogen bonds at the end of 4 ns for both CH4 and CO2 hydrate cases at temperatures of 240 and 260 K. The regular structure of the Hbonds in the region of hydrate crystal can be seen, whereas in the liquid water region the H-bonds are more randomly distributed. Also observable are the absence of H-bonds in the region of free gas formation at T = 260 K. Equivalent plots at two higher temperatures, 280 and 300 K (not shown) have complete H-bond breakage in the now larger region of the free gas bubbles formed. Figure 25 illustrates the time evolution on the number of hydrogen bonds at various temperatures during both the CH4 hydrate and CO2 hydrate dissociation. As noted above, the CO2 hydrate shows more rapid hydrogen bond breakage at all temperatures. 1984



Article

Figure 24. Spatial plots showing the hydrogen bonds at the end of 4 ns of NPT simulations for the CH4 and CO2 hydrate dissociation at temperatures 240 and 260 K.

Figure 25. Plots showing the number of hydrogen bonds over the length of 4 ns of NPT simulations for the CH4 and CO2 hydrate dissociation at temperatures 240, 260, 280, and 300 K.

shifts the curve to the left from the actual experimental curve (Figure 1). But because the overall shape of the shifted hydrate stability curves is not altered, our predicted hydrate dissociation rates are independent of the chosen water model. Thus the choice of TIP4P is sufficient for our purposes over the complete temperature range considered. Following this analysis, the intrinsic breakage rate at various temperatures is plotted versus 1/T as in Figure 26. This method captures the relative rates of hydrate dissociation for

simulations as the pressure driver is somewhat arbitrary but is based on Kim et al.’s experiments, which we were all near this value (and close to the hydrate stability curve under all conditions). Kim et al.’s fugacity is an effective pressure. Our sensitivity studies indicate that hydrate dissociation is more sensitive to temperature than to pressure (see section 4.3), as also found by Kim et al.7 experimentally. Here we have recognized by Conde and Vega20 that the equilibrium stability curve depends on the employed water model (TIP4P), which 1985



Article

reflect this microscopically, as all effects are ultimately determined by the input CO2 force field. Note as the activation energies represent relative rates as a function of temperature, these values can be considered as more reliable than the actual magnitudes of the intrinsic rate constant (i.e., less assumptions are required to establish a numerical value). Any quoted activated energies should be viewed as fitted values over the range of temperatures studied, and not extrapolated to temperatures beyond this range. As discussed by Ripmeester et al.,8 activation energies much larger than heats of formation may be a result of the large scale rearrangements of hydrate, water, and gas phases during the decomposition of a row of solid hydrate phase. As a point of reference, Horvat et al.24 quote heats of formation of CH4 hydrate of 54.49 and 57.98 kJ/mol for CO2 hydrate (i.e., very similar values). Our methane activation energy is comparable with English and Phelan’s18 activation energy value 58.3 kJ/mol from their planar MD model, and the methane hydrate heat of formation, Makogon.25 These values are somewhat lower than Kim et al.’s7 and Clarke and Bishnoi’s9 experimental activation energy of 81 kJ/mol. The MD simulated result of Myshakin et al.26 is consistent with their value. For CO2 hydrate dissociation, our overall activation energy value is consistent with the range of values calculated by Sarupria and Debenedetti27 from MD modeling (62.3 and 69.3 kJ/mol using two different evaluation methods). Our higher temperature activation energy is consistent with the CO2 hydrate heat of formation from Anderson.28 Clarke and Bishnoi obtained an experimental value of 102.9 kJ/mol, which is closer to our lower temperature CO2 activation energy. It is important to note that the choice of activation energy value is determined from the temperature dependence of the dissociation rate, as illustrated above. As the pressure dependence of dissociation is very weak relative to temperature, pressure results are averaged for this calculation, as was also done with laboratory dissociation data (Kim et al.7 for methane hydrate and Clarke and Bishnoi9,10 for CO2 hydrate). We are aware that our relative results for intrinsic rate constants are quite different from the macroscopic fitting results of Clarke and Bishnoi9,10 comparing CH4 and CO2 hydrate dissociation. They quote intrinsic rate constants for CO2 hydrate which are 4 orders of magnitude greater than that for CH4 hydrate, whereas our intrinsic rates have a 2 order of magnitude variation. However, this is based on their single choice for CO2 activation energy, which is significantly higher than that chosen for CH4 hydrate, which in turn affects their calculated intrinsic rate constant. Our differences in activation

Figure 26. Effect of temperature on the dissociation rate constants for the CH4 and CO2 hydrates. Values in parentheses are the activation energy in kJ/mol.

each (CH4 and CO2) hydrate at each temperature. These demonstrate Arrhenius-like behavior. Interestingly, the slopes of the CH4 hydrate dissociation and CO2 hydrate dissociation are similar, which indicate approximately equal activation energies. This observation is consistent with the rest of the results we have observed from our MD modeling, which showed relatively similar hydrate dissociation behavior for the two hydrates. The values of kd given in Table 2 are plotted in Figure 26 against the inverse of temperature. The plot shows a liner fit, suggesting the following Arrhenius-type equation for kd: kd = kd0e−Ea / RT

(1)

From the best linear fit in Figure 26, the activation energy was obtained as Ea = (R ln(k2/k1)/(1/T1 − 1/T2)), Here, k1 and k2 are the dissociation (hydrogen bonds breaking) rates at any two selected temperatures T1 and T2, respectively, along the fitted line. An optimum value of k0 was then obtained by substituting Ea into the above eq 1 at any given temperature. The optimum values obtained are given in Table 3. Some additional comments on these results should be made. We have demonstrated that the relationships (dissociations rates vs temperatures) for both gas hydrates are Arrhenius (i.e., ln (kd) vs 1/T is linear, Figure 26). However, two different slopes are noticeable in the CO2 hydrate dissociation. We obtained activation energies and intrinsic rates for these two dissociation regimes as well as a best fit overall value. Interestingly, the CO2 hydrate stability curve also shows a change in slope behavior due to liquid versus gas CO2 formation (Figure 1). It is possible that our MD simulations

Table 3. Simulated H-Bond Breaking Rates and Predicted Intrinsic Dissociation Rates of CH4 and CO2 Hydrates temp (K)

H-bond breaking rate khb (no./ps)

270 280 290 300

0.288 0.750 1.714 4.000

1.250 3.261 7.453 17.391

× × × ×

10−2 10−2 10−2 10−2

270 280 290 300

1.667 2.750 6.250 14.000

5.072 11.957 27.174 60.870

× × × ×

10−2 10−2 10−2 10−2

molecule escaping rate k = (khb/23) (molecule/ps)

molar surface rate kd =k(1/AS)(1/NA) (mol/nm2.ps)

CH4 Hydrate 1.441 3.769 8.595 20.055

× × × ×

10−26 10−26 10−26 10−26

14.410 37.690 85.951 20.055

× × × ×

103 103 103 104

5.849 13.788 31.336 70.194

× × × ×

10−26 10−26 10−26 10−26

58.495 13.788 31.336 70.194

× × × ×

103 104 104 104

CO2 Hydrate

1986

upscale intrinsic rate kd [mol/(m2 MPa s)]



Article

Table 4. Correlated Activation Energy and Intrinsic Dissociation Rate for the CH4 and CO2 Hydrate Dissociations CH4 hydrate

CO2 hydrate

kinetic parameters

T > 260 K

240 K < T < 260 K

300 K > T > 260 K

over all T range

activation energy, Ea (kJ/mol) intrinsic rate, k0d [mol/(m2 MPa s)]

59.8 5.100 × 1015

86.4 4.690 × 1021

56.9 7.767 × 1015

67.8 6.000 × 1017

clearly demonstrated gas bubble evolution mechanisms as well as quantifying the implications of the dissociation stages (Figure 2 schematic), including the drastic changes in overall system density, volume changes during gas phase evolution process, etc. Our simulations have established the essential similarities between CH4 hydrate and CO2 hydrate dissociation. The limiting behavior at lower temperature (no dissociation) and higher temperature (complete melting and formation of a gas bubble) have been illustrated for both hydrates. Due to the shift in the known hydrate stability curves between guest molecules, the intermediate behavior (e.g., 260 K) showed distinct differences however. Also, because of the more hydrogenbonding capability of CO2 in water, as reflected in its molecular parameters, higher solubility of dissociated CO2 in water was observed with a consequence of a smaller size of gas bubble formation. Additionally, we have proposed a novel method for analyzing hydrate dissociation based on H-bond breakage and used this to quantify the dissociation behaviors of both CH4 and CO2 hydrates. Activation energy, Ea, values from our MD studies were obtained (Table 4). These values were evaluated against several other published laboratory and MD values. Intrinsic rate constants were estimated and upscaled. We have also tried to indicate some macroscopic consequences of this analysis by proposing a kinetic reaction model consistent with macroscale fitted kinetic models we have developed previously.29 Further extensions of this work are possible. In particular, the role of partially filled cavities and the presence of various levels of salinity are of interest, as each represent factors in the Mallik gas hydrate reservoir we have studied at a macroscopic (field scale) modeling. Somewhat more challenging, or at least more computationally intensive, is to address hydrate formation modeling via molecular dynamics. Here we also compare the behavior of CH4 hydrate formation with CO2 hydrate formation using the same molecular parameters. The role of bounding sand (e.g., silica interfaces) on both hydrate dissociation and formation processes is also of interest, as all hydrate reserves are found encased in a porous media. As mentioned in the Introduction, the ultimate goal of this research is an improved understanding of mechanisms leading to CO2 replacement in CH4 hydrates, and how such a process could be economically achieved at the field scale (e.g., at the Mallik natural gas hydrate deposit). During primary production (reduction of pressure to cause hydrate dissociation), the CH4−CO2 replacement process may not be practical. But this process could be very attractive when primary production declines substantially. The main reason is that there is a significant amount of free CH4 gas in the reservoir at the end of primary production. That will allow a good mobility of CO2 and CH4 gases in the reservoir. On the basis of Figure 1, as CO2 partitions into different phases easily and has high solubility in water, CO2 mobility can be significantly reduced and trapped into the reservoir and subsequently form hydrate and at the same time CH4 can be produced as it remains in vapor phase. In our earlier paper,30 we postulated this idea as a

energies are smaller, and these affects are calculated intrinsic rate constant values. Their actual experimental dissociation rates for both hydrates are similar, as are our simulated dissociation rates. 5.4. Proposed Kinetic Reaction Model. On the basis of details MD simulation studies, we propose the following gas bubble kinetic reaction model for, e.g., CH4 hydrate dissociation in bulk water (CO2 hydrate dissociation would employ the same scheme but different parameters). kd

CH4 ·nH 2O → xCH4(d) + (1 − x)CH4(b) + nH 2O + ΔH kT

CH4(b) → CH4(g)

(1) (2)

Reaction 1 is the primary first-order hydrate dissociation reaction at rate of kd (1/time). Here dissociated guest CH4 molecules from hydrate, CH4·nH2O, partitioned into dissolved gas, CH4(d) of fraction x and remaining fraction (1 − x) into dispersed gas bubbles, CH4(b). The dissolved gas fraction, CH4(d), satisfies the water phase solubility requirements and does not participate in further gas dynamics. Reaction 2 is the secondary reaction to transfer, CH4(b), into free gas phase, CH4(g), at a rate of kT (1/time) and represents the ultimate connection of various isolated gas bubbles into a bulk gas phase. Here no heat of reaction for bubble coalescence is assumed. This step is not simulated in our current MD model due to the large spatial and time scales required but does provide a framework for the further evolution of gas. Our molecular dynamics simulations have generated comparative values for solubility, x, and dissociation rates, kd, for CH4 and CO2 hydrates. More complex multistage gas evolution models can be envisioned, but with many possibly short-lived intermediate components and (possibly ill-defined) reaction parameters. Because of the low viscosity of the continuous water phase, these intermediates have been neglected in our current simple model. This framework is quite similar to our macroscale gas bubbles evolution model in a porous media used to fit laboratory and field gas exsolution from hydrates,29 where the model parameters were obtained via sensitivity studies. The above MD gas bubble in bulk water represents a small gas bubble trapped within pore throat in a porous media flow system. In this sense, the fitted MD reaction model is subset (the first step) of our earlier multistep reaction model at short space and time scales.

6.0. CONCLUSIONS We have presented molecular dynamics simulations of gas hydrate dissociation comparing the behavior of CH4 and CO2 hydrates. These simulations were based on a theoretical hydrate crystal, generating two identical MD simulation boxes for CH4 and CO2 hydrates coexisting with water. The MD simulation system was first stabilized via a thorough EM, NVT, and NVE simulations before proceeding to NPT simulation for targeted dissociation pressure and temperature responses. Here we 1987



Article

(13) Yoon, J.-H; Chun, M.-K.; Lee, H. Generalized model for predicting phase behavior of clathrate hydrate. AIChE J. 2002, 48 (no. 6), 1317−1330. (14) Diamond, L. W.; Akinfiev, N. N. Solubility of CO2 in water from −1.5 to 100 °C and from 0.1 to 100 MPa: evaluation of literature data and thermodynamic modelling. Fluid Phase Equilib. 2003, 208, 265− 290. (15) Lu, H.; Seo, Y.; Lee, J.; Moudrakovski, I.; Ripmeester, J.; Chapman, N.; Coffin, R.; Gardner, G.; Poulman, J. Complex gas hydrate from the Cascadia margin. Nature 2007, 445, 303−306. (16) Palin, D. E.; Powell, H. M. Nature 1945, 156, 334; J. Chem. Soc. (London) 1947, 208; 1948, 571. (17) van der Waals, J. H.; Platteeuw, J. C. Clathrate solutions. Adv. Chem. Phys. 2007, DOI: 10.1002/9780470143483.ch1. (18) English, N. J.; Phelan, G. M. Molecular dynamics study of thermal-driven methane hydrate dissociation. J. Chem. Phys. 2009, 131, 074704. (19) Sloan, E. D. Physical/chemical properties of gas hydrates and application to world margin stability and climate change. In Gas Hydrate: Relevance to World Margin Stability and Climate Change; Henriet, J.-P., Mienert, J., Eds.; Geological Society: London 1998; Vol. 137, pp 31−50, Special Publications. (20) Conde, M. M.; Vega, C. Determining the three-phase coexistence line in methane hydrates using computer simulations. J. Chem. Phys. 2010, 133, 064507. (21) Matsumoto, M.; Saito, S.; Ohmine, I. Molecular dynamics simulation of the ice nucleation and growth process leading to water freezing. Nature 2002, 416, 409−413. (22) Maharajh, D. M.; Warzinski, J. The Temperature Dependence of the Diffusion Coefficients of Ar, CO2, CH4, CH3Cl, CH3Br, and CHCl2F in Water. Can. J. Chem. 1973, 51, 944−952. (23) Zeebe, R. E. On the molecular diffusion coefficients of dissolved CO2, HCO3-, and CO3-2 and their dependence on isotopic mass. Geochim. Cosmochim. Acta 2011, 75, 2483−2498. (24) Horvat, K.; Kerkar, P.; Jones, J.; Mahajan, D. Kinetics of the formation and dissociation of gas hydrates from CO2-CH4 mixtures. Energies 2012, 5, 2248−2262. (25) Makogon, Y. F. Hydrates of natural gas, W. J. Cieslewicz Translation; Pennwell: Tulsa, OK, 1981. (26) Myshakin, E. M.; Jiang, H.; Warzinski, R. P.; Jordan, K. D. Molecular dynamics simulations of methane hydrate decomposition. J. Phys. Chem. A 2009, 113, 1913−1921. (27) Sarupria, S.; Debenedeth, P. G. Molecular dynamics study of carbon dioxide hydrate dissociation. J. Phys. Chem. A 2011, 115, 6102− 6111. (28) Anderson, G. K. Enthalpy of dissociation and hydration number of methane hydrate from the Clapeyron equation. J. Chem. Thermodyn. 2004, 36, 1119−1127. (29) Uddin, M.; Wright, F.; Coombe, D. Numerical study of gas evolution and transport behaviours in natural gas-Hydrate reservoirs, J. Can. Pet. Technol. JCPT, 2011, vol. 50, no. 1. (30) Uddin, M.; Wright, F.; Coombe, D. Modelling of CO2 hydrate formation in geological reservoirs by injection of CO2 gas, J. Energy Resour. Technol. 2008 (JERT-07-1100).

macroscopic scheme, and which recognized the crucial role of gas bubble formation in this process.

■

AUTHOR INFORMATION

Corresponding Author

*M Uddin: phone, +1 780 450 5047; fax, +1 780 450 5242; email, mafi[email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors thank the Alberta Innovates − Technology Futures (AITF) and the Computer Modelling Group Ltd. (CMG) for their financial and technical support. Our extensive interactions and discussions with Fred Wright and Scott Dallimore of Geological Survey of Canada (NRCan) formed the motivation for this work, and they are gratefully acknowledged.

■

REFERENCES

(1) Dallimore, S. R.; Collett, T. S.; Uchida, T. In Scientific results from JAPEX/JNOG/GSC Mallik 2L-38 gas hydrate research well, Mackenzie Delta, Northwest Territories, Canada; Dallimore, S. R., Uchida, T., Collett, T. S., Eds.; Geological Survey of Canada, Bulletin 544, 1999. (2) Dallimore, S. R.; Collett, T. S. In Scientific results from the Mallik 2002 gas hydrate production research well program, Mackenzie Delta, Northwest Territories, Canada; Dallimore, S. R., Collett, T. S., Eds.; Geological Survey of Canada, Bulletin 585, 2005. (3) Dallimore, S. R.; Yamamoto, K.; Wright, J. F.; Bellefleur, G. In Scientific results from JOGMEC/NRCan/Aurora Mallik 2007−2008 gas hydrate production research well program, Mackenzie Delta, Northwest Territories, Canada; Dallimore, S. R., Yamamoto, K., Wright, J. F., Bellefleur, G., Eds.; Geological Survey of Canada, Bulletin 601, 2012. (4) Uddin, M.; Wright, F.; Dallimore, S. R.; Coombe, D. Gas hydrate production from the Mallik reservoir: numerical history matching and long-term production forecasting, In Scientific results from JOGMEC/ NRCan/Aurora Mallik 2007−2008 gas hydrate production research well program, Mackenzie Delta, Northwest Territories, Canada; Dallimore, S. R., Yamamoto, K., Wright, J. F., Bellefleur, G., Eds.;Geological Survey of Canada, Bulletin 601, 2012, p 261−289. (5) Uddin, M.; Wright, F.; Dallimore, S. R.; Coombe, D. Gas hydrate dissociations in Mallik hydrate bearing zones A, B, C by depressurization: effect of salinity and hydration number in hydrate dissociation. J. Pet. Sci. Eng. 2013, submitted for publication. (6) Uddin, M.; Wright, F.; Dallimore, S. R.; Coombe, D. Seismic correlated Mallik 3D gas hydrate distribution: effect of geomechanics in non-homogeneous hydrate dissociation by depressurization. J. Pet. Sci. Eng. 2013, submitted for publication. (7) Kim, H. C.; Bishnoi, P. R.; Heidemann, R. A.; Rizvi, S. S. H. Kinetics of methane hydrate decomposition. Chem. Eng. Sci. 1987, 42 (No. 7), 1645−1653. (8) Ripmeester, J. A.; Alireza, S.; Hosseini, B.; Englezos, P.; Alavi, S. Fundamentals of methane hydrate decomposition. CSUG/SPE 138112, presented at the Canadian Unconventional Resources and International Petroleum Conference, October 19−21, 2010, Calgary, AB, Canada. (9) Clarke, M. A.; Bishnoi, P. R. Determination of the activation energy and intrinsic rate constant of methane gas hydrate decomposition. Can. J. Chem. Eng. 2001, 79, 143−147. (10) Clarke, M. A.; Bishnoi, P. R. Determination of the intrinsic rate constant and activation energy of CO2 gas hydrate decomposition using in-situ particle size analysis. Chem. Eng. Sc. 2004, 59, 2983−2993. (11) Clarke, M. A.; Bishnoi, P. R. Measuring and modelling the rate of decomposition of gas hydraes formed from mixtures of methane and ethane. Chem. Eng. Sc. 2001, 56, 4715−4724. (12) Adisasmito, S.; Frank, R. J.; Sloan, E. D. Hydrates of carbon dioxide and methane mixtures. J. Chem. Eng. Data 1991, 36, 68−71. 1988


Kinetics of CH4 and CO2 Hydrate Dissociation and Gas Bubble

Recommend Documents