Heat Transfer Calculations for Decomposition of Structure I Methane

Publication Date (Web): May 17, 2013 ... Estimates of the properties at the macroscopic scale, like the equilibrium temperature of methane hydrate and...
0 downloads 0 Views 2MB Size
Download PDF
Article pubs.acs.org/JPCC

Heat Transfer Calculations for Decomposition of Structure I Methane Hydrates by Molecular Dynamics Simulation Vikesh Singh Baghel,† Rajnish Kumar,*,† and Sudip Roy*,‡ †

Chemical Engineering & Process Development and ‡Physical Chemistry Division, CSIRNational Chemical Laboratory, Pune-411008, India S Supporting Information *

ABSTRACT: Microcanonical ensemble molecular dynamics simulations of structure I methane hydrate is presented in this work to study the endothermic decomposition process. The mechanism of decomposition of methane hydrate as a function of time was explained at the molecular level. The initial temperature and pressure of the simulation were chosen so as to depict the natural gas hydrate in conditions of oceanic sediments. A more realistic strategy was developed to perform the microcanonical ensemble simulation of solid−liquid interface of hydrate and amorphous water. Two water models, SPC/E and TIP4P, were used for the simulations, and the results of the simulations were compared. Heat transfer calculations were performed on the adiabatic system, and an attempt has been made to fit the MD simulation results to the heat balance equations derived from the heat transfer calculations. Estimates of the properties at the macroscopic scale, like the equilibrium temperature of methane hydrate and rate of supply of hot water for sustained release of methane from solid hydrate phase, were determined. The equilibrium temperature obtained by the above method was found to be in agreement with the experimentally observed value. Both the SPC/E and TIP4P water models gave similar results.

I. INTRODUCTION Methane clathrate hydrates (methane hydrates/gas hydrates) belong to a class of inclusion compounds in which water molecules encapsulate the methane molecules.1 The three most common structures of gas hydrates are cubic structure I (SI), cubic structure II (SII), and hexagonal structure H (SH). The SI unit cell has two 512 (pentagonal dodecahedron) and six 51262 (hexagonal truncated trapezohedron) cages; similarly, the SII unit cell has sixteen 512 and eight 51264 cages, and SH has three 512, two 435663, and one 51268 cages. All the cages have one cavity each. From the thermodynamic point of view, the most favorable structure of methane hydrate is SI.1 In this paper, we have considered only SI methane gas hydrate. Though gas hydrates are considered as a nuisance in the oil and gas industries, research groups have found various novel applications of gas hydrates, like hydrogen storage and carbon dioxide sequestration.2−6 The discovery of naturally occurring methane gas hydrates has paved the way for the exploration of a new potential energy source in this era of depleting conventional fossil fuels. Methane gas hydrates are known to exist in the permafrost regions and in the oceanic sediments.7 Methane recovered from such a resource can go a long way in satisfying the need of a clean fuel (comparatively less environmental concern) having a high calorific value.7 Extraction of methane from natural gas hydrates, however, is not trivial, as its chemical and physical properties are very sensitive to pressure and temperature conditions, which are attributed mostly to © XXXX American Chemical Society

their complex structures and nonbonded interactions between the host and guest molecules. It is therefore, very difficult to predict pressure and temperature conditions for a safe and continuous supply of methane gas from a natural gas hydrate reservoir. It is important to note that uncontrolled release of methane from natural gas hydrates can increase the risks of disasters like earthquakes, tsunamis, and environmental hazards like global warming.8 Therefore, it is necessary to understand the mechanism of hydrate decomposition at the molecular level to get insight of structure collapse and related change in properties of gas hydrates. Also, as atomic and molecular level interactions between the guest and host molecules play a significant role, a systematic understanding of the decomposition process will provide the prospect to understand the dynamics of carbon dioxide sequestration and storage in natural sediments. A suitable and widely used technique to carry out such mechanistic study at the molecular level is classical molecular dynamics (MD) simulations. The Newton’s equations of motions are solved to find the trajectories of the particles, which can be processed and analyzed to get several macroscopic properties that are experimentally observed. MD simulations have successfully been used to study various biological,9 polymeric,10,11 and gas hydrate12−17 systems and have given a much deeper insight toward the mechanism and structure−property relations in these fields of study. Received: March 8, 2013 Revised: May 17, 2013

A

dx.doi.org/10.1021/jp4023772 | J. Phys. Chem. C XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry C

Article

Figure 1. Simulation box. Bulk water (WSOL) represents the water molecules that surround the SI fully occupied hydrate cages. Green molecules are the WSOL molecules, blue molecules are the Whyd molecules, gray represents the ML molecules, and yellow represents the Ms molecules.

