Prediction of the phase equilibria of methane hydrates using the direct phase coexistence methodology.

Prediction of the phase equilibria of methane hydrates using the direct phase coexistence methodology Vasileios K. Michalis, Joseph Costandy, Ioannis N. Tsimpanogiannis, Athanassios K. Stubos, and Ioannis G. Economou Citation: The Journal of Chemical Physics 142, 044501 (2015); doi: 10.1063/1.4905572 View online: http://dx.doi.org/10.1063/1.4905572 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/142/4?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Molecular dynamics simulation of CO2 hydrates: Prediction of three phase coexistence line J. Chem. Phys. 142, 124505 (2015); 10.1063/1.4916119 Note: A simple correlation to locate the three phase coexistence line in methane-hydrate simulations J. Chem. Phys. 138, 056101 (2013); 10.1063/1.4790647 Molecular dynamics simulations of vapor/liquid coexistence using the nonpolarizable water models J. Chem. Phys. 134, 124708 (2011); 10.1063/1.3574038 Molecular-dynamics evaluation of fluid-phase equilibrium properties by a novel free-energy perturbation approach: Application to gas solubility and vapor pressure of liquid hexane J. Chem. Phys. 124, 124111 (2006); 10.1063/1.2178321 Solid solutions Ne –n D 2 . Diagram of phase equilibrium Low Temp. Phys. 31, 947 (2005); 10.1063/1.2127875

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.129.120.3 On: Wed, 03 Jun 2015 12:58:40

THE JOURNAL OF CHEMICAL PHYSICS 142, 044501 (2015)

Prediction of the phase equilibria of methane hydrates using the direct phase coexistence methodology Vasileios K. Michalis,1 Joseph Costandy,1 Ioannis N. Tsimpanogiannis,2 Athanassios K. Stubos,2 and Ioannis G. Economou1,a) 1

Chemical Engineering Program, Texas A&M University at Qatar, P.O. Box 23847, Doha, Qatar Environmental Research Laboratory, National Center for Scientific Research NCSR “Demokritos,” Aghia Paraskevi, Attiki GR-15310, Greece 2

(Received 2 November 2014; accepted 18 December 2014; published online 22 January 2015) The direct phase coexistence method is used for the determination of the three-phase coexistence line of sI methane hydrates. Molecular dynamics (MD) simulations are carried out in the isothermal–isobaric ensemble in order to determine the coexistence temperature (T3) at four different pressures, namely, 40, 100, 400, and 600 bar. Methane bubble formation that results in supersaturation of water with methane is generally avoided. The observed stochasticity of the hydrate growth and dissociation processes, which can be misleading in the determination of T3, is treated with long simulations in the range of 1000–4000 ns and a relatively large number of independent runs. Statistical averaging of 25 runs per pressure results in T3 predictions that are found to deviate systematically by approximately 3.5 K from the experimental values. This is in good agreement with the deviation of 3.15 K between the prediction of TIP4P/Ice water force field used and the experimental melting temperature of ice Ih. The current results offer the most consistent and accurate predictions from MD simulation for the determination of T3 of methane hydrates. Methane solubility values are also calculated at the predicted equilibrium conditions and are found in good agreement with continuum-scale models. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4905572]

I. INTRODUCTION

The clathrate hydrates are solid solutions of water with other relatively small molecules, which exhibit crystal structure.1 The water molecules, interconnected through hydrogen bonds, form a three-dimensional crystal lattice with welldefined cages inside which the “guest” molecules are found. The guest molecules provide the stability to the hydrate structures and destruction of the crystal occurs (i.e., hydrate dissociation) when the guest molecules are removed from the structure either by increasing the temperature or lowering the pressure. Note, however, that the existence of metastable empty cages has been reported recently.2–4 There are more than 1301 known guests that form hydrates, such as methane, ethane, propane, iso-butane, carbon dioxide, hydrogen sulfide, nitrogen, and hydrogen. The size, as well as the interactions of the guest molecule with water, defines the type of crystal structure. The three most common structures are sI, sII, and sH which primarily differ in the number and type of formed cages. Under conditions of interest to practical applications, the methane hydrates, on which the present study is focused, form sI structure.1 Each unit cell consists of 2 pentagon dodecahedron cages (512) and 6 tetrakaidecahedron cages (51262). This structure belongs to the space group Pm3n with a lattice constant of about 1.2 nm. Methane hydrates have also been shown experimentally to undergo structural transitions at higher pressures. In particular, a)Author to whom correspondence should be addressed. Electronic mail:

[email protected] 0021-9606/2015/142(4)/044501/12/$30.00

Chou et al.5 reported an sI-to-sII structural transition at approximately 1000 bar and subsequently an sII-to-sH structural transition at approximately 6000 bar. A detailed discussion of all the experimental studies that are related to structural transitions is presented in the review by Loveday and Nelmes.6 Clathrate hydrates, although extensively studied as is clearly indicated by the continuously increasing number of studies per year,7 remain at the center of interest for various scientific as well as technological reasons. Natural gas hydrates are of particular importance as they play a major role in various technological and natural processes. First and foremost, the oil and gas-industry is interested in them from the flow assurance point of view as they tend to form at a pressure and temperature range which often lies inside the operational limits, resulting in blocking of pipelines and posing a severe safety threat.8–10 Understanding thus the phase behavior of these systems, either on their own or in combination with various inhibitors, is of primary importance. The fact that the volumetric density of most gases in hydrate form corresponds to that of a condensed state provides the ground for the exploration of the idea of using hydrates for gas storage and transportation purposes.11–16 Also, differences in structural and thermodynamic properties of different molecules in the hydrate form can potentially be exploited for the design of alternative separation processes for energy17,18 or environmentally related19–21 gas mixtures. Furthermore, hydrate growth results in ion exclusion from the crystal structure, a process that can be used for water purification and desalination.22,23

142, 044501-1

© 2015 AIP Publishing LLC


044501-2

Michalis et al.

