Comparison of the Properties of Xenon, Methane, and Carbon Dioxide

Nov 11, 2009 - a systematic comparison of the structural and dynamical properties for these three ... For comparison, simula- ...... hydrates under Co...
1 downloads 0 Views 3MB Size
Download PDF
J. Phys. Chem. C 2010, 114, 5555–5564

5555

Comparison of the Properties of Xenon, Methane, and Carbon Dioxide Hydrates from Equilibrium and Nonequilibrium Molecular Dynamics Simulations† H. Jiang‡ and K. D. Jordan* National Energy Technology Laboratory, U.S. Department of Energy, P.O. Box 10940, Pittsburgh, PennsylVania, and Department of Chemistry, UniVersity of Pittsburgh, Pittsburgh, PennsylVania 15260 ReceiVed: July 5, 2009; ReVised Manuscript ReceiVed: October 3, 2009

Molecular dynamics simulations are used to characterize the hydrates of Xe, methane, and CO2, allowing for a systematic comparison of the structural and dynamical properties for these three hydrates. Although the host-guest interaction energy for the T ) 0 K structures is most attractive in the case of Xe, other structural and dynamical properties from the simulations indicate that, in fact, host-guest coupling is most important for the CO2 hydrate. Specifically, the host lattice of CO2 hydrate expands more with increasing temperature than do the lattices of the xenon and methane hydrates, and the translational and rotational dynamics of the water molecules are predicted to be most perturbed in the CO2 hydrate. The simulations predict that the CO2 and xenon hydrates have lower speed of sound values and lower themal conductivities than methane hydrate or the empty lattice. Introduction Clathrate hydrates are solid compounds in which water cages enclathrate guest atoms or small nonpolar molecules.1 Methane hydrate, in particular, has attracted considerable interest due to its formation in natural gas pipelines and to the discovery that methane hydrate deposits on the bottom and beneath ocean floors and in permafrost regions contain most of the natural gas on Earth.2 Because methane is about a factor of 20 more potent greenhouse gas than CO2,2 large-scale decomposition of these deposits would have a serious environmental impact. Recovery of methane from hydrate deposits for commercial use may also be stimulated by the injection of CO2, which could also form hydrate from the water released when the methane hydrate decomposes and thus preserve the mechanical stability of the affected sediments.3,4 Numerous experimental5-12 and computational12-22 studies have been carried out for methane hydrate, and its structural, dynamical, and thermal properties are all well characterized. Much less work has been done on CO2 hydrate. Although the thermal expansivity23-26 and decomposition temperature27-31 of CO2 hydrate have been determined experimentally, we are unaware of any experimental or theoretical results for its dynamical or thermal transport properties. The radial distribution functions of CO2 hydrate have been calculated by Chialvo et al.22 and Ferdows et al.32 using molecular dynamics and Monte Carlo simulations, respectively. In the present study we use equilibrium and nonequilibrium MD simulations to calculate the structural, dynamical, and thermal proprties of methane and CO2 hydrate, as well as of Xe hydrate which has also been the subject of several experimental5,24,33-39 and theoretical38,40-43 studies. All three of these hydrates have the type I host lattice structure shown in Figure 1, consisting of two small (dodecahedral) 512 and six large (tetrakaidecahedral) 51262 water cages in the unit cell.1 †

Part of the “Barbara J. Garrison Festschrift”. * To whom correspondence should be addressed. E-mail: [email protected]. ‡ Present address: Department of Chemical Engineering, University of Michigan, Ann Arbor, MI 48109.

Figure 1. Structures of (a) tetrakaidecahedral (51262), (b) dodecahedral (512) cages, and (c) a unitcell of type I clathrate. Water molecules and guest species are represented by red and gray balls, respectively. Hydrogen bonds between water molecules are represented by red sticks.

The present study will thus provide information on the dynamical and thermal properties of CO2 hydrate and will also permit a systematic comparison of the properties of methane, CO2, and Xe hydrate. This comparison is particularly interesting as the methane molecule lacks a permanent dipole or quadruple moment, the carbon dioxide molecule is linear and has a large quadruple moment, and Xe, while lacking a permament moment, is much more polarizable than methane or CO2 and thus would be expected to have stronger polarization and dispersion interactions with the water molecules. For comparison, simulations have also been carried out on the hypothetical guest-free hydrate with the type I lattice structure. Computational Methods Molecular Models. In the simulations, nonpolarizable, rigid monomer models were employed for the H2O, CH4, and CO2 molecules. Although inclusion of explicit polarization can be important for some properties of hydrates,17-19,44 it significantly increases the computational time (compared to two-body potentials), and previous studies have show that simulations with two-body force fields do capture the qualitative trends17-19 which is the emphasis of the present study. Moreover, the use of rigid molecular models further reduces the computational cost as it permits longer time steps in the simulations and is not expected

10.1021/jp9063406  2010 American Chemical Society Published on Web 11/11/2009

5556

J. Phys. Chem. C, Vol. 114, No. 12, 2010

Jiang and Jordan

TABLE 1: Parameters for the Model Potentials Used in this Work Lennard-Jones parameters models water

46

methane14 xenon48 carbon dioxide51 water-methanea water-xenon b water-carbon dioxide50

a

atoms

partial charges

σ (Å)

ε (kJ/mol)

O H C H Xe C O C-O Owater-Cmethane Owater-Xe Hwater-Xe Owater-Ccarbon_dioxide Owater-Ocarbon_dioxide Hwater-Ccarbon_dioxide Hwater-Ocarbon_dioxide

-0.8476 0.4238 -0.5600 0.1400 0 0.5888 -0.2944 -

3.16556 3.64000 3.90110 2.79180 3.00000 2.89590 3.39450 3.43933 3.07455 3.11560 3.02970 2.81420 2.19420

0.65017 1.36500 1.89196 0.23983 0.68724 0.40598 0.94206 1.80837 0.11732 0.29300 0.53210 0.29973 0.39233

Estimated using combination rules.55 b Obtained from a least-squares fit of the ab initio data.49

TABLE 2: Binding Energies (kJ/mol) of Methane, Xenon, and Carbon Dioxide in the 512 and 51262 Water Cages guest molecule