Methane hydrate decomposition has been studied experimentally and theoretically in the past. Experimentally, the kinetics of the methane hydrate decomposition was measured by Kim et al.18 and they have concluded that the decomposition kinetics is dependent on the temperature, pressure, and particle surface area. They have correlated decomposition rates and respective temperatures via the Arrhenius equation, from which the resulting activation barrier is obtained as the heat of decomposition of the hydrate. To address the molecular level mechanism of decomposition, atomistic simulations were also performed with different sizes of simulation box and using different ensembles. Myshakin et al.14 simulated in the NPT ensemble the decomposition mechanism of the SI structure as a function of different cage occupancies with methane gas. They reported that because of empty cages the hydrate structure gets destabilized. English et al.19 used the same method but with a polarizable water force field to understand hydrate decomposition in the presence of liquid water and a mixture of water and methane. They have pointed out that the diffusion of methane molecules to bulk phase is the rate-determining step for decomposition. It is well-known that the hydrate decomposition is a nonisothermal process, as at the molecular level it involves dissociation and rearrangement of the hydrogen-bonding network. The hydrate decomposition process is endothermic,20 and the pressure of the system rises as the gas expands into the gaseous phase from the solid hydrate phase. Therefore, to simulate the real process of hydrate recovery from natural gas hydrate, simulations should be carried out under constant volume and constant energy (NVE) conditions. Considering the above facts, Alavi and Ripmeester21 performed NVE simulations and were able to observe the drop in temperature and rise in pressure during the decomposition of SI gas hydrates with 100% small (512) and large (51262) cage occupancy. They did a short NVT and NPT equilibration of 200 ps each with constrained hydrate molecules, followed by 5 ps of heating before performing the 450 ps NVE simulations. Performing NVE simulation is not trivial, because probable accumulation of numerical error arises due to the absence of a thermostat and barostat. Use of a thermostat and barostat helps by scaling the temperature and pressure and also takes care of the numerical errors in velocity and position of atoms that arise due to numerical integration of equation of motion. In this work, we propose a new strategy to perform NVE simulations to understand methane hydrate decomposition. To create an initial structure of solid−liquid interface for simulation,

we have used the experimental crystal structure of methane hydrate and bulk water system. Further, we have optimized the crystal structure of methane hydrate to the potential energy surface of SPC/E and TIP4P water models at 273 K and 10 MPa (experimentally reported temperature and pressure where the methane hydrate and liquid water falls in to the stable region of phase diagram).22 Interfacial liquid water was also optimized (equilibrated) so that the water molecules at the interface have minimized velocity gradients. Simulated annealing of the above system was done to increase its temperature in very small steps to the desired initial temperature of interest for studying the decomposition process. This procedure takes care of the proper optimization of the methane hydrate and bulk water system to the potential energy surface of water models used in this work. We believe that our system has a well-optimized solid−liquid interface and hence has correct initial velocity and positions for each atom to run NVE simulation. To further validate the application of the simulation procedure to macroscopic scale, an attempt was made to fit the data obtained from the MD simulations to equations developed from a macroscopic heat balance done on the simulation system. The equilibrium temperature obtained from the fit was found to be in agreement with the experimental results. To further extract useful macroscopic quantities, the rate of supply of hot water for a constant supply of methane gas was calculated, which is a very useful quantity for developing a technology for the methane extraction process from methane hydrates by thermal stimulation.23 The results obtained by using SPC/E and TIP4P water models were compared and found to be similar. The details of computational methods and heat transfer calculations are discussed in the following sections of this paper.

II. COMPUTATIONAL METHOD The unit cell of methane hydrate was obtained from the available literature.24,21 A 3 × 3 × 6 supercell of SI methane hydrate was placed in a simulation box with water on both the sides of hydrate, as shown in Figure 1. For tracking the molecules precisely, the water molecules forming hydrate cages (Whyd) were labeled differently from the bulk water outside the hydrate structure (WSOL), and the methane molecules present in the small cages (Ms) were labeled differently from that in large cages (Ml). The system consisted of 108 Ms molecules, 324 Ml molecules, 2484 Whyd molecules, and 3027 WSOL molecules. The final simulation box was of the dimension of 3.5286 × 3.5286 × 14.5793 nm3. The crystalline symmetry was preserved B

dx.doi.org/10.1021/jp4023772 | J. Phys. Chem. C XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry C

Article