The methane hydrate is a naturally occurring substance which is found under the oceanic floor and at the permafrost regions.24 In addition, the released methane, that would eventually reach the atmosphere, is a greenhouse gas more potent than carbon dioxide so that emissions from methane hydrate reservoirs could be assessed as a climate change risk.25,26 Additionally, the sudden dissociation of methane hydrates could give rise to geologic hazards such as the collapse of oceanic slopes.27,28 Last but not least, the methane hydrates draw attention as a potential energy source.29,30 The gradual decline of the conventional hydrocarbon reserves and the emergence of new technologies offer the possibility of access to unconventional deposits. It is estimated31–33 that the energy content of the methane hydrates around the world far exceeds the respective content of the known hydrocarbon reservoirs. Although exploitation of such reservoirs currently remains economically unfeasible and technologically challenging, it is difficult to ignore their future potential importance.34,35 Molecular simulation techniques play an increasingly important role in science today offering critical insight and improving our understanding of the properties of fluids and materials. Computational methods such as molecular dynamics (MD) and Monte Carlo (MC)36 simulations are nowadays used extensively in order to predict structural, thermodynamic, and transport properties of a wide range of materials. Both MD and MC have been used to study gas hydrates in parallel with experimental studies and equation of state approaches. Although the van der Waals–Platteeuw theoretical approach37 and its variations38–43 offers a satisfactory description of the phase equilibria for the most common hydrate systems, it is nevertheless intriguing to approach the hydrate growth and dissociation processes with molecular simulation methodologies as it can lead to a deeper understanding of the related microscopic phenomena. If the developed methodologies are proven to be accurate enough, they can be used as a tool for the design of new processes in a complementary role to the experimental and continuum-scale approaches. MD simulations can provide a fundamental understanding of processes at the molecular level. In particular, they can contribute in obtaining a better understanding of the mechanisms involved during the nucleation, growth, or dissociation of hydrate structures, as well as for the calculation of thermodynamic and transport properties. The study of the methane hydrate in particular holds a central role due to its major technological importance and the availability of experimental results for validation purposes. An accurate knowledge of the thermodynamic and kinetic properties of methane hydrate is a primary goal for a number of studies. MD studies considering methane hydrates have been reported as early as 1983 with the pioneering work of Tse et al.44 Subsequently, a significant number of MD studies examining methane hydrates have been reported focusing on topics such as hydrate nucleation,45–50 formation and growth,51,52 dissociation,53–58 inhibition,59–61 heat transfer effects,62–64 the effect of hydrate contact with solid surfaces,65–68 methane replacement by carbon dioxide,69–72 and the calculation of transport properties such as thermal conductivity,73–75 or thermodynamic properties such as three-phase equilibrium conditions.76–79

J. Chem. Phys. 142, 044501 (2015)

In the area of calculation of the three-phase equilibrium conditions for pure methane hydrate which is the main focus of the current study, there is a number of molecular simulation studies that have provided important insights of the mechanisms involved. Conde and Vega76,80 studied the methane hydrate three-phase equilibrium with the direct phase coexistence method. In this method using MD simulations, a three phase system consisting of a hydrate slab, a liquid water slab, and a gas methane slab are brought in contact and the evolution of the system is monitored under specified pressure and temperature. The authors showed that by varying the temperature for a given pressure, the three-phase coexistence temperature can be determined. Using a simple Lennard-Jones model for the methane molecule, they assessed the behavior of different water force fields and found that the determined coexistence temperature depends heavily on the predicted melting temperature of ice for each force field. In their simulations, they observed the formation of methane bubbles which was attributed to system size effects. Although their results exhibited significant uncertainty, they showed that the TIP4P/Ice81 offers the best predictions. Tung et al.77 used the same methodology but with a different initial configuration and by using the TIP4P/Ew82 force field for the water with the OPLS-All Atom (OPLSAA)83 force field for the methane, with geometric combing rules, and using in particular for the C–O interaction the parameters proposed by Cao et al.84 Jensen et al.85 calculated the three phase equilibrium conditions using the Monte Carlo simulation method in order to calculate the chemical potential of water and methane in the hydrate phase and in the liquid and gas phase, respectively. They used the TIP4P/Ice force field with the OPLS-United Atom (OPLS-UA)86 for methane. Their findings systematically over-predict the equilibrium temperature for a given pressure and disagree with Conde and Vega.76 Smirnov and Stegailov79 followed a different approach for the determination of the equilibrium conditions of the methane hydrate. They focused mainly on the cases where the hydrate dissociates, a process which is much faster than growth. They investigated TIP4P/Ice,81 TIP4P/2005,87 and SPC/E88 in combination with the Lennard-Jones model for methane proposed by Guillot and Guissani.89 Their results using TIP4P/Ice agree with Jensen et al.85 and while they disagree with Conde and Vega76 for the case of TIP4P/Ice, they agree for the case of TIP4P/2005. The main goal of this work is the determination of the three-phase coexistence line of sI methane hydrates using the direct phase coexistence methodology90 avoiding supersaturated conditions. Supersaturated conditions are not commonly encountered in nature, while on the other hand, methane concentration that corresponds to the solubility limit is believed to offer improved accuracy. Additionally, multiple independent runs per pressure and temperature have been used to compensate for the effect of stochasticity in the determination of the coexistence conditions. Finally, various insights are obtained in the study of hydrate growth and dissociation with MD simulations. The paper is organized with the following structure. Initially, the direct phase coexistence methodology is briefly described along with structure preparation, force field descrip-


044501-3

Michalis et al.

J. Chem. Phys. 142, 044501 (2015)

tion, and computational details. Following, the main findings of this work are presented and discussed for the prediction of the three-phase coexistence line of the methane hydrate. Finally, we end with the conclusions.

II. METHODOLOGY

In the direct phase coexistence methodology, different phases including a solid hydrate, liquid water and gas methane are brought in contact, and through MD NPT simulations, the system is allowed to evolve to the equilibrium state that corresponds to the given conditions. By scanning the temperature for a given pressure, it is possible to determine the coexistence temperature. The three-phase configuration used in this study consists of a solid hydrate slab, two liquid water slabs that surround the hydrate slab and one gas methane slab which lies between the water slabs, giving in total four slabs in a liquid water – solid hydrate – liquid water – gas methane (WHWM) arrangement. Each slab was initially separately equilibrated for each different pressure and subsequently the four slabs were connected with 0.1 nm buffer distance, followed by energy minimization via steepest descent to avoid bad contacts, while keeping the oxygen positions frozen. The slabs were connected with the interfaces normal to the xdirection. While equilibrating the water and methane slabs, their y and z dimensions were kept equal to the equilibrated y and z dimensions of the hydrate slab. The hydrate slab consists of a 2 × 2 × 2 sI hydrate supercell. The positions of the oxygen atoms within the unit cell of the sI hydrate were obtained from McMullan and Jeffrey.91 All the centers of the cages were filled with methane molecules, giving a total of 368 water and 64 methane molecules. Essentially, we examined the case of 100% cage occupancy. Finally, the positions of the hydrogen atoms were found by minimizing the energy of the system while fixing the oxygen atoms, resulting in a structure that respects the Bernal and Fowler rules.92 This approach has been used93 as an alternative to using a configuration with minimum dipole moment. Each water slab contained 368 water molecules while the methane slab contained 64 molecules for pressure values of 40 bar and 100 bar and 128 molecules for pressures of 400 bar and 600 bar. This variation in the number of the methane molecules was needed to compensate for the change of methane density and keep the length of the methane slab big enough to avoid bubble formation. A snapshot of the initial configuration at 100 bar is shown in Fig. 1(a). The TIP4P/Ice force field81 was used to represent the water and the OPLS-UA force field parameters86 were used in order to represent the methane as a single Lennard-Jones

