10608
J. Phys. Chem. B 2008, 112, 10608–10618
Molecular Dynamics Study of Methane Hydrate Formation at a Water/Methane Interface Junfang Zhang,† R. W. Hawtin,‡ Ye Yang,§ Edson Nakagava,| M. Rivero,| S. K. Choi,† and P. M. Rodger*,‡ CSIRO Petroleum, PriVate Bag 10, South Clayton, Victoria, 3169, Australia; Department of Chemistry, UniVersity of Warwick, CoVentry CV4 7AL, U.K., School of Biological Sciences, Nanyang Technical UniVersity, 60 Nanyang DriVe, Singapore 637551, and CSIRO Petroleum, P.O. Box 1130, Bentley, Western Australia, 6102, Australia ReceiVed: August 28, 2007; ReVised Manuscript ReceiVed: April 29, 2008
We present molecular dynamics simulation results of a liquid water/methane interface, with and without an oligomer of poly(methylaminoethylmethacrylate), PMAEMA. PMAEMA is an active component of a commercial low dosage hydrate inhibitor (LDHI). Simulations were performed in the constant NPT ensemble at temperatures of 220, 235, 240, 245, and 250 K and a pressure of 300 bar. The simulations show the onset of methane hydrate growth within 30 ns for temperatures below 245 K in the methane/water systems; at 240 K there is an induction period of ca. 20 ns, but at lower temperatures growth commences immediately. The simulations were analyzed to calculate hydrate content, the propensity for hydrogen bond formation, and how these were affected by both temperature and the presence of the LDHI. As expected, both the hydrogen bond number and hydrate content decreased with increasing temperature, though little difference was observed between the lowest two temperatures considered. In the presence of PMAEMA, the temperature below which sustained hydrate growth occurred was observed to decrease. Some of the implications for the role of PMAEMA in LDHIs are discussed. 1. Introduction Methane hydrate is a type of “clathrate”. Clathrates are complexes formed when molecules of one type are trapped inside a “host lattice” formed by another type of molecule.1 In the case of methane hydrate, the host lattice is a hydrogenbonded network of water molecules with methane molecules trapped in small polyhedral cages formed by the water network.2,3 In fact, many different compounds can form clathrate hydrates. For practical applications, the most important ones are components of natural gas, but in general they can range from small noble gas atoms to organics comparable in size with adamantine. The different components typically give rise to one of three types of gas hydrate that differ in the crystal structure of the host lattice and hence in the number and size of the cages found in the lattice. The type of the gas hydrate formed depends primarily on the size of the guest molecule. All these hydrate structures consist of at least two distinct types of water cage, and the crystal structure can be understood in terms of how these cages pack to fill 3D space. Methane hydrate is found in vast quantities in natural environments, such as permanently or seasonally frozen ground in polar regions, or marine sediment exposed to critical temperature/pressure environments. They are of an immediate and practical concern, due to the hazards they pose to oil and gas drilling and production operations.4 The formation of gas hydrate under certain temperature and pressure conditions can cause costly and hazardous blockages in natural gas pipelines.5 * To whom correspondence should be addressed. E-mail: p.m.rodger@ warwick.ac.uk. † CSIRO Petroleum, South Clayton, Victoria, Australia. ‡ University of Warwick. § Nanyang Technical University, Singapore. | CSIRO Petroleum, Bentley, Western Australia, Australia.
In fact, many components of natural gas in addition to methane will form gas hydrates at low temperatures and elevated pressures. Much experimental and computational research over the last ten years has centered on identifying “low dosage hydrate inhibitors” (LDHIs); these are additives that are active at low concentrations (typically about 1% by weight of water) and have the effect of delaying nucleation, slowing growth, or preventing aggregation of multiple crystals, rather than shifting the thermodynamic phase boundaries.6–13 Research to date has identified a number of compounds that can act as LDHIs, but no definitive underlying mechanism for their effect has yet been established. It is important to identify such an underlying mechanism, as this would considerably increase the efficiency with which new, more powerful, inhibitors can be developed. A number of molecular computations on the phase equilibrium for Structure I gas hydrates have been performed.14–19 Rodger has presented the results of a molecular dynamics computer simulation study of the factors influencing the stability and thermodynamic properties of gas hydrates.20 Much of the previous research on natural gas hydrate has focused on the equilibrium properties of gas hydrates, such as pressure measurements as a function of temperature, gas consumption versus temperature or structural properties.7,21–27 A number of studies have also been performed that are relevant to the mechanisms by which hydrate growth commences. Chau and Mancera performed Monte Carlo simulations to investigate the effect of pressure on the structural properties of the methane hydration shell in liquid water,28 while Mancera and Buckingham used molecular dynamics (MD) simulations to study hydrophobic hydration and its influence on ethane aggregation in water at different temperatures.29 Several different models for hydrate nucleation have been proposed,20,21,30 but all agree that the nucleation
10.1021/jp076904p CCC: $40.75 2008 American Chemical Society Published on Web 07/31/2008
MD of Methane Hydrate
J. Phys. Chem. B, Vol. 112, No. 34, 2008 10609
Figure 1. A configuration showing the design of the simulation cell. A film of water with dissolved methane is surrounded by a methane fluid, and the entire cell replicated periodically in three dimensions. Color scheme: O, red; H, white; C, cyan; N, dark blue; methane molecules are depicted by the small cyan spheres.
Figure 2. Hydrate clusters present in the initial configuration: (a) shows all clusters within the water film, with each distinct cluster depicted in a different color; (b) shows the largest cluster (blue) as a contiguous cluster rather than split by periodic boundaries as in panel a. A group of at least three near-complete dodecahedra cages can be seen at the center of the cluster in panel b with more partial cages present at the bottom of the cluster.
is a dynamic process. It is therefore desirable to study underlying dynamics at a molecular level. MD has proved useful in studying the dynamical processes accompanying nucleation in a number of other systems.31–33 More recently, it has also been used to simulate nucleation in methane hydrate directly,34 which has opened up the possibility of
Figure 3. Configuration of the simulation with PMAEMA at 220 K after 20 ns. The hydrate water molecules (identified using the local phase assignment) have been colored green and the hydrogen bonding network between them depicted as green cylinders; all other colors are as for Figure 1. Panel b shows the same configuration as in panel a but with the methane and nonhydrate water molecules removed for clarity.
10610 J. Phys. Chem. B, Vol. 112, No. 34, 2008 studying the influence of LDHIs on methane hydrate nucleation with greater detail than has previously been possible.35 MD simulations of known inhibitors in liquid water, and adsorbing onto a fixed hydrate surface, have been reported in the past.36–41 Configuration-biased Monte Carlo techniques have also been used to grow polymer on the hydrate surface.42 Water-methane-hydrate mixtures in their full dynamical complexity12,34,43 and hydrate nucleation and growth44 have been successfully modeled using MD simulations. Recently, Hawtin and Rodger45 have simulated the kinetic inhibitor poly(dimethylaminoethylmethacrylate), PDMAEMA, in a range of conformations and tacticities in a system of freely mobile and evolving hydrate and water. However, to the best of our knowledge, there are no comparable simulation studies showing the influence of temperature on the initial stages of hydrate nucleation and growth, nor on how this might be
Zhang et al. modified by the presence of additives. The purpose of this paper is to present just such a temperature study. In this paper, we present the results of a molecular dynamics study of liquid-water/methane mixtures in the presence of poly(methylaminoethylmethacrylate), PMAEMA. PMAEMA was chosen as a useful comparison with the earlier PDMAEMA simulations: PMAEMA has a polar hydrogen on the pendant amine groups and so can act as an H-bond donor to the water; PDAEMA has no such polar hydrogen atoms. Five different temperatures have been considered in long-time MD simulations at constant pressure and temperature (NPT MD). The remainder of this paper is structured as follows: in Section 2 we explain simulation details; in Section 3 the results and discussions are presented; and finally our analysis is summarized and concluded in Section 4.
Figure 4. Comparison of the hydrate content for the systems at temperatures of 220, 235, 240, 245, and 250 K. (a) with PMAEMA; (b) without PMAEMA; (c) with and without PMAEMA at 240 K.
MD of Methane Hydrate
Figure 5. Evolution of the hydrate content for the repeat simulations at temperature of 240 K: (a) with PMAEMA; (b) without PMAEMA.
2. Simulation Details In this study, the MD simulations were performed using GROMACS.46,47 The systems we have studied consist of CH4, H2O, and PMAEMA. Figure 1 presents an image of the simulation cell and thereby provides a schematic representation of the overall structure of the system. A film of water containing some dissolved methane (ca. 60 Å thick) was surrounded by methane fluid at 300 bar and the specified temperature. For the simulations with inhibitor, an oligomer of PMAEMA was introduced into the methane fluid phase and allowed to migrate to the aqueous phase by natural diffusion. The decision to simulate a methane/water interface, rather than hydrate formation in bulk water, was made to ensure that the mass transport limitations that control experimental hydrate formation with hydrophobic guests were included explicitly in the simulations. Simulations were performed for two systems: one with PMAEMA added and the other without the PMAEMA. A configuration from a pre-existing simulation of methane hydrate nucleation and growth44,45 was used as the initial point for simulations; this configuration was chosen from a period about 6 ns into the original simulation, during a period where a substantive cluster of methane hydrate (containing about 120 water molecules and associated methane; see Figure 2) had formed, and where there was little evidence of methane aggregation within the water layer. In detail, this configuration contained 1656 water and 1089 methane molecules in an orthorhombic cell of initial dimensions 3.33 × 2.96 × 15.15 nm3. The water was contained in an aqueous film ca. 6 nm thick. Methane was dispersed through this aqueous film at a mole fraction of 0.10; this mole fraction is significantly smaller than that of the hydrate phase (0.15), but much higher than concentration of methane found in bulk water. Some enrichment
J. Phys. Chem. B, Vol. 112, No. 34, 2008 10611 over the bulk concentration is to be expected in a thin water film in contact with methane at 300 bar, and this is likely to be the reason why methane hydrate is observed to nucleate preferentially at the water/methane interface, but the high concentration of aqueous methane must be expected to provide a large thermodynamic driving force for the formation of the hydrate phase in our simulations. We have also performed simulations at lower supersaturations (mole fractions of ca. 0.0535,48) and, apart from exhibiting longer induction times, these showed little difference in the nature or rate of the hydrate growth observed. The region surrounding this water film contained methane fluid (note that the critical temperature for methane is 190.6 K49). In the absence of the PMAEMA, this configuration was used without modification in isothermal-isobaric (NPT) simulations with the barostat being applied only in the z direction (i.e., normal to the water/methane interface). For the inhibited system, an octamer of PMAEMA was also included. The protocol for inserting the octamer into the methane hydrate interface was the same as that used in previous studies.44,45 Simulations were performed using an initial linear backbone conformation for the PMAEMA (180° for all backbone dihedral angles) and the conformation prepared in syndiotactic form. In previous research,45 tacticity was not found to influence the activity of PDMAEMA and so it was not considered as a factor in the present study. One octamer of PMAEMA was inserted into the methane gas phase in the simulation box, about 15 Å above the water surface. Methane molecules that overlapped the inhibitor were removed and the number of methane molecules in the hydrocarbon phase adjusted so that each system contained 1042 methane molecules. The insertion was always made into the methane fluid so as to avoid any disruption of the growing hydrate crystal due purely to the insertion process. To enhance comparability, all simulations at the different temperatures used the same initial configuration for the water molecules and for the aqueous methane molecules. Following previous studies,34,44,45 we used the simple point charge (SPC) water model50 and a united atom model for methane. The CHARMm force field (united atom + polar hydrogens) was used to model the inhibitor. Van der Waals force calculations were truncated at a distance of 1.2 nm. For calculating of long-range electrostatic forces, the Particle-Mesh Ewald (PME) method proposed by Darden51,52 was used. All bond lengths were constrained using the Shake algorithm.53 The temperature was fixed by using Berendsen thermostat54 at 220, 235, 240, 245, and 250 K, and the Berendsen pressure coupling algorithm was used to keep the pressure constant at 300 bar; this was applied only in the z direction (i.e., normal to the interface). Periodic boundary conditions were applied in three directions. The equations of motion were integrated with a 2 fs time step, though some simulations were repeated with a 1 fs time step and found to yield no significant differences from the larger time step. Each MD run was done in two stages. The first stage consisted of a 5 ps simulation with the steepest-descent method to perform energy minimization in order to reduce the thermal noise in the structures and potential energies. The second step consisted of a 30 ns production run. 3. Results and Discussion 3.1. Hydrate Content and Nucleation. The understanding of nucleation processes is built very firmly on classical nucleation theory. This uses a two-state model of the crystal nucleus, in which a crystalline core is surrounded by an
10612 J. Phys. Chem. B, Vol. 112, No. 34, 2008
Zhang et al.
Figure 6. Comparison of the hydrogen bond numbers for the systems at temperature 220, 235, 240, 245, and 250 K. (a) with PMAEMA; (b) without PMAEMA; (c) with and without PMAEMA at 220 and 235 K; (d) with and without PMAEMA at 240, 245, and 250 K.
interface, to explain why there is a free energy barrier to crystal formation. However, it is important to realize that a truly crystalline core is not a prerequisite for nucleation; rather, nucleation has occurred when a “crystallite” is more likely to grow than to shrink. Indeed, the available experimental evidence suggests that nucleation of hydrate structure, or possibly a precursor phase, precedes the formation of the stable thermodynamic phase: NMR studies on Xe hydrate55 and Raman studies of methane hydrate56 both indicate that the relative occupancy of the different cavities changes continuously during initial growth, and only attains the equilibrium ratio some time after nucleation, while X-ray studies of CO2 hydrate57 show evidence of the (metastable) type II structure during nucleation. This behavior is not unique to clathrate hydrates, and there is now a considerable body of evidence to show that the initial deposition of CaCO3 in nature occurs through the formation of amorphous nanoparticles, which subsequently aggregate and transform into the crystalline polymorphs.58–61 In view of this, to quantify the formation and growth of hydrate in our simulations we have developed a method (the “local phase assignment”) that identifies hydratelike structure in the molecular environment of each water molecule.40,62 The method is based on a set of three local order parameters designed to differentiate between hydrogen bond networks found in ice, hydrate, and liquid water structures; the analysis does not
differentiate between structure I and II hydrates, and so is appropriate for analyzing the initial formation in which elements of both structures may be present.57 Molecules are assigned a local phase only when this set of order parameters is consistent with the range of order observed in the respective bulk phases; the method correctly assigns 95% of the water molecules in simulations of bulk hydrate and of ice, and identifies a background level of 6% icelike and 6% hydratelike water molecules in liquid (SPC) water at 270 K but without any clustering of these solidlike molecules. Clusters of solidlike water molecules in our interfacial simulation can then be identified from the hydrogen bond network, defined using an O-O distance of 3.5Å. The outworking of this analysis is illustrated in Figure 3, which shows a representative example of the hydrate content in a snapshot of the configuration for the inhibited system at 220 K taken after 20 ns. The time evolution of the hydrate content for both inhibited and uninhibited systems is shown in Figure 4. The data for these curves were obtained by analyzing one configuration every 20 ps throughout each simulation. In general terms the data agrees well with that of Moon et al.,34,44 although there are small quantitative differences in that the onset of hydrate growth occurs at slightly lower temperatures in the present study (240 K instead of 250 K). This difference could be simply a manifestation of the inherently stochastic nature of crystal
MD of Methane Hydrate
J. Phys. Chem. B, Vol. 112, No. 34, 2008 10613
Figure 7. Comparison of the total internal energy for the systems at temperatures 220, 235, 240, 245, and 250 K. (a) with PMAEMA; (b) without PMAEMA.
nucleation, or might correlate with the different barostat methods used (isotropic in the work of Moon et al. compared with normal to the interfacial plane in the current work); however, given the consistent temperature trends reported below, the origin of these differences from Moon et al. are unlikely to be significant in the context of the present study. Three types of behavior are evident in Figure 4. At low temperatures (220 and 235 K), both the inhibited and uninhibited systems show immediate and sustained growth of the hydrate phase with no induction time evident. We conclude that the initial 120 water molecule hydrate cluster is larger than the critical cluster size under these subcoolings and supersaturation. At the higher temperatures (245 and 250 K), the hydrate content decreases significantly both with and without the inhibitor, reaching an apparent plateau of about 11 ( 2% prior to the end of each simulation. The final hydrate content is slightly lower in the presence of PMAEMA than it is without, but the difference is comparable with the observed fluctuations in the plateau regions The average plateau values are 181 (10.9%) and 187 (11.3%) water molecules for the inhibited and uninhibited systems at 245 K; 164 (9.9%) and 166 (10.0%) at 250 K. The third type of behavior is seen at the intermediate temperature (240 K) and is strikingly different. At this state point, the uninhibited system exhibits an induction time of 18 ns before sustained hydrate growth sets in; during the induction
period the hydrate content remains stable at about 15% of the available water. In contrast, the inhibited system shows a slow but steady decrease in the hydrate content throughout the duration of the 30 ns simulation. The difference between these two systems is seen clearly in Figure 4c. We note that crystal nucleation is a stochastic event, and so must consider the possibility that nucleation in the inhibited system is simply delayed by chance. Repeats of the simulations at 240 K were therefore performed, both with and without PMAEMA. These repeat simulations started from the same set of atomic positions as the first, but differed in that a different set of atomic velocities was allocated (again, chosen at random from a Boltzmann distribution at 240 K); the repeat simulations used GROMACS (version 3.3.1) and were run on the Altix 3700 BX2 (Itanium2 1.6 GHz). A total of two 30 ns simulations were performed for the uninhibited system and three in the presence of PMAEMA, and the results for the evolution of the hydrate content are depicted in Figure 5. In the absence of PMAEMA an induction time of about 20 ns followed by sustained growth of hydrate is seen in both systems (see Figure 5b), while a steady loss of hydrate structure is seen in all three simulations with PMAEMA present (see Figure 5a). We conclude that PMAEMA does actively inhibit the formation of the hydrate critical nucleus, but that the efficacy of this depends strongly on the degree of
10614 J. Phys. Chem. B, Vol. 112, No. 34, 2008
Zhang et al.
Figure 8. Radial distribution functions (RDFs) at 220, 240, and 250 K: (a) O-OH; (b) O-OL; (c) N-OH; (d) O-CH4; (e) N-CH4; (f) CH4-OH. In these labels, O is that carbonyl atom in PMAEMA, N is the amine nitrogen atom in PMAEMA, OH is the oxygen atom from hydrate water molecules, and OL is the oxygen atom from liquid water molecules.
subcooling. Results at 240 K presented in the rest of this paper are averaged over these repeat simulations. 3.2. Hydrogen Bonds. The hydrogen bonds (H-bonds) between all possible donors and acceptors were analyzed. Hydrogen bonds were defined to exist when both a hydrogendonor-acceptor angle (