along x and y directions by applying periodic boundary conditions in all the x, y, and z directions. Double-precision GROMACS25 (version 4.5.5) software package was used to perform simulations. Two well-known water models, SPC/E and TIP4P, were used. Alavi and Ripmeester21 showed that the tetrahedral rigid methane molecules of C−H bond length of 1.0895 Å and bond angle of 104.45° with nonbonded parameters taken from Murad-Gubbins potential26 showed satisfactory results to understand the decomposition process. In the current work, the same force field parameters for methane were used, and the standard Lorentz− Berthelot combination rule was applied. Potential energy functions were used as implemented in GROMACS25 with nonbonded Lennard-Jones potentials and Columbic interactions. The leapfrog algorithm was used to integrate the equations of motion with a time step of 1 fs for NVT and NPT simulations (equilibration steps) and 0.5 fs for NVE simulations (production run). The smaller simulation time step for NVE production run was chosen to reduce the error in the total energy calculations and hence to maintain a constant total energy. Shake algorithm was used to constrain the bond length between the atoms in rigid molecules (all molecules were taken as rigid) with a relative shake tolerance of 1 × 10−10. Short-range nonbonded interactions between the molecules were calculated within a cutoff distance (Rcut) of 15.0 Å. Particle mesh Ewald (PME) summations were used for calculating the Columbic long-range interactions with a relative Ewald tolerance of 1 × 10−6. Simulation Procedure. The pressure was chosen as 10 MPa to mimic the oceanic sediments conditions at a depth of 1000 m, where the gas hydrates are known to exist.27 The temperature was chosen as 273 K; at 10 MPa and 273 K, methane hydrate falls in the stable region of the phase diagram.22 The system was subjected to equilibration run at 273 K and 10 MPa. The equilibration was not trivial due to the presence of solid and liquid phases, as huge velocity gradients would occur at the liquid−solid interface, which could destabilize the system. The NVT simulations were performed for 2 ns with a position restrain (as implemented in GROMACS25) applied on the hydrate phase atoms (Whyd, Ms, and Ml) at 273 K with a Berendsen thermostat,28 which had a coupling time constant of 0.5 ps to prevent huge oscillations of temperature, as the system may be far from equilibrium. The output structure file of the NVT simulation was used as the input structure for the 1 ns NPT simulations with position restrain applied on the hydrate phase at 10 MPa pressure. By doing this, the WSOL molecules (bulk water) were equilibrated at the desired temperature and pressure. In the next step, the position restrain on the hydrate atoms (all atoms of Whyd, Ms, and Ml) were taken off and NPT simulations were run for 500 ps at the 273 K and 10 MPa with starting structure taken as the output structure of the previous NPT run. Due to the stepwise equilibration as reported above, the entire system was satisfactorily equilibrated at 273 K and 10 MPa pressure to the potential energy surface of the classical force field used for the simulation. A Nose−Hoover thermostat29,30 with coupling time constant of 0.5 ps and a Parrinello Rahman barostat31,32 with a coupling time constant of 2 ps were used for both the restrained and unrestrained NPT simulations. The output structure and velocities of atoms after equilibration were taken as input for the NVE production run. The typical simulation times for NVE runs were 2 ns. To study the decomposition process at the higher temperatures of 290, 300, and 310 K, simulated annealing was performed only to increase the temperature gradually from 273 K to higher temperatures with the velocities of the molecules and input structure

taken from the last step of equilibration as an input. The output structure and velocities of molecules after the simulated annealing step were taken as input for the NVE production runs. Both the Lennard-Jones and Columbic interaction potentials smoothly decayed to zero in the case of NVE simulations with the switch starting at a distance of 1 Å before the Rcut.25 A generalized simulation flowchart (Figure 2) shows the above-mentioned procedure stepby-step.

Figure 2. Simulation flowchart.

The methane molecules that came out of the open cages present at the hydrate−bulk water interface during equilibration were not considered for calculating the rate of hydrate decomposition. All the simulations were repeated five times, independent of previous run, to obtain properties with maximum and minimum errors. MD simulations typically provide a deep knowledge about the decomposition mechanism of methane hydrates at the molecular level. However, for developing a thermal stimulation technology for methane recovery from natural gas hydrate, it is necessary to use the knowledge gained from MD simulations to compliment macroscopic scale results. In this work, an attempt was made to fit the simulation results to the equations developed from macroscopic heat transfer equations (discussed in the following section) that were scaled so as to remove the effect of length and time. Although the limitations of scaling up of time and length scale are well-known, the results obtained from the fit gave an equilibrium temperature that is in agreement with the experimental value. This gave us the confidence to further calculate the rate of heating water to be supplied for sustained release of methane gas from the reservoir. The details of the heat transfer calculations are discussed in the next section.

III. HEAT TRANSFER CALCULATIONS The simulation strategy used in this work is designed to depict the methane extraction process by thermal stimulation in a reservoir. For a constant supply of methane gas from the reservoir, it is necessary to estimate the rate of hot water supply (represented as bulk water in the simulation box). The heat lost by the bulk water to the solid hydrate phase for endothermic decomposition should be replenished by the external hot water supply. A constant heat supply is necessary for sustained decomposition of solid hydrate phase. The rate of heat to be C

dx.doi.org/10.1021/jp4023772 | J. Phys. Chem. C XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry C

Article

Figure 3. Temperature and pressure profile of the system at different initial temperatures as a function of time. Plots a and b show temperature and pressure profiles, respectively, for the SPC/E water model, while c and d represent temperature and pressure profiles, respectively, for the TIP4P water model.

by hydrate is used for the supply of methane gas. Therefore, we have

supplied was equated to the rate of heat lost by the bulk water, which is given by eq 1 m w0CV (T0 − Tt ) dm ws = (T0 − Tref )CV tA s dt

− (1)

(3)

where dN/dt is the number of methane molecules moved out of the hydrate per unit time, L is the rate of accumulation of heat, and K is an arbitrary quantity that is independent of time. Equation 3 can be rearranged as

where mw0 is the initial mass of bulk water, mws is the additional mass of hot water to be supplied, CV is the specific heat of water at constant volume [taken as 84 J/(mol K)33 for the SPC/E water model and 75.3 J/(mol K)19 for the TIP4P water model], T0 is the initial temperature of the bulk water, t is the time, Tt is the temperature of the bulk water at time t, Tref is the reference temperature (which is taken as 273 K), and As is the cross sectional area of methane hydrate available for heat transfer from bulk water to the methane hydrate. Since mws = ρwVws, dVws m w0(T0 − Tt ) = dt (T0 − Tref )tA sρw

m w0CV dT dN L +K T= As dt dt As