FIG. 1. Three snapshots of a typical trajectory of the WHWM system at 294 K and 600 bar (run no. 4 of Fig. 2): (a) initial configuration at t = 0 ns, (b) intermediate step at t = 600 ns, and (c) final state at t = 1500 ns. The red and white lines represent the water molecules and the blue spheres represent the methane molecules.

sphere. The TIP4P/Ice is a rigid, non-polarizable model with one Lennard-Jones interaction site on the oxygen position. Positive charges are allocated on hydrogen sites while the negative charge is positioned at a distance of d OM = 0.1577 Å from the oxygen atom along the bisector of the H–O–H angle. The TIP4P/Ice model was developed in order to reproduce the melting temperature of hexagonal ice (Ih) predicting a melting temperature of 270(3) K,94 which according to the study of Conde and Vega76 offers an advantage to the calculation of the three-phase coexistence temperature for the hydrate systems. The Lorentz–Berthelot36 combining rules were used for the cross-interaction parameters. The atomic charges, Lennard-Jones parameters, and interatomic distances are summarized in Table I. All simulations were run in the isothermal–isobaric (NPT) ensemble using the GROMACS MD simulation package (version 4.6.3).95–97 The pressure and temperature coupling were implemented using the Berendsen scheme98 with time constants of 1 ps. The pressure coupling used was anisotropic with the same compressibility in all three dimensions, such that each side of the simulation box can fluctuate independently. The leap-frog integrator scheme was used with a time step of 2 fs. The total simulation time depended on pressure and temperature conditions, so that the total time was sufficient to provide clear conclusions about the final state of the system and was typically in the range of

TABLE I. Potential parameters of TIP4P/Ice (water)81 and OPLS-UA (methane)86 models. The interatomic distances (d), angle, charges (q), and the Lennard-Jones parameters σ and ε are given for all the atoms and the negative charge site of the water model (M). k B is the Boltzmann constant.

TIP4P/Ice OPLS-UA

d O H (Å)

H –O–H angle

q H (e)

q M (e)

d OM (Å)

σ (Å)

ε/k B (K)

0.9572

104.52

0.5897

−1.1794

0.1577

3.1668 3.73

106.1 148


044501-4

Michalis et al.

J. Chem. Phys. 142, 044501 (2015)

FIG. 2. Potential energy vs time for five independent NPT simulation runs of 1500 ns at 294 K and 600 bar. The inset figure focuses on the second stage of hydrate growth for the case of run no. 2.

1000–4000 ns. The long range Coulombic interactions were handled with the PME method.99 Both the Coulombic and the Lennard-Jones interaction potentials were calculated with a cutoff distance of 11 Å. Multiple independent runs at the same conditions of temperature and pressure were implemented using different random seeds for the initialization of the velocity. Periodic boundary conditions were applied in all directions.

III. RESULTS AND DISCUSSION

In Fig. 1, three different snapshots are presented that depict the initial (0 ns), an intermediate (600 ns), and the final state (1500 ns) of a typical run at 294 K and 600 bar. Figure 2 presents the time evolution of the potential energy for the case of 294 K and 600 bar for five independent runs for a simulation time of 1500 ns each. This temperature corresponds to hydrate growth conditions and is 7.7 K below the experimental equilibrium value.100 Growth is indeed observed in all the five cases which is also clear from the decrease of the potential energy up to the equilibrium value that corresponds to the transition of the whole system to the hydrate phase. The decrease of the potential energy in this case can be separated in two distinct stages which are characterized by different growth rates. The first one, starting from the beginning of the simulation until the onset of the second stage is the “normal” growth stage, where the methane molecules of the gas phase initially dissolve, subsequently slowly diffuse through the liquid water slabs, and eventually reach the water–hydrate interface, where hydrate growth takes place. The beginning of the second stage is characterized by a sudden decrease of the potential energy (as is clearly indicated in the inset of Fig. 2 for the case of run no. 2) and it ends with the full conversion of the system to the hydrate phase. After this point, the potential energy of the system remains practically constant. This accelerated growth rate is

due to the formation of a methane bubble which leads to locally supersaturated conditions. When enough methane and water molecules are removed from the gas and liquid phase, respectively, by being incorporated in the hydrate phase, then two separate possible phenomena may occur (1) either the width of the water phase is diminished allowing the presence of the methane molecules in the vicinity of the hydrate–water interface or (2) the width of the methane slab becomes small enough and comparable to the width of the fluctuations of the water–methane interface, thus allowing the methane slab to enter the liquid phase in the form of a bubble. This bubble can be either spherical or cylindrical depending on its position with respect to the periodic boundaries. An illustrative example of a formed methane bubble is given in Fig. 3. Similar behavior has been observed in a number of studies.49,76 It should be noted that under the conditions of the particular example, the methane is in the supercritical state (Pc = 45.992 bar, Tc = 190.564 K

FIG. 3. Snapshots of a trajectory of the WHWM system at 290 K and 400 bar where (a) a thin methane slab is shown at t = 972.6 ns, and (b) 0.7 ns later a methane bubble has been formed. The red and white lines represent the water molecules and the blue spheres represent the methane molecules.


044501-5

Michalis et al.

J. Chem. Phys. 142, 044501 (2015)

FIG. 4. Time evolution of the potential energy of five different runs at 283 K and 100 bar near equilibrium conditions.

(Ref. 101)) and the term “gas” is used to distinguish it from the liquid water slab. When supersaturated conditions are present, the hydrate growth rate is significantly higher than in the case of the first stage where the methane concentration in the water is close to or at the solubility limit. In any case, in all of our simulations, due to the choice of the initial configuration, the “normal” growth stage was always present and it offered enough indication for the existence of hydrate growth, without the need of supersaturated conditions. It is apparent from Fig. 2 that the “normal” growth rates of the five different runs are not equal, as indicated by the different slopes of the potential energy, a fact that must be attributed to the stochastic nature of the hydrate growth. Apart from the short-term potential energy fluctuations, there are fluctuations that expand over hundreds of nanoseconds. This case corresponds to alternating phenomena of hydrate growth and dissociation, although in the long term the system can be on average moving to one of the two final states. A direct consequence of this fact is that for the calculation of growth rates through MD simulations, averages should be used over an adequate number of runs. Additionally, the “induction” time for the bubble formation is not equal between different runs. The effect of stochasticity that is inherent to the direct phase coexistence methodology is more pronounced near the three-phase equilibrium conditions. Indicative results for the case of 100 bar and 283 K are presented in Fig. 4. For the given pressure of 100 bar, the experimental threephase coexistence temperature is T3,experimental = 286.2 K. The temperature of 283 K is very close to expected equilibrium temperature and in line with the work of Conde and Vega.76 These authors showed by comparing different water models that the deviation of the predicted T3 by each model from the experimental values is correlated with the deviation of the predicted melting temperature of ice by each model from the respective experimental value. In particular, in the case of TIP4P/Ice, the melting temperature is 270(3) K94 and thus