512 cage

51262 cage

methane xenon carbon dioxide

18.61 37.15 26.13

18.64 39.23 30.16

to introduce significant errors in the properties of interest in this study.17,18,45 The SPC/E model46 was employed for water. This model uses three point charges to represent the electrostatics and a single Lennard-Jones site, located on the oxygen atom, to account for repulsion and dispersion interactions. For methane, partial charges were employed on the carbon and hydrogen atoms, and a single Lennard-Jones site, located on the carbon atom, was employed. The partial charges and the Lennard-Jones parameters in the methane model were taken from ref 14, and the Lorentz-Berthelot combination rule47 was used to determine the parameters for the Lennard-Jones interactions between the water and methane molecules. The Xe-Xe interactions were described by a Lennard-Jones potential using parameters from ref 48, and the water-xenon interactions were described by Lennard-Jones potentials between the xenon atom and the oxygen and hydrogen atoms of the water molecules, with the parameters being determined by a least-squares fit to the ab initio water-xenon potential of Wen and Ja¨ger.49 For CO2 partial charges and Lennard-Jones sites were employed on the carbon and oxygen atoms, with the charges being taken from ref 50, and the Lennard-Jones parameters for the CO2-CO2 and H2O-CO2 interactions taken from refs 50 and 51, respectively. The parameters for the model potentials used in this work are summarized in Table 1. As a useful reference point, we list in Table 2 the potential energies for binding the guest species in isolated 51262 and 512 water cages. The binding energies were calculated by optimizing the position and orientation of the guest molecule in rigid cages, with the cage geometries and the hydrogen network topology52 being taken from the optimized lattice structure of the hypothetical guest-free clathrate hydrate. As seen from Table 2, the binding energies of the guests in both the small and large water cages increase along the sequence CH4 < CO2 < Xe, in agreement with Myoshi.53 Moreover, while the binding energies of methane in the two cages are essentially identical, those of Xe and CO2 are, respectively, 2.1 and 4.0 kJ/mol greater in the larger cage. It should be noted that the 512 and 51262 cages used

in the above analysis are, in each case, but one of a large number of near zero dipole moment cages, differing in their H-bonding arrangements, that might be sampled in the simulation, and which differ in energy.52 However, we expect the relative stabilities of the different host species in the two types of cages to be relatively independent of the H-bonding arrangements. Simulation Details. The equilibrium MD simulations were carried out using a cubic supercell comprised of eight (2 × 2 × 2) unit cells. All water cages, except in the case of the hypothetical guest-free hydrate, were fully occupied by guest molecules. In the starting structures, the positions of the oxygen atoms of the water molecules were taken from the X-ray diffraction structure of ethylene oxide hydrate,54 the orientations of the water molecules were initialized in a random fashion according to the Bernal-Fowler rules,55 and the guest molecules were placed at the centers of the cages. The lattice constants and the lengths of the simulation boxes were determined from 100 ps equilibrium NPT (P ) 1 atm) simulations, employing the Berendsen weak-coupling thermostat and barostat,56 with the coupling constants being τT ) 0.1 ps and τP ) 0.5 ps, respectively. In the NPT simulations, the first 50 ps were used for equilibration, and the averages were carried out over simulation times of 50 ps. Radial distribution functions, autocorrelation functions, and power spectra were calculated from 20 ps production runs in the NVE ensemble, starting with structures generated in the NPT simulations. A time step of 1 fs was used in both the NPT and NVE simulations. These simulations were carried out using the DL_POLY57 code. The speed of sound is one of the key factors controlling heat transport in solids, and for each hydrate system considered, the velocity of the longitudinal sound wave propagating along the Z direction of the simulation box was calculated from NVT simulations, using the continuum-media model,58

Vp )



C11 F0

where F0 is the zero-stress mass density and C11 is the compressional elastic constant, defined as

C11 )

σzz εzz

Properties of Xenon, Methane, and Carbon Dioxide Hydrates with εzz and σzz being the strain and stress in the Z direction, respectively. (Since the unit cell is cubic, it is arbitrary whether one considers the X, Y, or Z directions.) A reference system with eight (2 × 2 × 2) unit cells was constructed using the zero-stress lattice constant, which was determined from a 1 ns NPT simulation with the pressure tensor set to zero and maintained by the Berendsen weak-coupling thermostat and barostat with the same coupling constants as used for the NPT simulations described above. The strained systems were prepared by rescaling the simulation box of the reference system in the Z direction by a factor of (1 - εzz). Constant-volume (NVT) simulations, with the temperature maintained using Berendsen weak-coupling thermostat with τT ) 0.1 ps, were then carried out on the strained systems, and the stress σzz was estimated using data accumulated over 1 ns. σzz was calculated with εzz) 0.0001, 0.0002, 0.0003, and 0.0004, and the C11 elastic constant was determined using linear regression. These simulations were carried out using the GROMACS59,60 code, again with a time step of 1 fs. The thermal conductivities of the hydrates were calculated using nonequilibrium molecular dynamics (NEMD) simulations in the NVT ensemble,61-63 with supercells containing a heat source and a heat sink. The thermal conductivities were calculated from the resulting temperature gradients and Fourier’s law.64 In order to eliminate finite size effects, at each temperature considered, the simulations were carried out with increasingly large simulation boxes, and the bulk thermal conductivity was obtained by extrapolating the system to infinite size in the direction of the temperature gradient. The simulation boxes were built with (2 × 2 × n) unit cells of hydrate, with n, the number of unit cells in the Z direction, being 2, 4, 6, and 8. A constant energy flux of about 10-10 w/Å2 was imposed along the Z axis by use of the velocity rescaling procedure.63 The temperature of the whole system was maintained using the Berendsen weakcoupling thermostat with a time constant of τT ) 2.0 ps. Typically, the first 500 ps of each NEMD simulation was used to establish a steady temperature profile and was not used in the averaging. This was followed by 2.4 ns production runs. A time step of 2 fs was used in the simulations. The NEMD simulations were performed using a modified version of the DL_POLY57 code. In both the equilibrium and nonequilibrium MD simulations, the intermolecular van der Waals interactions were cut off at a distance of 10.0 Å, with long-range corrections47 being applied for the energy in all simulations and also for the pressure in the NPT simulations. The long-range electrostatics were handled by the smoothed particle mesh Ewald (SPME) method.65,66 The MD simulations of the three hydrates and of the empty lattice were carried out for T ) 30, 50, 100, 200, and 260 K, all at P ) 1 atm. At this pressure the CO2, Xe, and methane hydrates are stable up to about T ) 217, 200, and 193 K, respectively.1 Hence, all three of the hydrates are metastable at T ) 260 K as is methane hydrate at T ) 200 K, and the empty lattice at all temperatures. Thus, it is important to note that even at the highest temperature considered there was no evidence of decomposition of any of these systems on the time scale of the simulations. Results and Discussion Structural Properties. A. Lattice Constants and Thermal Expansion. The lattice constants of the three hydrates and of the empty hydrate lattice are reported in Figure 2, together with the experimental lattice constants of the Xe,24 methane,10 and CO225 hydrates. The simulations predict that the lattice constants