−

dT = A + BT dt

(4)

where A=

L m w0CV

B=

−KA s dN m w0CV dt

(5)

or (2)

⎞ ⎛ m w0CV ⎟ ⎜ K = − B ⎜ dN ⎟ ⎝ dt A s ⎠

where ρw is the density of water and Vws is the volume of water to be supplied. The density of water was taken as 998 kg/m3 for the SPC/E water model34 and 1001 kg/m3 for the TIP4P water model.35 To get a better understanding about the heat lost by the bulk water, the temperature profile of bulk water was calculated as follows: As the bulk water and the hydrate system were adiabatic, the heat lost by the bulk water should be equal to the heat gained by the solid methane hydrate for decomposition. By conservation of energy, the rate of heat lost by the bulk water along with the rate of heat required for methane extraction should be equal to the rate of accumulation of heat. The rate of heat required for methane extraction is assumed to be a function of the rate of methane gas coming out of the hydrate cages as the heat gained

(6)

On integrating eq 4 we have

∫T

T

−

0

dT = A + BT

∫0

t

dt

(7)

The rate of accumulation of heat in the hydrate system is mainly the latent heat of fusion of the hydrate. As our aim is to achieve a constant rate of methane supply by providing a constant heat to the system, the rate of accumulation of heat can be assumed to be independent of time. Hence, we have from eq 7 T=− D

⎞ A ⎛⎜ A + + T0⎟e−Bt ⎝B ⎠ B

(8)

dx.doi.org/10.1021/jp4023772 | J. Phys. Chem. C XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry C

Article

Figure 4. Four-body order parameter F4 of hydrate-forming Whyd water molecules as a function of time. The curves on the left represent the SPC/E water model, while the curves on the right represent the TIP4P water model. The inset on the left panel represents the schematic diagram for calculation of the four-body structural order parameter F4.

outermost hydrogen atoms. Outermost hydrogen atoms were selected by calculating distances between hydrogen atoms and the neighboring oxygen atom. The average value of F4 for hydrate system and water is 0.7 and −0.04, respectively.37 Hence, when the hydrate decomposes from hydrate phase to liquid phase, the value of F4 for Whyd molecules also decreases. Figure 4 shows a decreasing trend in the F4 value for both SPC/E and TIP4P water models at initial temperatures of 300 and 310 K. It can be seen in Figure 3 that for the simulation runs at initial temperatures of 300 and 310 K, the rate of temperature drop reduces as the temperature of the system reaches around 290 K, which is close to the experimentally observed equilibrium temperature of 285 K at 10 MPa pressure.22 This result is similar for both SPC/E and TIP4P water models; this behavior can be attributed to the reduced driving force for the endothermic decomposition of hydrates. Insufficient driving force prevents further decomposition as the temperature of the system reaches close to the equilibrium temperature, and therefore, the rate of decrease of temperature reduces. Thus, the decomposition process is found to be heat transfer limited, as also indicated by Ullerich et al.20 and Alavi and Ripmeester.21 To understand the mode of heat transfer during hydrate decomposition, the simulation box was divided into five regions (see Figure 5, regions 1 and 5 are bulk water, regions 2 and 4 are interfaces between bulk water and hydrate phase, and region 3 is the hydrate phase), and detailed analyses were done to measure the temperature variations in the bulk and at the gas/hydrate interface, as shown in Figure 5. The average kinetic energy of all the atoms present in each region was calculated. The temperature was calculated as 2KEavg/Ndfk, where KEavg is the average kinetic energy of the atoms in the region, k is the Boltzmann constant, and Ndf is the number of degrees of freedom. To represent the trend in the temperature drop, the running average of every 50 ps of simulation time (average of 100 frames, we saved frames at 0.5 ps intervals) was calculated and is represented in red in Figure 5. As the process of decomposition is endothermic, the hydrate phase takes the heat from the surrounding bulk water for decomposition. Figure 5 shows the temperature drop in all five

At steady state, the bulk water will achieve equilibrium temperature, Teq. So, as t → ∞, T → Teq, and we have Teq = −

A B

(9)

On substituting the values of A and B in eq 8, the temperature profile of bulk water is obtained as ⎛ KA s ⎛ dN ⎞ ⎞ ⎜ ⎟t ⎟ T = Teq + (T0 − Teq) exp⎜ ⎝ m w0CV ⎝ dt ⎠ ⎠

(10)

The simulation data obtained in this work is fitted in the eqs 2 and 10 as reported in the Results and Discussion section of this paper.

IV. RESULTS AND DISCUSSION Figure 3 shows the temperature and pressure variations in the system for both SPC/E and TIP4P water models at different initial temperatures, which have the same trend as observed by Alavi and Ripmeester.21 To observe the pressure changes clearly, running averages of every hundred pressure versus time data points obtained from the simulations were plotted. At 273 K and 10 MPa, the methane hydrate system does not undergo decomposition, and as expected, no temperature and pressure variation is observed during the simulation run, while at temperatures higher than the equilibrium temperature of 285 K and 10 MPa,22 the endothermic methane hydrate decomposition is accurately captured by the temperature drop and simultaneous rise in the pressure of the system. Decomposition studied at higher driving force36 shows a comparatively higher rate of temperature drop and pressure rise. To further confirm the decomposition at 300 and 310 K, the four-body structural order parameter, F4 for Whyd molecules was calculated as a function of time. F4 is defined as (1/n)∑ni=1cos 3ϕi, where n is the number of oxygen−oxygen pairs of water molecules present within 0.3 nm of distance from each other and ϕ is the torsion angle between them (inset, Figure 4). The angle ϕ is calculated for neighboring oxygen atoms with their E