T3,expected = T3,experimental − 3.15 K. The potential energy graphs in Fig. 4 correspond to different runs at the same pressure and temperature and show that the system for the given conditions near the equilibrium can either dissociate or form hydrate. Slightly different initial velocities in combination with the fluctuations imposed by the temperature and pressure coupling can lead to different final states (i.e., growth, dissociation). This behavior is persistent within a few Kelvin around the equilibrium temperature and it was found in all of our simulations. Additionally, it can be observed that the hydrate dissociation rate is different between the different runs in an analogous manner with that of hydrate growth that was discussed previously. It should be noted that Espinosa et al.,102 who studied the direct coexistence of a solid phase with a liquid-like phase of a pseudo hardsphere fluid, observed a similar stochastic behavior close to coexistence conditions and additionally showed that the range of this behavior depends on the system size. Nevertheless, a very large system for the case studied in this work is computationally intractable. This stochastic behavior poses a severe problem in the determination of the coexistence temperature if only one run is used. Figures 5(a) and 5(b) present two variations of temperature scans for the determination of T3 at 100 bar. For the first case, shown in Fig. 5(a), given the information provided by the potential energy curves, the calculated T3 is 280 K. For the second scenario, presented in Fig. 5(b), a different set of runs has been used, resulting in a calculated T3 value of 286 K. This discrepancy originates from the use of a limited (in this case just one) number of runs for each temperature and the stochastic nature of this kind of simulations. An additional implication which is caused by the stochasticity of the hydrate growth and dissociation processes is that the simulation time that is needed in order to draw a clear conclusion on the final state of the system varies significantly between different runs. This can be problematic particularly for temperatures near the equilibrium value.


044501-6

Michalis et al.

J. Chem. Phys. 142, 044501 (2015)

(a)

FIG. 5. Time evolution of the potential energy of two variations of temperature scans at 100 bar. Although the pressure and temperatures are the same for both figures (a) and (b), different independent runs have been selected in each scenario. In the case of figure (a), the calculated T3 is 280 K while in the case of figure (b), it is 286 K.

(b)

In Fig. 6 five different runs are presented for the case of 296 K and 400 bar. Under these conditions, all the runs lead to dissociation but the necessary simulation time ranges between 600 and 2500 ns. For example, it would be misleading to end the simulation of run no. 5 at a time in the range of 1100–1200 ns as this would draw the erroneous conclusion of hydrate growth. Indicative temperature scans, as in the case of 100 bar, shown in Fig. 5, are presented for the remaining of the pressures examined in this work in Figs. 7(a), 7(b), and 7(c) for 40, 400, and 600 bar, respectively. Although, as already discussed, single run temperature scans cannot be used safely for the determination of T3, a comparison between the different pressures reveals various trends. First, the higher the pressure, the less the required simulation time in order

to draw a conclusion about the final state of the system, although for temperatures very close to T3, this may not hold. Second, the hydrate growth rate increases by decreasing the temperature and this holds for all the pressures examined. An opposite conclusion can be drawn if a comparison is made between a limited number of runs; then due to stochasticity, one could find a scenario where the higher temperatures correspond to higher rates. On average, nevertheless, the rate is increased by decreasing the temperature. Finally, in a similar manner, dissociation rate increases as temperature increases. Based on the aforementioned observations and in order to determine the three-phase coexistence temperature of the methane hydrate system, we performed for each pressure a temperature scan of five different temperature values. For


044501-7

Michalis et al.

J. Chem. Phys. 142, 044501 (2015)

FIG. 6. Time evolution of the potential energy of five independent runs at 296 K and 400 bar. The dashed line denotes the value of the potential energy at three-phase equilibrium state.

each temperature, we performed five runs, resulting in 25 simulations per pressure. We examined four different pressure values, namely, 40, 100, 400, and 600 bar, which cover a broad range of the three-phase coexistence curve of the methane hydrate. Each temperature scan was centered on the expected three-phase coexistence temperature with a step of 2 K. In line with the work of Conde and Vega,76 we assumed that the expected T3 is equal to the experimental value corrected by the deviation of the prediction of the TIP4P/Ice water model for the melting temperature of ice Ih. This prediction is 270(3) K;94 and therefore, we considered a correction value of 3.15 K. An example of the averaging method for the determination of T3 for the case of 100 bar is presented in Table II. The final state of each of the 25 runs is characterized as either resulting in hydrate growth (g) or hydrate dissociation (d), and the appropriate averages are calculated as the arithmetic mean of the closest temperature that corresponds to different final states. For this case of 100 bar, the calculated coexistence temperature is T3 = 282.8 K with a standard deviation of 3.2 K. The experimental value of T3 for the given pressure is 286.2 K and so there is a difference of 3.4 K that is very close to the difference between the experimental and the expected value. The results for all the four different pressure values examined in this work are presented in tabulated form in Table III and plotted in Fig. 8 along with the experimental values,100 the expected values, the results by Conde and Vega,76 and the results by Tung et al.,77 who also used the direct phase coexistence methodology. Additionally, we have included results from Jensen et al.85 and Smirnov and Stegailov.79 The results of the present work are in excellent agreement with the expected values (T3,expected = T3,experimental − 3.15 K) in the examined range of 40–600 bar. The obtained results exhibit the capability of the specific methodology to offer a consistent prediction of the three-phase coexistence conditions with increased accuracy. The predictions by Conde and Vega, although in qualitative agreement, are not