J. Phys. Chem. C, Vol. 114, No. 12, 2010 5557

Figure 2. Lattice constants of the methane, Xe, and CO2 hydrates, as well as of the empty hydrate lattice as a function of temperature, from NPT simulations carried out at P ) 1 atm. Experimental lattice constants are plotted for comparison.

of the hydrates increase nearly linearly with increasing temperature. For temperatures less than about 150 K, the lattice constants of the guest-containing hydrates are smaller than those of the guest-free hydrate, with the calculated lattice constants being smallest for CO2 hydrate, next smallest for Xe hydrate, and largest for the empty lattice. At low temperatures, the relative ordering of the lattice constants calculated for the methane, Xe, and CO2 hydrates are in agreement with experimental data. The calculated lattice constants are ∼0.02 Å smaller than those obtained from experiment. This is primarily a consequence of the use of rigid-monomer models in the MD simulations.18 The smaller lattice constants and the generally larger rates of thermal expansion with increasing T of the guest-containing hydrates than for the guest-free hydrate can be understood as follows. At low temperatures, the guest species are localized near the centers of the cages,24 where the potential energies have their minima and attractive interactions dominate. As a result, the water molecules of the cages are pulled toward the cage centers, leading to reduced lattice constants. As the temperature increases, the guest species become more delocalized in the cages, with increasing probability for close approach to water molecules of the surrounding cages, repulsing the latter. At the temperature where the averaged repulsions and attractions between the encaged guest species and the water molecules in the surrounding cages balance one another, the lattice constant approaches that of the guest-free clathrate hydrate. At still higher temperatures, the repulsive effects dominate, causing the lattice to expand. The thermal expansion of the guest-containing hydrates thus can be regarded as resulting from a combination of both the intrinsic expansion of the empty lattice with increasing T and the perturbations due to the encaged guest species. As discussed above, the contribution of the guest species to the lattice expansion changes from negative to positive with increasing temperature, making the rate of expansion with increasing T larger for the guest-containing hydrates than for the guest-free hydrate. B. Radial Distribution Functions (RDFs). The spatial distributions of the guest species inside the water cages were investigated by calculating the RDFs between the center of mass of the encaged guest species and the oxygen atoms of the host water molecules. Panels a and b of Figure 3 report the T ) 30 K guest-water (GW) RDFs for the small (dodecahedral) and large (tetrakaidecahedral) cages, respectively. The locations of the maxima of the peaks in the RDFs are nearly identical for

5558

J. Phys. Chem. C, Vol. 114, No. 12, 2010

Figure 3. RDFs between the oxygen atoms of the host water molecules and the centers of mass of the guest species in the (a) small and (b) large water cages, respectively. The results are from simulations carried out at T ) 30 K and P ) 1 atm.

the three guests, although the peaks are appreciably broader for CO2 hydrate than for the Xe or methane hydrates, indicating that the CO2 molecules are more delocalized in the cages. The first peaks of the GW RDFs of carbon dioxide, xenon, and methane enclathrated in the small water cages have their first maxima around r ) 3.74 Å. On the other hand, in the large cages, the first peaks in the gGW RDFs of CO2 hydrate and methane hydrate occur near r ) 4.05 Å, whereas that of Xe hydrate occurs at r ) 3.93 Å, accompanied by a small shoulder near r ) 4.18 Å. The locations of the maxima of the first peaks in the T ) 30 K GW RDFs indicate that the encaged guest species are preferentially located near the centers of the water cages. The shoulder on the first peak of the RDF of xenon enclathrated in the large cages can be understood in terms of the noticeable asymmetry of tetrakaidecahedral cages. This feature is smeared out in the RDFs of the methane and CO2 hydrates due to the larger delocalization of the encaged methane and CO2 molecules. The first peaks in the gGW RDFs of the methane, xenon, and CO2 hydrates are considerably more asymmetrical for the tetrakaidecahedral cages than for the dodecahedral cages, which can also be attributed to the asymmetry of tetrakaidecahedral cages. The impact of the guest-host interactions on the lattice structures is examined in Figure 4a and b, where the T ) 30 K O-O and H-H RDFs of the host water molecules of the carbon dioxide, xenon, and methane hydrates, as well as of the guestfree hydrate are plotted. The calculated O-O RDFs of the methane, xenon, and CO2 hydrates are nearly identical to that of the guest-free hydrate. On the other hand, the first peak in

Jiang and Jordan

Figure 4. (a) Oxygen-oxygen and (b) hydrogen-hydrogen radial distribution functions of the host water lattice. The results are from simulations carried out at T ) 30 K and P ) 1 atm.

the calculated H-H RDF is appreciably lower in intensity for CO2 hydrate than for the Xe or methane hydrates, suggesting that the presence of guest CO2 molecules enhances the amplitude of the librational motion of the water molecules. This is consistent with the hydrogen atoms of the water molecules being attracted and released by the oxygen atoms of the CO2 molecules as the latter rotate inside the cages.24 In contrast, the encapsulation of methane or xenon does not significantly affect the librational motion of the water molecules. Dynamical Properties. A. Rattling Dynamics of the Encaged Guest Species. The rattling motions of the encaged Xe, CH4, and CO2 species were investigated using velocity autocorrelation functions (VACFs) and the associated Fourier spectra. The VACFs were calculated according to