dx.doi.org/10.1021/jp4023772 | J. Phys. Chem. C XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry C

Article

Figure 5. Temperature in the different regions of the simulation box as defined by the upper panel as a function of time. In the lower panel the upper and bottom curves are corresponding to SPC/E and TIP4P water models, respectively.

methane molecules during 1950−2000 ps confirms that the innermost cages in the simulation box are still intact. As the simulation time elapses, the rate of methane gas coming out of the hydrate cages, i.e., the kinetics of decomposition of methane hydrates, also decreases (see Figure 9, discussed below). This is possibly because of the decrease in temperature (see Figure 5) as a function of time in all the regions of the systems, which is unable to provide enough thermal energy for further decomposition. Now the question is whether the hydrate decomposition is row-by-row. To get a better understanding of this mechanism, we calculated the percentage of number of Whyd molecules moved out of cage from each of the six unit cells (rows) present along the z axis with respect to time (Figure 8). We have chosen 0.5 nm distance as a criterion for the percentage of water molecules going out of the cage; 0.5 nm signifies two factors at a time, the movement of water molecules by slightly greater than half of the cage diameter and crossing two solvation shells of water [as shown in the radial distribution of oxygen atoms in the hydrate cage (Figure S1 in the Supporting Information)]. This criterion guarantees that the water molecules have moved out of the cage. While methane hydrate dissociates, there is a possibility of positional drift of the system. Therefore, for calculating the displacement of molecules (water and methane, as depicted in Figures 8 and 9, respectively) the amount of drift is compensated by translating the hydrate phase at every 0.5 ps. Figure 8 shows that the decomposition of methane hydrate starts from the hydrate−bulk water interface and the decomposition of the unit cell located away from the interface starts only when at least 30− 50% of the interface unit cell had decomposed. This observation shows that the decomposition of hydrate cages is not entirely row by row. After decomposition starts in one row, the next row is

regions for both SPC/E and TIP4P water models. The initial temperature was 310 K and pressure was 10 MPa. Similar results were observed for SPC/E and TIP4P water models. The temperature drop in all the five regions reaches close to 290 K; at this point the rate of temperature drop decreases with time. When the hydrate decomposes, the molecules in the interfacial region (regions 2 and 4) move more (because decomposition starts at the interface), resulting in higher fluctuations in the temperature in this region as compared to region 3 (hydrate phase). These results gave a clear understanding of the heat transfer taking place in the system. However, the hydrate decomposition process is also known to be mass transfer limited.21 Alavi and Ripmeester21 have showed that significant mass transfer resistance is offered to the decomposition process as almost an entire row of methane molecules comes out of the hydrate phase. Figures 6 and 7 of our calculations show the mechanism of this decomposition process. The left sides of Figures 6 and 7 show the snapshots of methane molecules (both Ms and ML) at different simulation times for SPC/E and TIP4P water models, respectively, at an initial simulation temperature of 310 K and 10 MPa pressure. It can be seen that almost an entire row decomposes at about 300 and 500 ps. The right sides of Figures 6 and 7 depict the partial densities (averaged over 50 ps of each interval and all five trajectories) of methane molecules as a function of box length along the z direction at different simulation times. From these figures, we observe that the methane density decreases in the interface regions, which suggests that decomposition starts at the hydrate−bulk water interface. However, the structure does not decompose completely after the 2 ns of simulation, as shown in the snapshots (see Figures 6 and 7). The partial density of F

dx.doi.org/10.1021/jp4023772 | J. Phys. Chem. C XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry C

Article

Figure 6. The left panel shows the snapshots of methane molecules present in the simulation box at different simulation times for the SPC/E water model at an initial simulation temperature of 310 K and 10 MPa pressure. The right panel shows the partial density of methane molecules as a function of the z direction of the box length at different simulation times for the SPC/E water model.

Figure 7. The left panel shows the snapshots of methane molecules present in the simulation box at different simulation times for the TIP4P water model at an initial simulation temperature of 310 K and 10 MPa pressure. The right panel shows the partial density of methane molecules as a function of the z direction of box length at different simulation times for the TIP4P water model. G

dx.doi.org/10.1021/jp4023772 | J. Phys. Chem. C XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry C

Article

Figure 8. Percentage of hydrate-forming water molecules (Whyd) present in the six different crystal unit cells (rows along the z direction) of methane as a function of time: (left) SPC/E and (right) TIP4P water models. The inner row of hydrate-forming water molecules (Whyd) starts to dissociate when almost 30−50% of the Whyd molecules present at the hydrate−bulk water (WSOL) interface has dissociated.

Figure 9. Number of methane molecules moved out of the cages as a function of time at different initial temperatures: (left) SPC/E and (right) TIP4P water models.