quantitatively consistent with the slope of the experimental values. On the other hand, the deviation of the results of Tung et al. are consistently on average 5.8 K lower than the experimental values. The predicted melting temperature of ice using the TIP4P/Ew water model is 244(3) K94 which is approximately 29 K below the experimental value. This result is not in agreement with the conclusions of Conde and Vega76 that the shift in the predicted T3 follows the deviation of the predicted melting temperature of ice by the water model in use. A probable reason for this discrepancy could be the use by Tung et al. of modified cross-interactions’ parameters between the oxygen atom of water and the carbon atom of methane. The synergetic effect on the calculation of T3 between the predicted melting temperature of a water forcefield and the methane-water interaction has been also pointed out by Conde and Vega.80 It would be interesting in a future work to examine if such modifications of the cross-interaction parameters lead to different and probably increased solubility of methane in water and if this aspect is more important than the effect of the water model of choice on the calculation of the coexistence temperature. Jensen et al.,85 who use a different methodology that is based on thermodynamic integration and Monte Carlo simulations with TIP4P/Ice for the calculation of chemical potentials, systematically overpredict T3. On the other hand, the results of Smirnov and Stegailov79 that appear to be in agreement with Jensen et al. are based on MD calculations of direct phase coexistence of methane hydrate with liquid water and without any methane gas phase. It is clear from visual inspection of the trajectories of numerous runs that below a critical number of methane molecules in the gas phase in contact with the liquid water, methane is forming a gaseous bubble (either spherical or cylindrical due to periodicity). This bubble is moving inside the liquid water and as soon as it comes in contact with the hydrate slab the hydrate growth becomes extremely fast. Locally, at the hydrate–water interface, the presence of the methane bubble corresponds to supersaturated conditions


044501-8

Michalis et al.

J. Chem. Phys. 142, 044501 (2015)

FIG. 7. Indicative temperature scans of time evolution of the potential energy for different pressure values of (a) 40, (b) 400, and (c) 600 bar.


044501-9

Michalis et al.

J. Chem. Phys. 142, 044501 (2015)

TABLE II. Statistical averaging of 25 simulation runs for the determination of the three-phase coexistence temperature at a pressure of 100 bar. The calculated temperature is T3 = 282.8 K with a standard deviation of 3.2 K. The final state of each run is denoted as (g ) for hydrate growth or (d) for hydrate dissociation.

TABLE III. Experimental and calculated three-phase coexistence temperatures along with standard deviations and deviation of calculated from experimental. Experimental values are taken from the polynomial fit reported by Moridis.100 P (bar)

T (K)

Number 1

Number 2

Number 3

Number 4

Number 5

279 281 283 285 287 T3 (K)

g g g g d 286

g g g d d 284

g g d d d 282

g g d d d 282

g d d d d 280

well above the solubility limit. We have carried out simulations of hydrate in contact with liquid water artificially supersaturated with methane. We used a 2 × 2 × 2 hydrate slab with 368 water and 64 methane molecules and a supersaturated liquid water slab with 368 water and 32 methane molecules. In this case, we observed hydrate growth at temperature values well above the experimental T3. For the case of 400 bar, it was found that dissociation occurred at temperatures above 333 K, which is 36 K above the experimental value. This fact can provide a possible explanation for the overprediction of T3 in Conde and Vega’s results at certain pressures. For example, at 400 bar, they presented hydrate growth at 300 K76 while in our simulations, where we tried to avoid bubble formation for as long as possible, only dissociation occurs at this temperature. Additionally, it can provide a possible explanation to the overprediction of T3 by Smirnov and Stegailov who do not use a methane gaseous phase in their methodology. In this case, it is possible that during the dissociation, the methane molecules from the hydrate locally increase the methane concentration in the liquid water over the solubility limit, contrary to our case where the methane molecules diffuse through the liquid water and enter the gas phase.

40 100 400 600

T3,experimental (K)

T3,this work (K)

∆T (K)

277.4 286.2 297.2 300.7

273.0(2.4) 282.8(3.2) 293.4(0.9) 297.0(0.0)

−4.4 −3.4 −3.8 −3.7

Although widely known, it is worthwhile mentioning that the direct contact of different phases via a finite interface renders the use of dispersion corrections inapplicable.36 The dispersion corrections are customarily used in MD simulations due to the unavoidable calculation of the attractive dispersion potential up to a cutoff value as the size of the simulation box is finite. This truncation of potentials, such as the Lennard-Jones, is compensated by the use of dispersion corrections, also named tail corrections, and it is based on the assumption that the simulated system is isotropic and homogeneous. This assumption means that the radial distribution function of the system for distances greater than the cutoff value is equal to one and thus an analytic expression for the corrections can be calculated. The dispersion corrections apply to energy and pressure, however, do not affect the trajectory of the system. In the case of the direct phase coexistence method, as described in this work, it is clear that dispersion corrections cannot be used as the system is anisotropic and inhomogeneous. In particular, the use of tail corrections in combination with pressure coupling produces misleading results which are manifested mainly in erroneous values of the gas phase density. In Fig. 9 the density profile is presented for a system of gaseous methane in equilibrium with liquid water at 100 bar and 280 K with and without dispersion corrections. The results presented here are for

FIG. 8. Experimental98 and calculated values from this work and from the literature for the three-phase coexistence temperature of the methane hydrate system (Conde and Vega,76 Tung et al.,77 Jensen et al.,85 and Smirnov and Stegailov79). All authors used TIP4P/Ice except Tung et al. who used TIP4P/Ew. The expected values presented in the figure are defined as T3,expected = T3,experimental − 3.15 K.


044501-10

Michalis et al.

J. Chem. Phys. 142, 044501 (2015)

FIG. 9. Density of methane in a methane– water–methane configuration at 280 K and 100 bar, with and without dispersion corrections. The dashed line depicts the experimental value of the pure methane reported by NIST.101

a methane–water–methane configuration with 64–368–64 molecules, respectively. The MD calculated bulk density of the methane in the case where dispersion corrections are used is approximately 139 kg/m3, while in the case without dispersion corrections, the calculated density is 78 kg/m3. The experimental density101 for these conditions of the pure methane is 85.5 kg/m3. Although the particular conditions of this example are within the hydrate growth regime, similar behavior is observed in a wide temperature and pressure range. Effectively, the dispersion corrections in combination with pressure coupling are causing compression of the gas phase, an effect which is more pronounced for lower values of pressure according to our observations. Without the use of dispersion corrections, which is the case for all the simulations presented in this work for the determination of the three phase coexistence temperature, there is a persistent error in the calculated densities of the

order of approximately 5%–10% for the gas phase and 1% for the liquid water phase. Nevertheless, it appears that the magnitude of these deviations does not affect the calculation of the three-phase hydrate equilibrium conditions. The solubility of methane in the aqueous phase, along the three-phase equilibrium curve, as calculated by the current MD simulations, has also been estimated for the cases of the four pressures that were considered in the current study. The calculated MD results are shown in Fig. 10 where we plot the solubility of methane as a function of T3 as calculated in this work. Shown also are the solubility calculations by two different continuum-scale models1,103 as a function of the equilibrium temperature as calculated by those models. A detailed discussion on the performance of such thermodynamic models in predicting methane solubility in water has been reported recently by Tsimpanogiannis et al.104 One should recall that the T3 of the MD simulations