VACF(τ) ) 〈VF(t)VF(t + τ)〉/〈VF(t)VF(t)〉 where F V (t) and F V (t + τ) are molecular center-of-mass (COM) velocities at time t and t + τ, respectively, and 〈...〉 represents the average over guest species and time t. The Fourier transforms of the VACFs give the corresponding density of states (DOS) associated with the molecular COM motions

DOS(ω) )

∫0∞ VACF(τ) cos(ωτ) dτ

Figure 5a reports the calculated VACFs of carbon dioxide, xenon, and methane, enclathrated in the small hydrate cages,

Properties of Xenon, Methane, and Carbon Dioxide Hydrates

J. Phys. Chem. C, Vol. 114, No. 12, 2010 5559

Figure 5. (a) Velocity autocorrelation functions and (b) associated spectra of the guest species in the small hydrate cages. The results are from simulations at T ) 30 K and P ) 1 atm.

Figure 6. (a) Velocity autocorrelation functions and (b) associated spectra of the guest species in the large hydrate cages. The results are from simulations at T ) 30 K and P ) 1 atm.

and Figure 5b reports the corresponding power spectra. The VACFs of all three guests in the small cages oscillate between positive and negative values with the oscillation being damped out with increasing time. The oscillations are damped much more strongly for CO2 hydrate than for the Xe or methane hydrates, indicating that the coupling between the rattling motions of the encaged guest molecules and the vibrational motions of the host lattice is strongest for CO2 hydrate. The Fourier spectrum for methane rattling in the small cages has a well-defined peak near 83 cm-1 and a long tail extending to lower frequencies, while that of xenon hydrate has a pronounced peak at 35 cm-1, a minor peak at 43 cm-1, shoulders at 54 and 23 cm-1, and a long tail extending to high frequencies. The Fourier spectrum of carbon dioxide in the small cages is characterized by broad structure spanning the range 10-120 cm-1, with weak bumps at 51 and 67 cm-1. While this structure is very different from that in the spectra calculated for methane hydrate and xenon hydrate, similar broad structure has been observed in the inelastic neutron scattering spectrum of N2 molecules in hydrate cages.39 The VACFs and the associated Fourier spectra of carbon dioxide, xenon, and methane in the large water cages are reported in Figure 6. As for the small cages, the VACFs for all three guest species display damped oscillatory behavior, with the oscillations being damped out much most rapidly for CO2, suggesting a stronger guest-host coupling for CO2 hydrate. The spectrum for the rattling of methane in the large cages displays peaks at 32 and 57 cm-1 and that for Xe rattling in the large cage displays peaks at 17 and 29 cm-1. In contrast, the spectrum of carbon dioxide in the large cage displays a single

broad peak with its maximum near 31 cm-1 and a tail extending to zero frequency. The low-frequency tail is absent in the spectra of Xe hydrate and methane hydrate. The doublets in the spectra of Xe and methane in the large cages are due to the asymmetry of the large cage. The peaks in the calculated rattling spectra of the Xe and CH4 hydrates are very close to those observed in inelastic neutron scattering measurements7,9,37-39,67 and in previous MD simulations of the spectra.13-18,38,40 B. Rotational Dynamics of the Encaged Carbon Dioxide Molecules. Inelastic neutron scattering (INS) and nuclear magnetic resonance (NMR) experiments by Tse et al.8 and Garg et al.,6 respectively, indicate that methane molecules in hydrate cages rotate nearly freely and rapidly even at temperatures as low as 4-5 K. Using the orientational order parameters calculated from molecular dynamics simulations, Tse et al.14 demonstrated that the enclathrated methane molecules have no preferred orientation relative to the surrounding water cages. Ikeda et al., based on Raman spectra68 and neutron powder diffraction24,25 measurements of CO2 hydrate, concluded that the CO2 molecules rotate anisotropically in the tetrakaidecahedral cages and that the surrounding water molecules are affected by the motion of the carbon dioxide molecules. They also suggested that, at low temperatures, the CO2 molecules rotate rapidly in the tetrakaidecahedral cages but that their rotational motion is suppressed in the dodecahedral cages.25 In this work, the rotational autocorrelation functions (RACF) and the associated Fourier spectra from MD simulations are used to to investigate the rotational motions of the enclathrated

5560

J. Phys. Chem. C, Vol. 114, No. 12, 2010

Figure 7. (a) Rotational autocorrelation functions and (b) the associated spectra of methane molecules in the hydrate small and large hydrate cages, from NVE simulations at T ) 30 and 260 K. Results are also included for the free CH4 molecule.

carbon dioxide and methane molecules. The rotational autocorrelation function was calculated as69

RACF(t) )

〈

N

∑

〉

1 Fu (0)u F (t) i N i)1 i

where F ui(t)(i ) 1 - N) is an instantaneous unit vector lying along a C-O bond of CO2 or a C-H bond of methane, N ) 2 and 4 for CO2 and methane melocules, respectively, and 〈...〉 denotes the ensemble average. The calculated RACFs of the enclathrated molecules are compared to those of the freely rotating molecules, with the latter being calculated using 50 ps NVE simulations of a box of 64 noninteracting molecules with the initial coordinates and velocities of the molecules and the sizes of the simulation boxes being taken from 100 ps NPT simulations in which the molecules were permitted to interact and which were carried out at T ) 30 and 260 K and P ) 1 atm. The RACFs of the freely rotating methane molecules and of the methane molecules enclathrated in the small and large cages are plotted in Figure 7a, and the corresponding spectra are reported in Figure 7b. The RACF for free-rotating methane molecules at T ) 260 K is characterized by a sharp declination from one to zero within the first 0.2 ps, followed by a fast rise to 0.37 in the subsequent 0.2 ps, with small oscillations about this value at still longer times. The RACF for free-rotating