The rate of methane release was calculated for different initial temperatures. Any methane molecule that had displaced more than 0.8 nm in either the x, y, or z direction from its starting position was considered to be out of the cage.21 The number of methane molecules moved out of the cages as a function of time at different initial temperatures for both SPC/E and TIP4P water models is shown in Figure 9. It shows that the rate increases with increase in temperature as more kinetic energy is provided to the molecules to overcome the thermodynamic, mass transfer, and heat transfer barriers.

also affected. This is because the water molecules are making the cage by forming hydrogen bonds to its neighboring molecules located in the nearest rows. So if one row starts to decompose, the row next to it also gets affected. It also clearly shows that the decomposition process is mass transfer limited. Both SPC/E and TIP4P water models showed similar results. The heat transfer and mass transfer resistances affect the kinetics of the methane extraction process, which is defined by the rate at which methane gas comes out of the methane hydrate. H

dx.doi.org/10.1021/jp4023772 | J. Phys. Chem. C XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry C

Article

Figure 10. Temperature profile of bulk water (WSOL) as a function of time. The upper curve represents the values for an initial temperature of 300 and 310 K, respectively: (left) SPC/E and (right) TIP4P water models. The blue lines represent the equilibrium temperature as observed from the temperature profile of WSOL molecules.

Now, we correlate the molecular level understanding with the macroscopic properties in the following way. The average temperature profile of the bulk water (WSOL) molecules present on both the sides of methane hydrate was calculated as a function of time. Figure 10 shows the temperature profile of WSOL molecules at initial temperatures of 300 and 310 K for SPC/E and TIP4P water models. The fluctuations in the temperature as a function of time, which is averaged out by calculating the running average of every 100 data points, are shown in Figure 10. The running average data was fitted to eq 10 to calculate the equilibrium temperature. For the initial simulation temperature of 300 K (and 10 MPa), the equilibrium temperatures of methane hydrate−bulk water system for SPC/E and TIP4P water models were found to be 288.74(±3.37) and 289.5(±4.24) K, respectively. While for initial simulation temperature of 310 K (and 10 MPa), the equilibrium temperatures for SPC/E and TIP4P water models were found to be 292.94(±1.85) and 289.26(±2.52) K, respectively. The experimentally reported equilibrium temperature at 10 MPa pressure for methane hydrate is 285 K.22 Since, the volume of our simulation box was kept constant; the pressure of the system would increase as methane gas moves out from solid hydrate phase and goes to liquid water phase. This is because methane molecule will be less soluble in liquid water, being hydrophobic in nature. This will give rise to a repulsive interaction between water and methane and, hence, the rise in pressure. The resulting equilibrium temperature at the elevated pressure should therefore be higher than 285 K, as observed in our calculations. The equilibrium temperatures calculated by both the SPC/E and TIP4P water models fall in the same temperature range, and hence, in this case, too, both water models showed similar results. The equilibrium temperature gave us confidence to further estimate the rate of hot water supply for constant methane

extraction rate as per eq 2. In order to apply the MD simulation data to the macroscopic scale, the length scale and time scale should be considered. As we know that the heat transfer rate is inversely proportional to the cross sectional area available for heat transfer, the rate of supply of hot water for sustained supply of methane gas should be calculated per unit interface area to remove the effect of the length scale. Also, as time is on both the right- and left-hand sides of eq 2, the effect of time scale is removed. The temperature of WSOL at time t (Tt) was taken as the running average of the last 100 data points of Figure 10. The calculated values of rate of supply of hot water at 300 K was found to be 1.54(±0.36) and 1.75(±0.15) m3/s m2 for SPC/E and TIP4P water models, respectively. At 310 K, the rate of supply of hot water was found to be 0.96(±0.14) and 1.17(±0.26) m3/s m2 for SPC/E and TIP4P water models, respectively. As expected, the rate of supply of hot water at 310 K is observed to be lower than the rate of supply for water at 300 K in our calculations. However, these calculated rates of hot water supply should be validated by conducting suitable laboratory experiments. To further remove the effects of length and time scales, experiments should be carried out to check the ratio of the hot water supply at 300 K to the water supply at 310 K. This ratio comes out to be 1.67(±0.61) and 1.58(±0.36) for SPC/E and TIP4P water models from our calculations.

V. CONCLUSIONS Double-precision GROMACS25 software package is used to study the decomposition mechanism of fully occupied SI methane hydrate by employing NVE MD simulations using SPC/E and TIP4P water models. A novel strategy for a gas hydrate system under NVE simulating conditions is developed and discussed in detail. This simulation procedure tunes the I

dx.doi.org/10.1021/jp4023772 | J. Phys. Chem. C XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry C

■