FIG. 10. Methane solubility in water as a function of the three-phase coexistence temperature. Circles denote MD calculations and lines denote continuum-scale models.1,103


044501-11

Michalis et al.

consistently deviates by 3.5 K form the experimental values for each given pressure. We observe that good agreement between the MD simulations and the two thermodynamic models is obtained. Such behavior should be expected for the TIP4P/Ice water model based on the findings by Docherty et al.105 and the discussion presented by Conde and Vega.76 In particular, Docherty et al. reported that in order to describe correctly the excess chemical potential of methane in water, a positive deviation from the energetic Lorentz–Berthelot combining rule needs to be introduced for the cross-interaction parameters of water and methane for the case of TIP4P/200587 model for water in combination with the Lennard-Jones parameters reported by Guillot and Guissani89 for methane. In that particular case, the deviation is 7%. As discussed by Conde and Vega,76 such an approach produces, essentially, the same cross-interaction energy parameter as in the case of TIP4P/Ice with Guillot and Guissani methane model using standard Lorentz–Berthelot rules. Additional discussion related to methane dissolved in water can be found in the work of Reed and Westacott.106

J. Chem. Phys. 142, 044501 (2015)

TIP4P/Ice water force field of the melting temperature of ice from the respective experimental value. Additionally, emphasis was given on the problems that arise from the inherently erroneous use of dispersion corrections for anisotropic and inhomogeneous systems in combination with pressure coupling. Finally, the solubility of methane in the liquid water was calculated along the three-phase coexistence line and good agreement was found with two different continuum scale models. ACKNOWLEDGMENTS

This publication was made possible by NPRP Grant No. 6-1547-2-632 from the Qatar National Research fund (a member of the Qatar Foundation). The statements made herein are solely the responsibility of the authors. We are grateful to the High Performance Computing Center of Texas A&M University at Qatar for generous resource allocation. 1E.

IV. CONCLUSIONS

The direct phase coexistence method was implemented in this work for the determination of the three-phase coexistence conditions of the methane hydrate. Four different pressures were examined (40, 100, 400, and 600 bar) that cover a wide range of the three phase equilibrium conditions and the respective coexistence temperatures T3 were calculated. The determination of T3 was based on statistical analysis of multiple independent runs per pressure, totaling up to 25 simulations per pressure with a time length of each simulation in the range of 1000–4000 ns. This averaging procedure was deemed necessary in order to compensate for the inherent stochasticity of the method. It was shown that the stochastic behavior of the system, if not treated with caution, can lead in the case of a small number of runs to erroneous results. It was additionally observed that the supersaturated conditions that can occur during this kind of simulations can affect the calculated T3. For this reason in all the simulations for the determination of T3 in this study, the instantaneous growth of methane bubble was avoided by using an adequate number of methane molecules. By avoiding methane supersaturation, the hydrate growth and dissociation rates are significantly lower. This fact leads to the need of relatively long simulations for the elucidation of the final state of the system (i.e., growth, dissociation). The aforementioned procedure resulted in the determination of T3 with the direct phase coexistence method with a high level of accuracy. It was found that for the TIP4P/Ice water model that was used in this study in combination with the OPLS-UA force field for methane, the MD results follow the same trend as the experimental values and they have excellent agreement with the expected T3 which is consistently 3.15 K below the experimental. This difference was expected in accordance with the work of Conde and Vega and it is attributed to the difference in the prediction by the

D. Sloan and C. A. Koh, Clathrate Hydrates of Natural Gases, 3rd ed. (CRC Press, Taylor and Francis Group, 2008). 2L. C. Jacobson, W. Hujo, and V. Molinero, J. Phys. Chem. B 113, 10298 (2009). 3A. Falenty, T. C. Hansen, and W. F. Kuhs, in Proceedings of the International Conference on Gas Hydrates (ICGH8), 28 July–1 August 2014, Beijing, China, Paper No. T1-123. 4M. M. Conde, C. Vega, G. A. Tribello, and B. Slater, J. Chem. Phys. 131, 034510 (2009). 5I. M. Chou, A. Sharma, R. C. Burruss, J. F. Shu, H. K. Mao, R. J. Hemley, A. F. Goncharov, L. A. Stern, and S. H. Kirby, Proc. Natl. Acad. Sci. U. S. A. 97, 13484 (2000). 6J. S. Loveday and R. J. Nelmes, Phys. Chem. Chem. Phys. 10, 937 (2008). 7E. D. Sloan, Am. Mineral. 89, 1155 (2004). 8E. G. Hammerschmidt, Ind. Eng. Chem. 26, 851 (1934). 9C. A. Koh, Chem. Soc. Rev. 31, 157 (2002). 10M. A. Kelland, Energy Fuels 20, 825 (2006). 11N. I. Papadimitriou, I. N. Tsimpanogiannis, and A. K. Stubos, Colloids Surf. A 357, 67 (2010). 12A. A. Khokhar, J. S. Gudmundsson, and E. D. Sloan, Fluid Phase Equilib. 150, 383 (1998). 13W. L. Mao, H. K. Mao, A. F. Goncharov, V. V. Struzhkin, Q. Z. Guo, J. Z. Hu, J. F. Shu, R. J. Hemley, M. Somayazulu, and Y. S. Zhao, Science 297, 2247 (2002). 14W. L. Mao and H. K. Mao, Proc. Natl. Acad. Sci. U. S. A. 101, 708 (2004). 15N. I. Papadimitriou, I. N. Tsimpanogiannis, A. T. Papaioannou, and A. K. Stubos, J. Phys. Chem. C 112, 10294 (2008). 16N. I. Papadimitriou, I. N. Tsimpanogiannis, C. J. Peters, A. T. Papaioannou, and A. K. Stubos, J. Phys. Chem. B 112, 14206 (2008). 17N. Dabrowski, C. Windmeier, and L. R. Oellrich, Energy Fuels 23, 5603 (2009). 18D. L. Zhong, N. Daraboina, and P. Englezos, Energy Fuels 106, 425 (2013). 19S. P. Kang and H. Lee, Environ. Sci. Technol. 34, 4397 (2000). 20Y. T. Seo, I. L. Moudrakovski, J. A. Ripmeester, J. W. Lee, and H. Lee, Environ. Sci. Technol. 39, 2315 (2005). 21A. Adeyemo, R. Kumar, P. Linga, J. Ripmeester, and P. Englezos, Int. J. Greenhouse Gas Control 4, 478 (2010). 22H. Kubota, K. Shimizu, Y. Tanaka, and T. Makita, J. Chem. Eng. Jpn. 17, 423 (1984). 23K. N. Park, S. Y. Hong, J. W. Lee, K. C. Kang, Y. C. Lee, M. G. Ha, and J. D. Lee, Desalination 274, 91 (2011). 24J. B. Klauda and S. I. Sandler, Energy Fuels 19, 459 (2005). 25K. A. Kvenvolden, Proc. Natl. Acad. Sci. U. S. A 96, 3420 (1999). 26D. Archer, B. Buffett, and V. Brovkin, Proc. Natl. Acad. Sci. U. S. A. 106, 20596 (2009). 27N. Kukowski, J. Greinert, and S. Henrys, Mar. Geol. 272, 141 (2010). 28M. Maslin, M. Owen, R. Betts, S. Day, T. Dunkley Jones, and A. Ridgwell, Philos. Trans. R. Soc., A 368, 2369 (2010). 29Y. F. Makogon, S. A. Holditch, and T. Y. Makogon, J. Pet. Sci. Eng. 56, 14 (2007).