Jiang and Jordan methane molecules at T ) 30 K is similar to that at T ) 260 K, with the initial declination and the long-time oscillations being somewhat slower. The Fourier spectrum of the freely rotating methane molecules exhibits a sharp central peak at zero frequency and well-defined, structured wings peaking around (23 and 86 cm-1 at T ) 30 and 260 K, respectively. At T ) 260 K, the RACFs of the enclathrated methane molecules in both small and large cages coincide almost exactly with that of freely rotating methane molecules in the first 0.3 ps; however, after 0.45 ps the RACFs of the encaged methane molecules decay exponentially with a time constant of about 0.9 ps. Correspondingly, the Fourier spectrum of enclathrated methane molecules at T ) 260 K resembles that of the freely rotating methane molecules but has a slightly broader central peak and much smoother side wings. At T ) 30 K, the RACFs of the enclathrated methane molecules are characterized by a decay in the first 0.7 ps, which nearly coincides with that of the freely rotating methane. After 0.7 ps the RACFs of the enclathrated methane molecules are characterized by a slower, exponential-like decay with a time constant of about 2 ps. In the corresponding Fourier spectrum of the enclathrated methane molecules the central peak is significantly broadened, and the wings collapse. The above analysis suggests that at high temperatures (∼260 K) the enclathrated methane molecules rotate almost freely inside the cages, in accordance with experimental observations. However, at low temperatures, the rotation of the enclathrated methane molecules is somewhat more hindered. In addition, the simulations show that the spectra of the methane molecules in the large and small cages are nearly indistinguishable. The RACFs of the enclathrated and free-rotating carbon dioxide molecules are plotted in Figure 8a, and the corresponding spectra are reported in Figure 8b. The T ) 260 K RACF of the freely rotating carbon dioxide molecules decreases rapidly from 1 to -0.42 in the first 1 ps, increases to zero during the next 1 ps, and then oscillates between (0.15 for the remainder of the simulation. The RACF of the freely rotating carbon dioxide molecules at T ) 30 K is similar to that at T ) 260 K but with the declination and oscillations being slower. The spectrum of the freely rotating carbon dioxide molecules is characterized by well-defined wings, peaking at (4 and 14 cm-1 at T ) 30 and 260 K, respectively. The RACFs and spectra of the enclathrated CO2 molecules deviate remarkably from those of the free rotating CO2 molecules. The T ) 260 K RACFs of the enclathrated CO2 molecules rapidly decay to about 0.2 in 1.2 ps and then slowly decay to zero with a 2.7 ps time constant. The spectra of the enclathrated CO2 molecules at T ) 260 display a single broad peak centered at zero frequency, with long tails extending to about (50 cm-1. The T ) 30 K RACFs of the enclathrated CO2 molecules decay extremely slowly. The spectra of enclathrated carbon dioxide molecules at T ) 30 are qualitatively similar to that for T ) 260 K, with the peak being significantly narrower. Interestingly, the RACFs and spectra for the enclathrated CO2 molecules are similar for the two size cages. These results indicate that the rotational motion of the CO2 molecules is strongly affected by their interactions with the water molecules of the surrounding cages at both high and low temperatures, with the magnitude of the perturbations being greater at low temperatures, where the enclathrated carbon dioxide molecules essentially undergoing wagging rather than rotational motion inside the cages. The prediction that the enclathrated carbon dioxide molecules rotate anisotropically is consistent with the conclusion of Ikeda et al. made on the basis

Properties of Xenon, Methane, and Carbon Dioxide Hydrates

J. Phys. Chem. C, Vol. 114, No. 12, 2010 5561

Figure 8. (a) Rotational autocorrelation functions and (b) the associated spectra of CO2 molecules in hydrate cages from NVE simulations at T ) 30 and 260 K. Results are also included for the free CO2 molecules.

Figure 9. DOS associated with the (a) translational and (b) librational motions of the water molecules in the methane, xenon, and CO2 hydrates as well as in the guest-free hydrate. The inset in (a) shows an enlarged part of the translational spectra in the low-frequency region. The simulation results are for T ) 30 K and P ) 1 atm.

of their inelastic neutron scattering measurements.24,25 However, their suggestion that the enclathrated carbon dioxide molecules rotate with significantly different speeds in large and small cages at low temperatures is not supported by our simulation results. The simulations described above are classical, and it is likely that inclusion of nuclear quantum effects could prove important for the rotational dynamics, especially for T ) 30 K. Nonetheless, the finding that CO2 is much more hinered in its rotation than is CH4 in the hydrate cages is likely to persist in a quantum treatment. C. Dynamics of the Host Lattice. Figure 9a reports the DOS associated with the translational motions of the host water molecules in the CO2, CH4, Xe, and guest-free hydrates at T ) 30 K. In general, the water translational DOSs of the different systems are quite similar, in each case having two broad features peaking near 70 and 300 cm-1. The most interesting region of the spectra is that below 40 cm-1 (region A), which is enlarged and reported in the inset of Figure 9a. In this region, the spectrum of the host lattice of CO2 hydrate decreases slowly with decreasing frequency, down to ω ) 0. This is in contrast to the spectra for the guest-free hydrate and for Xe hydrate and methane hydrate, which go to zero at finite (∼12-17 cm-1) frequencies. As noted above, the spectrum associated with the rattling of the guest CO2 molecules in the large cages (see Figure 6) also extends down to ω ) 0. These results are consistent with a strong coupling between the low frequency vibrations of the CO2 guest molecules and those of the water lattice.

The translational spectrum of the water lattice of xenon hydrate extends about 5 cm-1 lower in frequency than the corresponding spectra of the guest-free hydrate or of methane hydrate. The host translational spectrum of xenon hydrate has three well-defined, low-frequency peaks, centered at 17, 23, and 33 cm-1, very close in frequencies to the peaks (17, 29, and 35 cm-1) in the spectrum for the rattling motion of Xe in the hydrate cages. These values agree well with those obtained from previous INS experiments37-39,67 and computer simulations40,41 of Xe hydrate. The 23 and 33 cm-1 peaks in the host spectrum of xenon hydrate also fall close to the frequencies of the two lowest frequency features in the spectrum of the guest-free hydrate (which has peaks at 25, 33, 45, and 70 cm-1). The spectrum of the host lattice of methane hydrate has peaks at 26, 45, and 62 cm-1, while the three peaks for the rattling motion of the methane molecule in the hydrate cages occur at 33, 57, and 84 cm-1. There are also significant differences in the host translational spectra of the three hydrates in the range 250-400 cm-1 (region B). Specifically, for each hydrate species, the vibrational structure in this region is shifted to the blue of that of the empty lattice, with the shift being significantly greater for the Xe and CO2 hydrates than for the methane hydrate. In the librational region (540-1000 cm-1) the spectra of the three hydrate systems are similar to that of the empty hydrate lattice, with that for CO2 hydrate differing the most from that