methane hydrate and water system to match the experimental equilibrium data. The temperature drop and pressure rise during the methane hydrate decomposition process was observed by performing NVE simulations. As expected, the decomposition starts from the hydrate−bulk water interface; however, the decomposition of the second and subsequent rows of the unit cell (from the interface toward solid hydrate phase) is not instantaneous, and it only dissociates once 30−50% of the interface unit cell has already decomposed. The bulk water supplies the heat required for the melting of solid hydrate phase and subsequent removal of methane. The rate of decomposition decreases as the temperature of the bulk water approaches close to the equilibrium temperature. In the absence of a continuous supply of external hot water, the decomposition process is more like an adiabatic system, and further hydrate decomposition will cease once the temperature of bulk water reaches the equilibrium temperature of the gas hydrate. The decomposition process is thus found to be mass transfer and heat transfer limited. Basic heat transfer calculations were done on the adiabatic system, and equations useful at the macroscopic scale were derived. The results obtained from the MD simulations were fitted to heat transfer equations with assumptions of a constant supply of heat and hence constant methane production rate. The equilibrium temperatures thus calculated were found to be close to experimentally observed values. The rate of supply of hot water for a constant supply of methane was reported, which needs to be further validated by laboratory experiments. All the results reported in this work showed similar behavior for both SPC/E and TIP4P water models. Since, the melting temperatures of ice for SPC/E and TIP4P water models are very similar (SPC/E, 215 K; TIP4P, 232 K) at 1 atm of pressure38 and our calculations were done at high pressure (i.e., 10 MPa), the melting phase equilibria is similar for both the water models. As TIP4P is a four-site model, it requires a longer simulation time and memory than the three-site SPC/E water model. Therefore, the SPC/E water model for studying the adiabatic decomposition of methane hydrates by MD simulations may remain as a good choice as water force field.

■

REFERENCES

(1) Sloan, E. D.; Koh, C. A. Clathrate Hydrates of Natural Gases, 3 ed.; CRC Press, Taylor & Francis Group: Boca Raton: FL, 2008. (2) Kumar, R.; Klug, D. D.; Ratcliffe, C. I.; Tulk, C. A.; Ripmeester, J. A. Low-Pressure Synthesis and Characterization of Hydrogen-Filled Ice Ic. Angew. Chem. 2013, 125 (5), 1571−1574. (3) Kumar, R.; Linga, P.; Ripmeester, J.; Englezos, P. Two-Stage Clathrate Hydrate/Membrane Process for Precombustion Capture of Carbon Dioxide and Hydrogen. J. Environ. Eng. 2009, 135 (6), 411−417. (4) Kumar, A.; Sakpal, T.; Linga, P.; Kumar, R. Influence of Contact Medium and Surfactants on Carbon Dioxide Clathrate Hydrate Kinetics. Fuel 2013, 105 (0), 664−671. (5) Florusse, L. J.; Peters, C. J.; Schoonman, J.; Hester, K. C.; Koh, C. A.; Dec, S. F.; Marsh, K. N.; Sloan, E. D. Stable Low-Pressure Hydrogen Clusters Stored in a Binary Clathrate Hydrate. Science 2004, 306 (5695), 469−471. (6) Lu, H.; Wang, J.; Liu, C.; Ratcliffe, C. I.; Becker, U.; Kumar, R.; Ripmeester, J. Multiple H2 Occupancy of Cages of Clathrate Hydrate under Mild Conditions. J. Am. Chem. Soc. 2012, 134 (22), 9160−9162. (7) Makogon, Y. F. Natural Gas HydratesA Promising Source of Energy. J. Nat. Gas Sci. Eng. 2010, 2 (1), 49−59. (8) Maslin, M.; Owen, M.; Betts, R.; Day, S.; Dunkley Jones, T.; Ridgwell, A. Gas Hydrates: Past and Future Geohazard? Philos. Trans. R. Soc., A 2010, 368 (1919), 2369−2393. (9) Pandey, P. R.; Roy, S. Headgroup Mediated Water Insertion into the DPPC Bilayer: A Molecular Dynamics Study. J. Phys. Chem. B 2011, 115 (12), 3155−3163. (10) Pahari, S.; Choudhury, C. K.; Pandey, P. R.; More, M.; Venkatnathan, A.; Roy, S. Molecular Dynamics Simulation of Phosphoric Acid Doped Monomer of Polybenzimidazole: A Potential Component Polymer Electrolyte Membrane of Fuel Cell. J. Phys. Chem. B 2012, 116 (24), 7357−7366. (11) Pandey, P. R.; Roy, S. Early Stages of Unwinding of Zwitterionic α-Helical Homopolymeric Peptides. Chem. Phys. Lett. 2011, 514 (4−6), 330−335. (12) Bai, D.; Chen, G.; Zhang, X.; Wang, W. Microsecond Molecular Dynamics Simulations of the Kinetic Pathways of Gas Hydrate Formation from Solid Surfaces. Langmuir 2011, 27 (10), 5961−5967. (13) Walsh, M. R.; Koh, C. A.; Sloan, E. D.; Sum, A. K.; Wu, D. T. Microsecond Simulations of Spontaneous Methane Hydrate Nucleation and Growth. Science 2009, 326 (5956), 1095−1098. (14) Myshakin, E. M.; Jiang, H.; Warzinski, R. P.; Jordan, K. D. Molecular Dynamics Simulations of Methane Hydrate Decomposition. J. Phys. Chem. A 2009, 113 (10), 1913−1921. (15) Iwai, Y.; Nakamura, H.; Arai, Y.; Shimoyama, Y. Analysis of Dissociation Process for Gas Hydrates by Molecular Dynamics Simulation. Mol. Simul. 2009, 36 (3), 246−253. (16) Qi, Y.; Ota, M.; Zhang, H. Molecular Dynamics Simulation of Replacement of CH4 in Hydrate with CO2. Energy Convers. Manage. 2011, 52 (7), 2682−2687. (17) Geng, C.-Y.; Wen, H.; Zhou, H. Molecular Simulation of the Potential of Methane Reoccupation during the Replacement of Methane Hydrate by CO2. J. Phys. Chem. A 2009, 113 (18), 5463−5469. (18) Kim, H. C.; Bishnoi, P. R.; Heidemann, R. A.; Rizvi, S. S. H. Kinetics of Methane Hydrate Decomposition. Chem. Eng. Sci. 1987, 42 (7), 1645−1653. (19) English, N. J.; Macelroy, J. M. D. Atomistic Simulations of Liquid Water Using Lekner Electrostatics. Mol. Phys. 2002, 100 (23), 3753− 3769. (20) Ullerich, J. W.; Selim, M. S.; Sloan, E. D. Theory and Measurement of Hydrate Dissociation. AIChE J. 1987, 33 (5), 747−752. (21) Alavi, S.; Ripmeester, J. A. Nonequilibrium Adiabatic Molecular Dynamics Simulations of Methane Clathrate Hydrate Decomposition. J. Chem. Phys. 2010, 132 (14), 144703−8. (22) Goel, N. In Situ Methane Hydrate Dissociation with Carbon Dioxide Sequestration: Current Knowledge and Issues. J. Petrol. Sci. Eng. 2006, 51 (3−4), 169−184. (23) Gabitto, J. F.; Tsouris, C. Physical Properties of Gas Hydrates: A Review. J. Thermodyn. 2010, 2010.