044501-12 30R.

Michalis et al.

Boswell and T. S. Collett, Energy Environ. Sci. 4, 1206 (2011). Pinero, M. Marquardt, C. Hensen, M. Haeckel, and K. Wallmann, Biogeosciences 10, 959 (2013). 32A. V. Milkov, Earth-Sci. Rev. 66, 183 (2004). 33E. B. Burwicz, L. H. Rupke, and K. Wallmann, Geochim. Cosmochim. Acta 75, 4562 (2011). 34M. R. Walsh, S. H. Hancock, S. J. Wilson, S. L. Patil, G. J. Moridis, R. Boswell, T. S. Collett, C. A. Koh, and E. D. Sloan, Energy Econ. 31, 815 (2009). 35G. J. Moridis, T. S. Collett, R. Boswell, M. Kurihara, M. T. Reagan, C. Koh, and E. D. Sloan, SPE Reservoir Eval. Eng. 12, 745 (2009). 36M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids (Oxford Science Publications, 1989). 37J. H. van der Waals and J. C. Platteeuw, Adv. Chem. Phys. 2, 1 (1958). 38W. R. Parrish and J. M. Prausnit, Ind. Eng. Chem. Process Des. Dev. 11, 26 (1972). 39G. D. Holder, S. P. Zetts, and N. Pradhan, Rev. Chem. Eng. 5, 1 (1988). 40A. L. Ballard and E. D. Sloan, Fluid Phase Equilib. 194, 371 (2002). 41A. L. Ballard and E. D. Sloan, Fluid Phase Equilib. 218, 15 (2004). 42Y. Y. Bozhko, O. S. Subbotin, V. M. Fomin, V. R. Belosludov, and Y. Kawazoe, J. Eng. Thermophys. 23, 9 (2014). 43S. Dufal, A. Galindo, G. Jackson, and A. J. Haslam, Mol. Phys. 110, 1223 (2012). 44J. S. Tse, M. L. Klein, and I. R. McDonald, J. Phys. Chem. 87, 4198 (1983). 45F. Jimenez-Angeles and A. Firoozabadi, J. Phys. Chem. C 118, 11310 (2014). 46S. Sarupria and P. G. Debenedetti, J. Phys. Chem. Lett. 3, 2942 (2012). 47M. R. Walsh, G. T. Beckham, C. A. Koh, E. D. Sloan, D. T. Wu, and A. K. Sum, J. Phys. Chem. C 115, 21241 (2011). 48G. J. Guo, Y. G. Zhang, C. J. Liu, and K. H. Li, Phys. Chem. Chem. Phys. 13, 12048 (2011). 49M. R. Walsh, C. A. Koh, E. D. Sloan, A. K. Sum, and D. T. Wu, Science 326, 1095 (2009). 50C. Moon, R. W. Hawtin, and P. M. Rodger, Faraday Discuss. 136, 367 (2007). 51J. Vatamanu and P. G. Kusalik, J. Chem. Phys. 126, 124703 (2007). 52C. Moon, P. C. Taylor, and P. M. Rodger, J. Am. Chem. Soc. 125, 4706 (2003). 53T. Yagasaki, M. Matsumoto, Y. Andoh, S. Okazaki, and H. Tanaka, J. Phys. Chem. B 118, 1900 (2014). 54S. A. Bagherzadeh, P. Englezos, S. Alavi, and J. A. Ripmeester, J. Chem. Thermodyn. 44, 13 (2012). 55S. Alavi and J. A. Ripmeester, J. Chem. Phys. 132, 144703 (2010). 56E. M. Myshakin, H. Jiang, R. P. Warzinski, and K. D. Jordan, J. Phys. Chem. A 113, 1913 (2009). 57N. J. English, J. K. Johnson, and C. E. Taylor, J. Chem. Phys. 123, 244503 (2005). 58P. M. Rodger, Annals N.Y. Acad. Sci. 912, 474 (2000). 59B. J. Anderson, J. W. Tester, G. P. Borghi, and B. L. Trout, J. Am. Chem. Soc. 127, 17852 (2005). 60D. M. Duffy, C. Moon, and P. M. Rodger, Mol. Phys. 102, 203 (2004). 61C. Moon, P. C. Taylor, and P. M. Rodger, Can. J. Phys. 81, 451 (2003). 62V. S. Baghel, R. Kumar, and S. Roy, J. Phys. Chem. C 117, 12172 (2013). 63N. J. English and G. M. Phelan, J. Chem. Phys. 131, 074704 (2009). 64N. J. English and J. S. Tse, Phys. Rev. Lett. 103, 015901 (2009). 65A. Phan, D. R. Cole, and A. Striolo, J. Phys. Chem. C 118, 4860 (2014). 66S. A. Bagherzadeh, P. Englezos, S. Alavi, and J. A. Ripmeester, J. Phys. Chem. B 116, 3188 (2012). 67S. Liang, D. Rozmanov, and P. G. Kusalik, Phys. Chem. Chem. Phys. 13, 19856 (2011). 68R. T. Cygan, S. Guggenheim, and A. F. K. van Groos, J. Phys. Chem. B 108, 15141 (2004). 69Y. Iwai, H. Nakamura, and M. Hirata, Mol. Simul. 38, 481 (2012). 31E.

J. Chem. Phys. 142, 044501 (2015) 70Y.