5562

J. Phys. Chem. C, Vol. 114, No. 12, 2010

Jiang and Jordan

TABLE 3: Longitudinal Sound Speeds (km/s) of the Methane, Xenon, and Carbon Dioxide Hydrates and of the Guest-Free Hydrate

calculated measured

T (K)

guest-free hydrate

methane hydrate

xenon hydrate

carbon dioxide hydrate

30 150 260 2635 27770

4.7 ( 0.1 4.0 ( 0.1 3.8 ( 0.2 -

4.3 ( 0.1 4.2 ( 0.1 4.0 ( 0.1 3.4 ( 0.7 3.7 ( 0.1

3.3 ( 0.1 3.0 ( 0.1 2.8 ( 0.1 2.9 ( 0.1 -

3.5 ( 0.1 3.3 ( 0.1 3.1 ( 0.1 -

of the empty lattice. This in accord with our previous conclusion, based on the consideration of the radial distribution functions, that the CO2 guest molecules impact the librational motion of the water molecules more than do the Xe or CH4 guest species. Heat Transport Properties. The analysis of the structural, spectral, and dynamical properties of the carbon dioxide, xenon, and methane hydrates indicates that, in all three cases, the presence of the guest species perturbs the host lattice, with the perturbations being more pronounced for the carbon dioxide and xenon hydrates than for methane hydrate. In an earlier study19 we found that for methane hydrate guest-host coupling reduces the thermal conductivity by about 20% at T ) 30 K and that the impact of guest-host coupling on the thermal conductivity diminishes rapidly with increasing temperature. Considering that the guest-host coupling is even stronger in the carbon dioxide and xenon hydrates, we expect that the heat transport may be less effective in these crystals than in methane hydrate. This is verified by comparing the calculated sound speeds and thermal conductivities of CO2 and xenon hydrates with those of methane hydrate and the hypothetical clathrate hydrate with empty water cages. A. Sound Speed. The sound speeds of carbon dioxide, xenon, and methane hydrates and of the guest-free hydrate were calculated for T ) 30, 150, and 260 K using the approach described in the Simulation Details section. The results of the calculations are reported in Table 3. The measured sound speeds of methane5,70 and xenon5 hydrate have been included for comparison. The calculations predict that, as expected, sound speeds of all four systems decrease with increasing temperature. The sound speed of the guest-free hydrate is calculated to be 4.7, 4.0, and 3.9 km/s at T ) 30, 150, and 260 K, respectively. The sound speed of methane hydrate is calculated to be ∼9% smaller than that of the guest-free hydrate at T ) 30 K but nearly identical to that of the guest-free hydrate at T ) 260 K. The calculated sound speed of methane hydrate at T ) 260 K is 4.0 ( 0.1 km/s, which agrees reasonably well with the experimental values of 3.4 ( 0.4 at T ) 263 and 3.7 ( 0.1 km/s 277 K reported by Kiefte et al.5 and Waite et al.,70 respectively. The calculated sound speeds of carbon dioxide and xenon hydrates are significantly (>20%) lower than those of the guest-free and methane hydrates over the temperature range considered, with the sound speeds of carbon dioxide hydrate being slightly larger than those of xenon hydrate. The calculated sound speed of xenon hydrate at T ) 260 K is 2.8 km/s, which is in excellent agreement of the value (2.9 km/s for T ) 266 K) measured by Kiefte et al.5 The theoretical prediction of a lower sound speed of xenon hydrate than for methane hydrate is thus in accord with experiment. B. Thermal ConductiWities. Figure 10 summarizes our calculated thermal conductivities for the CO2, xenon, and methane hydrates, as well as of the guest free hydrate. For comparison, experimental thermal conductivities36 of Xe hydrate and methane hydrate11 are also reported. As noted in ref 19,

Figure 10. Calculated thermal conductivities of the guest-free, methane, xenon, and CO2 hydrates as a function of temperature.

for temperatures greater than 50 K, the calculated thermal conductivities of methane hydrate and the empty hydrate lattice are quite similar. Over the temperature range considered, the calculated thermal conductivities of Xe hydrate and CO2 hydrate are about 15% smaller than that of methane hydrate, which is consistent with the stronger perturbation of the lattice in the Xe and CO2 cases. In going from T ) 50 to 30 K, the thermal conductivities of the Xe, CO2, and methane hydrate decrease while that of the empty lattice increases. The different behavior of the thermal conductivities of the hydrates and the empty lattice at low temperatures is a result of guest-host coupling.19,71 The present results for Xe hydrate differ significantly from those of Inoue et al. whose NEMD simulations at T ) 100 K with the TIP4P force field gave a 12-fold lower thermal conductivity for Xe hydrate than for the empty hydrate lattice.41 The very different behavior predicted in the two sets of MD simulations appears to be the result of incomplete convergence of the calculations of Inoue et al., as has been discussed in ref 19. The calculated thermal conductivities are appreciably larger than the measured values obtained by Krivchikov et al.36 As shown in our earlier investigation,19 this is largely a consequence of the neglect of polarization in the force field used for the simulations. As seen from Figure 10, there is structure in the measured thermal conductivity vs T curves of both Xe and methane hydrate, which is not recovered in the calculations. For example, the measured thermal conductivity of Xe hydrate rapidly decreases from T = 60 to 100 K and then gradually rises, whereas the calculated thermal conductivity does not experience a rapid fall off between T ) 60 and 100 K and decrease weakly at still at higher temperatures, in contrast with experiment. A similar disagreement between the calculated and measured thermal conductivities exists for methane hydrate. Moreover, for T < 70 K, the experiments give higher values of the thermal conductivity for Xe hydrate than for CH4 hydrate, in contrast to the simulations. These differences could be due to quantum effects which are absent in the simulations. We are unaware of any experimental measurements of the thermal conductivity of CO2 hydrate. It is interesting to note that the calculated values of the thermal conductivities for methane and CO2 hydrate differ by about 0.2 W/m · K. If this difference is real, it could have important implications on the proposed process for utilizing CO2 as a means to stimulate methane recovery from hydrate.3,4

Properties of Xenon, Methane, and Carbon Dioxide Hydrates Conclusions In this work we made a systematic comparative studies of structural, dynamical, spectroscopic, and transport properties of the CO2, methane, and Xe hydrates using classical molecular dynamics simulations. The simulations indicate that of the three guest species, the CO2 molecules are the most delocalized in the water cages. Moreover, their rattling and rotational motions are significantly more hindered than for the CH4 or Xe guests, with the DOS spectra for the rattling motion of the CO2 molecules being spread over a large (∼100 cm-1) frequency range. The host lattice of CO2 hydrate expands more with increasing temperature than do the lattices of xenon and methane hydrate, and the dynamics of the water molecules is more perturbed in CO2 hydrate than in the xenon and methane hydrates. Specifically, both the translational and librational motions of the water molecules are significantly impacted by the encaged CO2 molecules. The impact on the water librations is attributed to the hydrogen atoms of the water molecules being temporarily attracted and released by the oxygen atoms of the carbon dioxide molecules as the latter rotate inside the cages. Interestingly, CO2 hydrate and Xe hydrate are calculated to have similar speeds of sound and thermal conductivity values, both of which are lower than the corresponding values of methane hydrate. Although guest-host coupling in the CO2, xenon, and methane hydrates significantly reduces the thermal conductivities at low (T e 50 K) temperatures, such couplings are predicted to be relatively unimportant at higher temperatures. Finally, the molecular level calculations indicate that at temperatures relevant for commercial applications the thermal conductivities of methane and CO2 hydrates may differ by about 0.2 W/m · K. This cannot be validated, owing to the lack of experimental thermal conductivity for CO2. Such validation would be most timely considering the current interest in utilizing CO2 to displace methane from natural deposits of methane hydrate. Acknowledgment. This technical effort is part of the National Energy Technology Laboratory’s ongoing research in gas hydrates under Contract No. DE-AC26-04NT41817, subtask 41817.606.08.102. The calculations were carried out on computers in the University of Pittsburgh’s Center for Molecular and Materials Simulations. References and Notes (1) Sloan, E. D.; Koh, C. A. Clathrate hydrates of natural gases, 3rd ed.; CRC Press: Boca Raton, FL, 2008. (2) Buffett, B. A. Annu. ReV. Earth Planet. Sci. 2000, 28, 477. (3) Goel, N. J. Petrol. Sci. Eng. 2006, 51, 169. (4) White, M. D.; McGrail, B. P. Numerical Simulation of Methane Hydrate Production from Geologic Formations via Carbon Dioxide Injection; Offshore Technology Conference, 2008, Houston, TX. (5) Kiefte, H.; Clouter, M. J.; Gagnon, R. E. J. Phys. Chem. 1985, 89, 3103. (6) Garg, S. K.; Gough, S. R.; Davidson, D. W. J. Chem. Phys. 1975, 63, 1646. (7) Sears, V. F.; Powell, B. M.; Tse, J. S.; Ratcliffe, C. I.; Handa, Y. P. Physica B 1992, 180, 658. (8) Tse, J. S.; Ratcliffe, C. I.; Powell, B. M.; Sears, V. F.; Handa, Y. P. J. Phys. Chem. A 1997, 101, 4491. (9) Tse, J. S.; Shpakov, V. P.; Belosludov, V. R.; Trouw, F.; Handa, Y. P.; Press, W. Europhys. Lett. 2001, 54, 354. (10) Takeya, S.; Kida, M.; Minami, H.; Sakagami, H.; Hachikubo, A.; Takahashi, N.; Shoji, H.; Soloviev, V.; Wallmann, K.; Biebow, N.; Obzhirov, A.; Salomatin, A.; Poort, J. Chem. Eng. Sci. 2006, 61, 2670. (11) Krivchikov, A. I.; Gorodilov, B. Y.; Korolyuk, O. A.; Manzhelii, V. G.; Conrad, H.; Press, W. J. Low Temp. Phys. 2005, 139, 693. (12) Rosenbaum, E. J.; English, N. J.; Johnson, J. K.; Shaw, D. W.; Warzinski, R. P. J. Phys. Chem. B 2007, 111, 13194.

J. Phys. Chem. C, Vol. 114, No. 12, 2010 5563 (13) Tse, J. S.; Klein, M. L.; McDonald, I. R. J. Phys. Chem. 1983, 87, 4198. (14) Tse, J. S.; Klein, M. L.; McDonald, I. R. J. Chem. Phys. 1984, 81, 6146. (15) Itoh, H.; Kawamura, K. Ann. N.Y. Acad. Sci. 2000, 912, 693. (16) Greathouse, J. A.; Cygan, R. T.; Simmons, B. A. J. Phys. Chem. B 2006, 110, 6428. (17) English, N. J.; MacElroy, J. M. D. J. Comput. Chem. 2003, 24, 1569. (18) Jiang, H.; Jordan, K. D.; Taylor, C. E. J. Phys. Chem. B 2007, 111, 6486. (19) Jiang, H.; Myshakin, E. M.; Jordan, K. D.; Warzinski, R. P. J. Phys. Chem. B 2008, 112, 10207. (20) English, N. J.; Tse, J. S. Phys. ReV. Lett. 2009, 103, 015901. (21) Zele, S. R.; Lee, S. Y.; Holder, G. D. J. Phys. Chem. B 1999, 103, 10250. (22) Chialvo, A. A.; Houssa, M.; Cummings, P. T. J. Phys. Chem. B 2002, 106, 442. (23) Hester, K. C.; Huo, Z.; Ballard, A. L.; Koh, C. A.; Miller, K. T.; Sloan, E. D. J. Phys. Chem. B 2007, 111, 8830. (24) Ikeda, T.; Mae, S.; Yamamuro, O.; Matsuo, T.; Ikeda, S.; Ibberson, R. M. J. Phys. Chem. A 2000, 104, 10623. (25) Ikeda, T.; Yamamuro, O.; Matsuo, T.; Mori, K.; Torii, S.; Kamiyama, T.; Izumi, F.; Ikeda, S.; Mae, S. J. Phys. Chem. Solids 1999, 60, 1527. (26) Udachin, K. A.; Ratcliffe, C. I.; Ripmeester, J. A. J. Phys. Chem. B 2001, 105, 4200. (27) Mohammadi, A. H.; Anderson, R.; Tohidi, B. AlChE J. 2005, 51, 2825. (28) Fan, S. S.; Chen, G. J.; Ma, Q. L.; Guo, T. M. Chem. Eng. J. 2000, 78, 173. (29) Nakano, S.; Moritoki, M.; Ohgaki, K. J. Chem. Eng. Data 1998, 43, 807. (30) Ohgaki, K.; Makihara, Y.; Takano, K. J. Chem. Eng. Jpn. 1993, 26, 558. (31) Adisasmito, S.; Frank, R. J.; Sloan, E. D. J. Chem. Eng. Data 1991, 36, 68. (32) Ferdows, M.; Ota, M. Chem. Eng. Technol. 2005, 28, 168. (33) Schober, H.; Itoh, H.; Klapproth, A.; Chihaia, V.; Kuhs, W. Eur. Phys. J. E 2003, 12, 41. (34) Stueber, D.; Jameson, C. J. J. Chem. Phys. 2004, 120, 1560. (35) Jameson, C. J.; Stueber, D. J. Chem. Phys. 2004, 120, 10200. (36) Krivchikov, A. I.; Gorodilov, B. Y.; Korolyuk, O. A.; Manzhelii, V. G.; Romantsova, O. O.; Conrad, H.; Press, W.; Tse, J. S.; Klug, D. D. Phys. ReV. B 2006, 73, 064203/1. (37) Gutt, C.; Baumert, J.; Press, W.; Tse, J. S.; Janssen, S. J. Chem. Phys. 2002, 116, 3795. (38) Itoh, H.; Chazallon, B.; Schober, H.; Kawamura, K.; Kuhs, W. F. Can. J. Phys. 2003, 81, 493. (39) Chazallon, B.; Itoh, H.; Koza, M.; Kuhs, W. F.; Schober, H. Phys. Chem. Chem. Phys. 2002, 4, 4809. (40) Tse, J. S.; Klein, M. L.; McDonald, I. R. J. Chem. Phys. 1983, 78, 2096. (41) Inoue, R.; Tanaka, H.; Nakanishi, K. J. Chem. Phys. 1996, 104, 9569. (42) Ikeda-Fukazawa, T.; Yamaguchi, Y.; Nagashima, K.; Kawamura, K. J. Chem. Phys. 2008, 129, 224506. (43) Tanaka, H. Chem. Phys. Lett. 1993, 202, 345. (44) Myshakin, E. M.; Jiang, H.; Warzinski, R. P.; Jordan, K. D. J. Phys. Chem. A 2009, 113, 1913. (45) Castillo-Borja, F.; Vazquez-Roman, R.; Bravo-Sanchez, U. Mol. Simul. 2008, 34, 661. (46) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. (47) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: New York, 1989. (48) Vrabec, J.; Stoll, J.; Hasse, H. J. Phys. Chem. B 2001, 105, 12126. (49) Wen, Q.; Ja¨ger, W. J. Phys. Chem. A 2006, 110, 7560. (50) Zhang, Z. G.; Duan, Z. H. J. Chem. Phys. 2005, 122. (51) Duan, Z. H.; Zhang, Z. G. Geochim. Cosmochim. Acta 2006, 70, 2311. (52) Yoo, S.; Kirov, M. V.; Xantheas, S. S. J. Am. Chem. Soc. 2009, 131, 7564. (53) Miyoshi, T.; Imai, M.; Ohmura, R.; Yasuoka, K. J. Chem. Phys. 2007, 126. (54) McMullan, R. K.; Jeffrey, G. A. J. Chem. Phys. 1965, 42, 2725. (55) Bernal, J. D.; Fowler, R. H. J. Chem. Phys. 1933, 1, 515. (56) Berendsen, H. J. C.; Postma, J. P. M.; Van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684. (57) Smith, W.; Forester, T. R. J. Mol. Graphics 1996, 14, 136–141. (58) Ibach, H.; Lu¨th, H. Solid-state physics: an introduction to principles of materials science, 3rd ed.; Springer: Berlin; New York, 2003.

5564

J. Phys. Chem. C, Vol. 114, No. 12, 2010

(59) Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R. Comput. Phys. Commun. 1995, 91, 43. (60) Lindahl, E.; Hess, B.; van der Spoel, D. J. Mol. Modeling 2001, 7, 306. (61) Muller-Plathe, F. J. Chem. Phys. 1997, 106, 6082. (62) Oligschleger, C.; Schon, J. C. Phys. ReV. B 1999, 59, 4125. (63) Jund, P.; Jullien, R. Phys. ReV. B 1999, 59, 13707. (64) Incropera, F. P. Fundamentals of heat and mass transfer, 6th ed.; John Wiley: Hoboken, NJ, 2007. (65) Ewald, P. P. Ann. Phys. 1921, 64, 253. (66) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. J. Chem. Phys. 1995, 103, 8577.

Jiang and Jordan (67) Schober, H.; Itoh, H.; Klapproth, A.; Chihaia, V.; Kuhs, W. F. Eur. Phys. J. E 2003, 12, 41. (68) Ikeda, T.; Mae, S.; Uchida, T. J. Chem. Phys. 1998, 108, 1352. (69) Danten, Y.; Tassaing, T.; Besnard, M. J. Phys. Chem. A 2000, 104, 9415. (70) Waite, W. F.; Helgerud, M. B.; Nur, A.; Pinkston, J. C.; Stern, L. A.; Kirby, S. H.; Durham, W. B. Ann. N.Y. Acad. Sci. 2000, 912, 1003. (71) Tse, J. S.; White, M. A. J. Phys. Chem. 1988, 92, 5006.

JP9063406