ASSOCIATED CONTENT

S Supporting Information *

The O−O radial distribution function and total energy of the system for different initial temperature as a function of time for SPC/E and TIP4P water models. This material is available free of charge via the Internet at http://pubs.acs.org.

■

Article

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected] (S.R.); [email protected] (R.K.). Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors would also like to acknowledge the financial support from Department of Science and Technology (project code GAP294026) for this work. V.S.B. would like to thank CSIR, India for providing the scholarship for his graduate studies. S.R. gratefully acknowledges the Center for Excellence in Scientific Computing, NCL (project code CSC0129), for financial support and computational time. The authors sincerely thank Dr. Saman Alavi for his valuable guidance and providing the starting structure file that was suitably modified for use in the GROMACS software package. J

dx.doi.org/10.1021/jp4023772 | J. Phys. Chem. C XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry C

Article

(24) Mak, T. C. W.; McMullan, R. K. Polyhedral Clathrate Hydrates. X. Structure of the Double Hydrate of Tetrahydrofuran and Hydrogen Sulfide. J. Chem. Phys. 1965, 42 (8), 2732−2737. (25) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4 (3), 435−447. (26) Murad, S.; Gubbins, K. E. Molecular-Dynamics Simulation of Methane Using a Singularity-Free Algorithm. Abstr. Pap. Am. Chem. Soc. 1978, 175 (Mar), 3−3. (27) Collett, T. S. Energy Resource Potential of Natural Gas Hydrates. AAPG Bull. 2002, 86 (11), 1971−1992. (28) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81 (8), 3684−3690. (29) Nose, S. A Unified Formulation of the Constant Temperature Molecular Dynamics Methods. J. Chem. Phys. 1984, 81 (1), 511−519. (30) Hoover, W. G. Canonical Dynamics: Equilibrium Phase−Space Distributions. Phys. Rev. A 1985, 31 (3), 1695−1697. (31) Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52 (12), 7182−7190. (32) Nosé, S.; Klein, M. L. Constant Pressure Molecular Dynamics for Molecular Systems. Mol. Phys. 1983, 50 (5), 1055−1076. (33) Park, S.-M.; Nguyen, P. H.; Stock, G. Molecular Dynamics Simulation of Cooling: Heat Transfer from a Photoexcited Peptide to the Solvent. J. Chem. Phys. 2009, 131 (18), 184503−10. (34) Mark, P.; Nilsson, L. Structure and Dynamics of the TIP3P, SPC, and SPC/E Water Models at 298 K. J. Phys. Chem. A 2001, 105 (43), 9954−9960. (35) Jorgensen, W. L.; Jenson, C. Temperature Dependence of TIP3P, SPC, and TIP4P Water from NPT Monte Carlo Simulations: Seeking Temperatures of Maximum Density. J. Comput. Chem. 1998, 19 (10), 1179−1186. (36) Linga, P.; Kumar, R.; Englezos, P. Gas Hydrate Formation from Hydrogen/Carbon Dioxide and Nitrogen/Carbon Dioxide Gas Mixtures. Chem. Eng. Sci. 2007, 62 (16), 4268−4276. (37) Moon, C.; Hawtin, R. W.; Rodger, P. M. Nucleation and Control of Clathrate Hydrates: Insights from Simulation. Faraday Discuss. 2007, 136 (0), 367−382. (38) Vega, C.; Sanz, E.; Abascal, J. L. F. The Melting Temperature of the Most Common Models of Water. J. Chem. Phys. 2005, 122 (11), 114507−9.

K

dx.doi.org/10.1021/jp4023772 | J. Phys. Chem. C XXXX, XXX, XXX−XXX