Qi, M. Ota, and H. Zhang, Energy Convers. Manage. 52, 2682 (2011). Y. Geng, H. Wen, and H. Zhou, J. Phys. Chem. A 113, 5463 (2009). 72A. Svandal, T. Kuznetsova, and B. Kvamme, Fluid Phase Equilib. 246, 177 (2006). 73H. Jiang, E. M. Myshakin, K. D. Jordan, and R. P. Warzinski, J. Phys. Chem. B 112, 10207 (2008). 74N. J. English, Mol. Phys. 106, 1887 (2008). 75E. J. Rosenbaum, N. J. English, J. K. Johnson, D. W. Shaw, and R. P. Warzinski, J. Phys. Chem. B 111, 13194 (2007). 76M. M. Conde and C. Vega, J. Chem. Phys. 133, 064507 (2010). 77Y. T. Tung, L. J. Chen, Y. P. Chen, and S. T. Lin, J. Phys. Chem. B 114, 10804 (2010). 78A. R. Finney and P. M. Rodger, Phys. Chem. Chem. Phys. 13, 19979 (2011). 79G. S. Smirnov and V. V. Stegailov, J. Chem. Phys. 136, 044523 (2012). 80M. M. Conde and C. Vega, J. Chem. Phys. 138, 056101 (2013). 81J. L. F. Abascal, E. Sanz, R. García Fernández, and C. Vega, J. Chem. Phys. 122, 234511 (2005). 82H. W. Horn, W. C. Swope, J. W. Pitera, J. D. Madura, T. J. Dick, G. L. Hura, and T. Head-Gordon, J. Chem. Phys. 120, 9665 (2004). 83W. L. Jorgensen, D. S. Maxwell, and J. Tirado-Rives, J. Am. Chem. Soc. 118, 11225 (1996). 84Z. T. Cao, J. W. Tester, K. A. Sparks, and B. L. Trout, J. Phys. Chem. B 105, 10950 (2001). 85L. Jensen, K. Thomsen, N. von Solms, S. Wierzchowski, M. R. Walsh, C. A. Koh, E. D. Sloan, D. T. Wu, and A. K. Sum, J. Phys. Chem. B 114, 5775 (2010). 86W. L. Jorgensen, J. D. Madura, and C. J. Swenson, J. Am. Chem. Soc. 106, 6638 (1984). 87J. L. F. Abascal and C. Vega, J. Chem. Phys. 123, 234505 (2005). 88H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, J. Phys. Chem. 91, 6269 (1987). 89B. Guillot and Y. Guissani, J. Chem. Phys. 99, 8075 (1993). 90A. J. C. Ladd and L. V. Woodcock, Chem. Phys. Lett. 51, 155 (1977). 91R. K. McMullan and G. A. Jeffrey, J. Chem. Phys. 42, 2725 (1965). 92J. D. Bernal and R. H. Fowler, J. Chem. Phys. 1, 515 (1933). 93S. Sarupria and P. G. Debenedetti, J. Phys. Chem. A 115, 6102 (2011). 94R. G. Fernandez, J. L. F. Abascal, and C. Vega, J. Chem. Phys. 124, 144506 (2006). 95D. Van der Spoel, E. Lindahl, B. Hess, G. Groenhof, A. E. Mark, and H. J. C. Berendsen, J. Comput. Chem. 26, 1701 (2005). 96B. Hess, C. Kutzner, D. van der Spoel, and E. Lindahl, J. Chem. Theory Comput. 4, 435 (2008). 97S. Pronk, S. Pall, R. Schulz, P. Larsson, P. Bjelkmar, R. Apostolov, M. R. Shirts, J. C. Smith, P. M. Kasson, D. van der Spoel, B. Hess, and E. Lindahl, Bioinformatics 29, 845 (2013). 98H. J. C. Berendsen, J. P. M. Postma, W. F. Vangunsteren, A. Dinola, and J. R. Haak, J. Chem. Phys. 81, 3684 (1984). 99U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee, and L. G. Pedersen, J. Chem. Phys. 103, 8577 (1995). 100G. J. Moridis, SPE J. 8, 359 (2003). 101E. W. Lemmon, M. O. McLinden, and D. G. Friend, “Thermophysical properties of fluid systems,” NIST Chemistry WebBook, NIST Standard Reference Database Number 69, edited by P. J. Linstrom and W. G. Mallard (National Institute of Standards and Technology, Gaithersburg MD, 2014). 102J. R. Espinosa, E. Sanz, C. Valeriani, and C. Vega, J. Chem. Phys. 139, 144502 (2013). 103Z. H. Duan and S. D. Mao, Geochim. Cosmochim. Acta 70, 3369 (2006). 104I. N. Tsimpanogiannis, I. G. Economou, and A. K. Stubos, Fluid Phase Equilib. 371, 106 (2014). 105H. Docherty, A. Galindo, C. Vega, and E. Sanz, J. Chem. Phys. 125, 074510 (2006). 106S. K. Reed and R. E. Westacott, Phys. Chem. Chem. Phys. 10, 4614 (2008). 71C.


Phase Coexistence in a Dynamic Phase Diagram.

A Small-Volume Apparatus for the Measurement of Phase Equilibria.

Phase Equilibria and Crystallography of Ceramic Oxides.

Phase equilibria in the neodymium-cadmium binary system.

Azithromycin hydrates-implications of processing-induced phase transformations.

Phase equilibria and stability characteristics of chlorpropamide-urea solid dispersions.

Dynamic phase coexistence in glass-forming liquids.

Prediction of prostate cancer recurrence using quantitative phase imaging.

Fermi-surface reconstruction and complex phase equilibria in CaFe2As2.

Pressure and phase equilibria in interacting active brownian spheres.

Mechanical instability of monocrystalline and polycrystalline methane hydrates.

Development of the Transferable Potentials for Phase Equilibria Model for Hydrogen Sulfide.

Mobile phase versus stationary phase approaches to the direct injection of biological fluids in liquid chromatography.

A generalized phase resetting method for phase-locked modes prediction.

Experimental verification of methane-carbon dioxide replacement in natural gas hydrates using a differential scanning calorimeter.

Ld phase coexistence modulation induced by GM1.

A possible four-phase coexistence in a single-component system.

Methodology of isfahan tobacco use prevention program: first phase.

Phase coexistence in single-lipid membranes induced by buffering agents.

On the determination of phase boundaries via thermodynamic integration across coexistence regions.

Polyaromatic hydrocarbons do not disturb liquid-liquid phase coexistence, but increase the fluidity of model membranes.

Inhibited phase behavior of gas hydrates in graphene oxide: influences of surface and geometric constraints.

Time-resolved X-ray diffraction and Raman studies of the phase transition mechanisms of methane hydrate.

Evaluation of the Phase-Dependent Rhythm Control of Human Walking Using Phase Response Curves.