Methane hydrate: shifting the coexistence temperature to higher temperatures with an external electric field

This is an Accepted Manuscript of an article published by Taylor & Francis in “Molecular simulation” on 25th May 2016, available online: http://www.tandfonline.com/doi/abs/10.1080/08927022.2016.1139704 Methane hydrate: shifting the coexistence temperature to higher temperatures with an external electric field D.P. Luisa, J. López-Lemusb, M. Ll. Maspochc, E.A. Franco-Urquizaa and H. SaintMartind aCONACYT Research Fellow-Centro de Ingeniería y Desarrollo Industrial, Querétaro, México; bFacultad de Ciencias, Universidad Autónoma del Estado de México, Toluca, México; cCentre Català del Plàstic, Universitat Politècnica de Catalunya, Terrassa, Spain; dInstituto de Ciencias Físicas, Universidad Nacional Autónoma de México, Morelos, México ABSTRACT In the present work, we used molecular dynamic simulations of the equilibrium NPT ensemble to examine the effect of an external electric field on the three-phase coexistence temperature of methane gas, liquid water and methane hydrate. For these simulations, we used the TIP4P/Ice rigid water model and a singlesite model for methane. The simulations were implemented at two pressures, 400 and 250 bar, over temperatures ranging from 285 to 320 K and from 280 to 315 K, respectively. The application of an external electric field in the range of 0.1–0.9 Vnm−1 caused the effect of the thermal vibrations of the water molecules to become attenuated. This resulted in a shift of the three-phase coexistence temperature to higher temperatures. Electric fields below this range did not cause a difference in the coexistence temperature, and electric fields above this range enhanced the thermal effect. The shift had a magnitude of 22.5 K on average. 1. Introduction The clathrate hydrates are non-stoichiometric inclusion compounds in which water molecules form cavities, or ‘cages’, that trap molecules of gas, called guest molecules. These ‘cages’ are held together by the hydrogen bonding of the water molecules [1]; the three main crystal structures in which such cages form have been labelled sI, sII and H.[1,2] Alarge amount of methane is presumably stored in hydrates at the bottom of the ocean; in some investigations, it has been estimated that the energy that could be extracted from this source amounts to twice that of all other fossil fuels combined.[3]Moreover, the combustion of methane is cleaner because it produces only water and carbon dioxide; however, both the latter and methane are greenhouse gases. Thus, apart from technologies for safe and efficient extraction, the exploitation of oceanic methane hydrates would require the development of technologies to expediently condense H2O and reinsert CO2 into the hydrate from which the CH4 was extracted. Gas hydrates can also form during gas/oil transport, in which they pose a major problem because they can plug transportation pipelines.[4] In fact, flow assurance for the oil and gas industry is considered the most important application of gas hydrate research, which includes the study of the molecular pathways and the chemical and physical concepts underlying gas hydrate formation and decomposition.[5] Thus, significant efforts have been made in the past 5 years to better understand the mechanisms of natural gas hydrate crystallisation and dissociation, especially through the use of numerical simulations. Several studies have focused on the capabilities of various techniques and molecular models to describe the experimental conditions required for the coexistence of liquid water–methane hydrate,[6–13] whereas the current ability to perform simulations within microseconds [14,15] has allowed for the direct study of the nucleation and crystallisation processes,[16–25] the dissociation processes,[26–29] and even how these processes are affected by various factors, such as bubble formation [30] or the presence of NaCl in an aqueous solution.[31] It is expected that this knowledge will contribute to the development of better criteria for choosing among the various chemical inhibitors of hydrate formation that are available [32] for specific situations or even to modification of these inhibitors and the design of new ones. Physical methods can also be used to prevent hydrate formation or to destroy an already existing clog; a more practical alternative than heating or depressurising is to apply an external electric field.[33] This technique is routinely applied in the food industry to defrost or freeze food, for instance, freezing food under a static electric field; in this case, the primary objective is to improve the quality of frozen foods by decreasing the size of the ice crystals that form.[34,35] Moreover, electric fields of up to 50Vnm−1 can be experimentally produced by applying potentials of 1.0–5.0 kV onto tips with radii of 10–100 nm.[36] In simulations, static electric fields have been applied during the crystallisation of supercooled water, inducing dramatic electro-freezing [37,38]; in addition molecular dynamics (MD) simulations using the TIP4P [39] andST2 [40–43]water models have focused on the alignment of water molecule clusters that can be produced through the application of a constant field [44– 46] and the consequent orientations of the dipole moments. This has been observed to occur in simulations at T = 243 K, with electric fields in the 1– 5Vnm−1 range. Shevkunov and Vegiri [45] demonstrated the destruction of a normal solidor liquidlike cluster state using computer simulation of water clusters Figure 1. (Colour online) Snapshot of the initial configuration used for all simulations in this work. Figure 2. (Colour online) Graphs of potential energy as a function of time for determining the coexistence temperature at a fixed pressure; in (a), the temperature obtained was 287.5 K (+/−2.5 K) at a pressure of 250 bar, and in (b), the temperature obtained was 297.5 K (+/−2.5 K) at a pressure of 400 bar. (N = 40) in the presence of an external electric field in the 0.5–7.0Vnm −1 range. They showed that this process has the characteristics of a first-order phase transition. Aragones et al. [47] showed that an external electric field can modify the phase diagram of water; they developed Monte Carlo simulations with a specific field of 0.3Vnm−1 and observed a displacement of the phase boundaries as an effect of the electric field, e.g. the melting points of ice III and ice V were found to increase by approximately 15K. English and MacElroy [48] studied the effects of electric fields on the structures of methane hydrates via MD techniques, they used the rigid/polarisable TIP4P-FQ [49] water potential in simulations of systems consisting of a methane hydrate crystallite surrounded by a liquid phase under the influence of electromagnetic fields in the 5–7500GHz range at root mean square electric field intensities of up to 2Vnm−1. They reported that applying the field disrupted the hydrogen bonds that held the water cages together and ultimately led to the dissociation of the hydrate at intensities above 1Vnm−1; the most effective frequency range for crystal break-up was approximately 50– 100 GHz. Meanwhile, Luis et al. [8,50] tested various water models, such as SPC/E,[51] TIP4P [39] and TIP5P,[52] to study the dissociation of sImethane hydrate and found that in a system consisting of a methane hydrate crystallite without interfaces, dissociations occurred when an external electric field greater than 1.5Vnm−1 was applied at a pressure of 20 bar and a temperature of 248 K. For calculations of phase equilibrium conditions, Jensen et al. [7] used the TIP4P/Ice [53] model and a united-atom Lennard-Jones (LJ) potential for methane to determine the three-phase equilibrium conditions by using the Monte Carlo simulation method to calculate the chemical potentials of water and methane in the hydrate phase and in the liquid and gas phases, respectively. They found that the deviation in the coexistence temperature between the simulated and experimental data increased with pressure, and their calculations showed good agreement between the simulated and experimental three-phase coexistence curves. Recently, the methane hydrate three-phase equilibrium via a direct coexistence method [54–56] using MD simulations has been studied, for instance, Conde and Vega [9,10] used the TIP4P,[39] TIP4P/Ice [53] and TIP4P/2005 [57] water models showed that this method reproduced the experimentally determined three-phase coexistence temperature for the hydrate– liquid water–methane gas system. They showed that the TIP4P/Ice model yielded the best prediction. In addition, Tung et al. [58] used the same methodology but with a different initial configuration, and using the TIP4P-Ew [39] water potential and the OPLS-AA force field for methane, they found that factors such as the solubility of methane in water, the mass transport of methane via diffusion and the affinity of methane for incomplete water cages at the interface are three key factors that drive the growth rate of methane hydrate. More recently, Michalis et al. [59] used the direct phase coexistence method in a system that consisted of a solid hydrate slab, two liquid water slabs surrounding the hydrate slab and one gaseous methane slab. Using the TIP4P/Ice water model in combination with the OPLS-UA force field for methane, they found that the statistical analysis of multiple independent runs per pressure and long-term simulations were necessary to compensate for the inherent stochasticity of the method. The results followed the same trend as the experimental values and were highly similar to the experimental findings. In another study, Smirnov and Stegailov [12] used the TIP4P/Ice, [53] TIP4P/2005 [57] and SPC/E [51] water models in combination with the LJ model for methane proposed by Guillot and Guissani [60], but they followed a different approach to determine the coexistence equilibrium temperatures and pressures of the methane hydrate. They focused mainly on cases in which the hydrate melts, a process that takes only a few nanoseconds and thus occurs more rapidly than growth. Their results obtained using TIP4P/Ice agree with those of Jensen et al. and although they disagree with the results of Conde and Vega for the TIP4P/Ice model, they agree in the case of the TIP4P/2005 model. The TIP4P/Ice model combined with a united-atom LJ potential for methane is an adequate choice for classical MDsimulations of methane hydrates. [14,15] This combination was recently used [61] to explore the dissociation of the sI methane hydrate, induced by high-intensity (2–5Vnm−1) electric fields. The same combination of models and the same direct coexistence method, used in Refs. [9,59], were employed in this work to extend the studies on their ability to reproduce the experimental hydrate liquid water gas coexistence temperatures under various different pressures, and to explore the shift produced by an external constant electric field, using intensities in the range from 0.1 to 0.9Vnm−1. 2. MD methods In the simulations in this study we used the rigid/nonpolarisable TIP4P/Ice water potential. [53] This water mode has an LJ interaction site located on the oxygen atom, positive charges located at the positions of the H atoms, and a negative charge located at a distance dOMfrom the oxygen along the H–O–H bisector. For the methane, we used a single-site model proposed by Guillot and Guissani [60] consisting of an LJ 12-6 interaction. The potential parameters for these models are given in Table 1. The water potential was chosen because of its good performance in reproducing of the three-phase coexistence line of a binarymixture of water and methane, as observed by Conde and Vega [9]. The simulations were implemented using the GromacsMDpackage, version 4.1.5. [62–64] For the long-range Coulombic (C) interaction, the PME algorithm was used with a cut-off radius of 0.9 nm, an LJ interaction was implemented with a cut-off radius of 0.9 nm, and the Lorentz–Berthelot mixing rules were implemented in both cases. The simulations were implemented as NPT MD simulations using three-dimensional periodic boundary conditions. First, to determine the coexistence temperatures, we used equilibrium MD simulations, and we analysed the evolution of the potential energy as a function of time in simulations near the experimental coexistence temperature. We plotted the potential energy vs. time. The temperature was maintained using aNosé–Hoover thermostat [65,66]. We used a ParrinelloRahman barostat tomaintain a constant pressure. The system was equilibrated for 50 ps, and the three different sides of the simulation box were allowed to fluctuate independently. The time step used in the simulations was 2 fs. Our system consisted of three phases (hydrate, liquid water and gaseous methane). The hydrate crystallite structure was based on the results of the X-ray study performed by McMullan and Jeffrey [67] to investigate the type I structure of ethylene oxide hydrate, which provides the positions of the oxygen atoms in every water molecule and the mass centres of the unitary cells for the methane molecules. The orientations of the hydrogen atoms in the water molecules were randomly set such that they satisfied the Bernard–Fowler rules and the dipole moment of the system was nearly zero. The unit cell, with a side length of 1.203 nm, was replicated 2 × 2 × 2 times to form a cubic cell with a side length of 2.406nm that included 368 molecules of water and 64 molecules of methane; therefore, the ‘cages’ of water in the methane hydrate were completely occupied by methane molecules. The initial configuration for our simulations consisted of such a cubic methane hydrate crystallite surrounded, on its xy side, by a cubic cell with a side length of 2.406nm containing 368 water molecules and, on its −xy side, by a cubic cell with a side length of 2.406nm containing 64methane gas molecules. Therefore, the dimensions of this configuration were 2.406nm×2.406nm×7.218 nm. This configuration was equilibrated for 50 ps to obtain our initial configuration for all simulations; see Figure 1. First, our simulations were performed with the temperature of the system changing from 280 to 310Kat a pressure of 400 bar and from 280 to 315K at a pressure of 250 bar, without any external electric field, to calculate the coexistence temperatures at these pressures. Once the coexistence temperatures were calculated, we performed a series of simulations at 400 and 250 bar, with temperatures varying from near the coexistence temperature to temperatures far above the coexistences temperature; these simulations included an external electric field. For these simulations, the expression for the total potential energy includes the following terms [68]: U total =_ a<b U ab − Ee_ i q i z i , (1) where E is the field strength, z i is the z Cartesian coordinate of partial charge i, and U ab includes both the C and LJ contributions. The strength of the applied external electric field was in the range of 0.05–1.2Vnm −1 . Figure 3. (Colour online) Experimental and calculated values from this work and from the literature for the three-phase coexistence temperature (Conde and Vega [9], Jensen et al. [7], Michalis et al. [59]). The experimental results were taken from Sloan et al. [1] and all authors used the TIP4P/Ice water model. Figure 4. (Colour online) Number of hydrogen bonds as a function of time for a system at 400 bar of pressure and different temperatures. Table 1. Potential parameters for the methane molecule model [60] and the TIP4P/Ice water model [53] used in this work. 3. Results and discussions 3.1. Calculation of the coexistence temperature without an external electric field This section presents the results of our calculations of the three phase coexistence temperature using the direct coexistence method. For the initial temperatures, we selected 295K for a pressure of 400 bar and 280K for a pressure of 250 bar because these points correspond to stable hydrates in the phase diagram. Moreover, at these points, the simulations of the three phases evolve in such a way that the hydrate phase occupies the entire system and form within a relatively short time; afterward, a series of simulations at increasing temperatures was performed to observe the evolution over time of the (LJ + C) potential energy for each temperature. A gradual increase in potential energy indicated melting, whereas a decrease indicated crystallisation. The resulting equilibrium phase coexistence temperature was calculated as the average of the lowest temperature at which the hydrate melted and the highest temperature at which the system froze. In Figure 2(a), we show the time evolution of the potential energy for the TIP4P/Ice water model at a pressure of 250 bar for the reproduction of the coexistence temperature, resulting in a value of 287.5 K, which is close to the experimental temperature (291.1 K) for that pressure. In Figure 2(b), we show the performance of the TIP4P/Ice model at a pressure of 400 bar. This result is in good agreement with the results of Conde and Vega. For the TIP4P/Ice system, the coexistence temperature is 297.5Kat 400 bar, which is close to the experimental coexistence temperature (297.2 K). For this reason, we decided to use the system with the TIP4P/Ice water model and a single LJ centre for methane to study hydrate formation under an external electric field. The direct coexistence method possesses an inherent degree of stochasticity,[59,69] which is why our result at 400 bar is slightly different from the results of Conde and Vega (−5K) and Michalis et. al. (+4.1 K). This can be observed in Figure 3, where the logarithm of the pressure is plotted as a function of temperature. The experimental data were taken from Ref. [1] To analyse the crystallisation process, we measured the number of hydrogen bonds throughout the simulation time. In order for a hydrogen and an acceptor were considered hydrogen bonded, two criteria had to be satisfied: (1) the distance between the hydrogen and the acceptor had to be less than 0.35nm , and (2) the hydrogen-donor–acceptor angle had to be less than 30 ◦ ; when the system crystallises, this number increases and tends to approach the maximum possible value (two hydrogen bonds per water molecule), whereas if the system does not crystallise, the number of hydrogen bonds remains constant or tends to decrease slightly. Figure 4 shows the increase or decrease in the number of hydrogen bonds for a system at 400 bar. For temperatures below the coexistence temperature, the number of hydrogen bonds increases to reach almost 1472, and for temperatures above the coexistence temperature, the number of hydrogen bonds remains near 1380. The increment in the number of hydrogen bonds is directly associated with a decrement in the potential energy of the system; see Figure 2. The order that emerges in the incrementing of the number of hydrogen bonds has no important consequences for the total dipole moment of the system. In Figure 5(a), we have plotted both the components and the norm of the total dipole moment vector of the system during a simulation of a system at 400 bar and 290 K. In (b), we present the results for a system at 400 bar and 310 K. In both cases, the components of the total dipole moment remain nearly constant, with a norm close to zero. This means that even when the system has crystallised or melted, the total dipole moment does not exhibit any significant changes. 3.2. Calculation of the coexistence temperature under an external electric field The main issue to be analysed in this work is the effect of an applying external electric field during the formation ofmethane hydrate, that is, how this field can modify the three-phase coexistence temperature of water, methane gas and hydrate. We performed two series of runs at two fixed pressures (250 and 400 bar) at temperatures similar to and higher than the coexistence temperature that included an external electrostatic field the magnitude of which varied from 0.05 to 1.1Vnm−1. We began increasing the temperature from a point slightly above the coexistence temperature and raised it to temperatures considerably higher than the coexistence temperature, all under the applications of an external electric field. At the pressure of 400 bar, the first temperature at which the system did not crystallise without the electric field was 300 K; see Figure 2(b). However, we observed a sudden crystallisation when we applied an external field. This crystallisation occurred only for a certain range of electric fields. Figure 6(a) shows the time evolution of the potential energy of a system at 400 bar and 300 K; for electric field magnitudes between 0 and 0.1Vnm−1, the system evolved such that the hydrate phase melted. Then, there was a range of electric field amplitudes (0.2–0.9Vnm−1) in which the system evolved such that the hydrate phase grew to occupy the entire system. Finally, for electric field magnitudes above 0.9Vnm−1 the system suffered a rapid melting that lasted for only the first 2 ns of the simulation. It could be hypothesized that the observation of this effect might be a result of using the direct coexistence method, which possesses an inherent degree of stochasticity,[59,69] However, this stochasticity applies only in a small region near the coexistence temperature; however, when we continued increasing the temperature at the same pressure in the presence of an electric field, the same effect was observed at 305K even though the range of magnitudes at which the system crystallised was reduced to 0.3–0.8Vnm−1; see Figure 6(b). Even at 310 K; see Figure 6(c). And at 315 K; see Figure 6(d)], this effect was observed. In the case at 310 K, the range was 0.3–0.7Vnm−1, and for the case of the system at 315 K, the range was 0.3–0.7Vnm−1. Finally, the effect was no longer present when we explored the evolution of potential energy of the system at the temperature of 320 K; see Figure 7. In this case, all simulations that included an electric field at 320K and 400 bar exhibited an increasing potential energy, with consequent melting of the system. To analyse the effect of the electric field on the number of hydrogen bonds formed during crystallisation, we measured this number throughout the simulation. The electric field attenuated the thermal effect on the molecules, thereby promoting the formation of hydrogen bonds. Figure 8 shows that for certain values of the electric field, there was an increment of the number of hydrogen bonds, whereas if the field was not applied, the system suffered a decrement in this number. The effect of the electric field on the system includes an alignment of the water dipole moments in the direction of the field; this effect of alignment competes with the thermal effect, enabling the formation of more hydrogen bonds. In Figure 9, we present a plot of the total dipole moment of the system as a function of time for a system that crystallises at 400 bar, 310K and 0.7Vnm−1. When the system crystallises under the influence of an external electric field, the system retains a net dipole moment aligned to the direction of the field. The same procedure was followed for a pressure of 250 bar, we began increasing the temperature from a point slightly above the coexistence temperature and raised it to temperatures considerably higher than the coexistence temperature, all under the applications of an external electric field. At 250 bar, the first temperature at which the system did not crystallise in the absence of an electric field was 290 K; see Figure 2(a). However, the same effect was observed when we included the electric field in our simulations, that is, a range of electric fields existed in which the effect of the thermal vibrations was attenuated by the effect of the electric field, allowing crystallisation of the system to occur. At 290 K, this range was from 0.1 to 0.9Vnm−1; see Figure 10(a). The same phenomenon was also observed in our simulations at 295, 300, 305 and 310 K; see Figure 10(b)–(e), respectively. The range in which this phenomen on was observed was 0.2–0.8Vnm−1 at 295 K, 0.3–0.9Vnm−1 at 300 K, 0.4– 0.9Vnm −1 at 305 K, and 0.4–0.8Vnm −1 at 310 K. Finally, at 315 K, all simulations exhibited an increasing potential energy, with the consequent melting of the system. If we consider the most efficient electric field magnitude and calculate the coexistence temperature as the average of the highest temperature at which the system crystallised and the lowest temperature at which the system melted, then the simulations suggest that an external electric field can shift the coexistence temperature from 297.5K (+/−2.5 K) to 317.5K (+/−2.5 K) at 400 bar and from 287.5K (+/−2.5 K) to 312.5K (+/−2.5 K) at 250 bar, that is 22.5K on average. This result is comparable with the results of Aragones et al. [47] for ices III and V which increase the melting point by about 15K by means of a field of 0.3Vnm −1 . Figure 11 shows the proposed electric field–temperature phase diagrams at 250 and 400 bar, which include the region in which the hydrate exists in the system under the application of different electric fields. We can see a small variation of the range of field and of magnitude of the temperature shift in both pressures studied in this work, so maybe this effect could depend on pressure. Figure 5. (Colour online) Graphs of the components and norm of the total dipole moment of the system. In (a), the system crystallises at 400 bar and 290 K. In (b), the system melts at 400 bar and 310 K. Figure 6(Colour online) Evolution of the potential energy at a fixed pressure of 400 bar for different electric field magnitudes and temperatures. (a) Shows the results at 300 K; the range in which the system crystallises is 0.2–0.9Vnm −1 . (b) Shows the results at 305 K; the range in which the system crystallised is 0.3–0.8Vnm −1 . (c) Shows the results at 310 K; the range in which the system crystallised is 0.3–0.7Vnm −1 . (d) Shows the results at 305 K; the range in which the system crystallised is 0.3–0.7Vnm −1 . Figure 7. (Colour online) Graphs of the evolution potential energy at a fixed temperature and pressure of 320 K and 400 bar, respectively, under the application of electric fields of different magnitudes. Figure 8. (Colour online) Graphs of the evolution of the number of hydrogen bonds at a fixed temperature and pressure of 310 K and 400 bar, respectively, under the application of electric fields of different magnitudes. Figure 9. (Colour online) Graphs of the components and norm of the total dipole moment of the system. In this case, the system crystallises at 400 bar, 310 K and 0.7Vnm −1 . Figure 10. (Colour online) Evolution of the potential energy at a fixed pressure of 250 bar for different electric field magnitudes and temperatures. (a) Shows the results at 290 K; the range in which the system crystallises is 0.1–0.9Vnm −1 . (b) Shows the results at 295 K; the range in which the system crystallised is 0.2–0.8Vnm −1 . (c) Shows the results at 300 K; the range in which the system crystallised is 0.3–0.9Vnm −1 . (d) Shows the results at 305 K; the range in which the system crystallised is 0.4–0.9Vnm −1 . (e) Shows the results at 310 K; the range in which the system crystallised is 0.4–0.8Vnm −1 . (f) Shows the results at 315 K; where the potential energy of the system only increases over time and, in consequently, the hydrate is melted throughout the entire range of electric fields. The results between the maximum and minimum values of the electric field at which the system crystallises are not shown for simplicity. Figure 11. Proposed electric field–temperature phase diagram at 250 and 400 bar with errors of (+/−2.5 K). 4. Conclusions The direct phase coexistence method was implemented in this work to study the effect of an external electric field on the threephase coexistence temperature of a binary mixture of water and methane. Two different pressures were examined: 250 and 400 bar. The effects of an external electric field in the range of 0.1–0.9Vnm −1 , which includes the ordering of the water molecules in the direction of the field, caused the effect of the thermal vibration of the water molecules to become attenuated. This resulted in a shifting of the three-phase coexistence temperature to higher temperatures. Electric fields with magnitudes below this range did not cause a differences in the coexistence