On the Thermodynamic Stability of Clathrate Hydrates VI: Complete

Subscriber access provided by Duke University Libraries

Article

On the Thermodynamic Stability of Clathrate Hydrates VI: Complete Phase Diagram Hideki Tanaka, Takuma Yagasaki, and Masakazu Matsumoto J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.7b10581 • Publication Date (Web): 06 Dec 2017 Downloaded from http://pubs.acs.org on December 12, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry B is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

On the Thermodynamic Stability of Clathrate Hydrates VI: Complete Phase Diagram

Hideki Tanaka*, Takuma Yagasaki, and Masakazu Matsumoto Research Institute for Interdisciplinary Science, Okayama University, Okayama, 700-8530, Japan

ABSTRACT We develop a method to evaluate the thermodynamic stability of clathrate hydrates relative to host water and/or guest species. This enables to investigate complete phase behaviors of clathrate hydrates in the whole space of the thermodynamic variables, not only temperature and pressure but also composition, with only the intermolecular interactions as input parameters. A complete phase diagram of clathrate hydrate is settled with this method, specifically the region enclosed by the hydrate/water and hydrate/guest phase boundaries where a clathrate hydrate is the only stable phase. The method is applied to methane clathrate hydrate, which results in an excellent agreement in dissociation pressure with the experimental observations. It is found that the hydrate/water phase boundary is significantly affected by the phase transition of water from ice to liquid. This transition limits the stable area of the clathrate hydrate terminated at the dissociation temperature, which otherwise exhibits unphysical divergence. It is essential to choose a pressure as an independent variable so as to calculate the accurate phase equilibria in composition space. The present method 1 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 42

establishes a pressure-temperature-composition relation for a single stable phase of clathrate hydrate as well as recovering its temperature dependence of the dissociation pressure.

I. INTRODUCTION Clathrate hydrate is a nonstoichiometric compound composed of host water and small guest molecules, the latter of which are mostly hydrophobic or only weakly polarized. Because guest molecules are contained in the cages, which are usually 12, 14, 15, and/or 16-hedra made of water molecules, the size of guest molecules is bounded to butane or neo-pentane.1-2 Methane clathrate hydrate in offshore and permafrost regions is expected to be an energy resource.2 Methane clathrate hydrate can also be used for storage of methane. It is essential to elucidate both the thermodynamic stability and the phase transition mechanism of clathrate hydrates from a microscopic view point for their efficient usage. This material also serves to study the hydrophobic hydration in which the structure around a solute is believed to have some resemblance to the clathrate hydrate structure yet its extremely low solubility in liquid water (10-4 to 10-5 in mole fraction) contrasts the high occupancy in clathrate hydrate (more than 0.1).3-5 Computer simulations have provided a wealth of information on properties of clathrate hydrates at the molecular level. Molecular dynamics (MD) simulations have been carried out to investigate not only time dependent processes but also thermodynamic and structural properties of clathrate hydrates.6-10

Various static

properties, such as the pressure dependence of the guest occupancy, have been evaluated by grand canonical Monte Carlo (MC) simulations.11-15 An equilibrium condition for a clathrate hydrate to be formed or dissociated has been 2 ACS Paragon Plus Environment

Page 3 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


conducted on the basis of the van der Waals and Platteeuw (vdWP) theory.16

This

theory enables to draw a phase diagram in temperature and pressure space.16 It relies on a semi-open system of the clathrate hydrate in equilibrium with both an aqueous (or ice) phase and a guest-rich phase. In such a system, the chemical potential of water in the clathrate hydrate is the same as that in the aqueous phase and the chemical potential of the guest species in the clathrate hydrate equals to that in the guest-rich phase. The temperature dependence of the dissociation pressure and the cage occupancies at dissociation pressures for various clathrate hydrates can be predicted with some empirical parameters such as the chemical potential difference between water (ice) and hypothetical empty clathrate hydrate. There are several restrictions in the original vdWP theory: (1) neglect of influence of a large guest species on the host lattice, (2) neglect of guest-guest interactions, (3) restriction of the cage occupation by at most one guest molecule, and (4) constraint to a constant volume even at a high pressure. We have developed methodologies to draw the phase diagram of clathrate hydrates computationally based on the vdWP theory.17-23

In these studies, the empirical

parameters in applying the vdWP theory are replaced by those obtained from molecular mechanical calculations using intermolecular interaction models, which are capable of reproducing the phase diagrams of pure water (ice) and the guest species. The restrictions mentioned above have been removed in the series of our studies (except for (2)) to deal with various types of clathrate hydrates and to yield a better agreement with experimental observations by (1) allowing the shift of vibrational frequency in the host lattice due to the guest molecules,17-19 (2) introducing the occupancy-dependent free energy calculated from a mean field treatment,22-23 (3) extending the theory to multiple-occupation,20 and (4) conversion to a constant pressure system.21 3 ACS Paragon Plus Environment


Clathrate hydrates are multi-component systems and thus compositions should be considered as thermodynamic variables in addition to temperature and pressure. Here, we consider a phase diagram of a single guest component clathrate hydrate in temperature-composition space at a fixed pressure. The Gibbs phase rule claims that the boundary at a given pressure is represented by a line in the phase diagram for a two-component and two-phase system. A typical result of our new theory is demonstrated in Figure 1, in which two phase boundaries are drawn for clathrate hydrate; one is the hydrate/water boundary (left blue line in Figure 1) and the other is the hydrate/guest boundary (right blue line). There must be a finite region enclosed by the two boundaries, no matter how narrow it is. In this region, the clathrate hydrate is the only stable phase. For a two-component and single-phase system, the cage occupancy varies with the composition according to the Gibbs phase rule. The range of the composition becomes smaller with increasing temperature and disappears at a certain temperature, called a dissociation temperature (pressure), where the three phases (hydrate, water, and gas) coexist. The original vdWP theory can estimate only the dissociation temperature (pressure) because of the assumption of the semi-open system.2 However, clathrate hydrates can form in water-rich and guest-rich conditions which are somewhere away from the dissociation temperature (pressure) in the phase diagram. Thus, it is of particular interest to consider the phase boundary of clathrate hydrate which is in equilibrium with either the aqueous (ice) or guest fluid phase. Moreover, information on guest-contents at various temperatures and pressures is indispensable in industrial use of clathrate hydrate as an energy resource or a storage medium.

4 ACS Paragon Plus Environment

Page 4 of 42

Page 5 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 1. Schematic temperature-composition diagram of methane clathrate hydrate at 10 MPa. The vertical line (red dotted) is the line boundary of methane clathrate hydrate assuming the full occupancy and the hatched region enclosed by the blue lines is the single stable phase of methane clathrate hydrate conducted in the present study.

Although methane clathrate hydrate is a nonstoichiometric compound, its phase behavior in temperature-composition space was first proposed assuming the full occupation of the available cages regardless of temperature and pressure.24 The region where methane clathrate hydrate is the only stable phase has long been represented as a vertical line in the temperature-composition phase diagram as shown in Figure 1. Later, a region for methane clathrate hydrate was estimated to be unimodal on the basis of experiments.25-26 A theoretical calculation was carried out to estimate the region where a clathrate hydrate is the only stable phase using a water interaction model which is not so sophisticated but has an essential element of water to form four hydrogen bonds.27 The



Page 6 of 42

key issue of this study is the application of the Gibbs-Duhem equation to calculate the chemical potentials of water at various filling ratios, which were obtained by semi-grandcanonical MC simulations.28-29 In response to this pioneering work, it is desirable to construct a simple theoretical method to evaluate the phase behaviors of clathrate hydrates using more accurate intermolecular interaction parameters of water and guest species for drawing a qualitatively correct phase diagram in a wide range of pressure-temperature-composition space. It is impossible to draw the phase diagram of clathrate hydrate in composition space by either the original vdWP theory16 or the methodologies proposed in our previous studies.17-23 In the present study, we develop a statistical mechanical theory to predict phase behaviors in temperature-pressure-composition space. This is performed by conversion from a semi-grand canonical to a canonical system, thereby changing the independent variable from the chemical potential of guest species to the composition along with that from the volume to the pressure.

Heavy MD or MC simulations are

not required in the new method. We estimate the location of the hydrate/water and hydrate/guest phase boundaries of methane clathrate hydrate up to 200 MPa where hexagonal ice (ice Ih) or liquid water can be in equilibrium with clathrate hydrate of either structure I or II (CS-I or CS-II).2, 30 The shape of the region enclosed by the two boundaries is found to be not a simple unimodal but significantly dependent on the pressure. It is also found that the area of the single stable phase of methane clathrate hydrate becomes smaller in the temperature-composition space as applying the higher pressure and that the magnitude of the area cannot be estimated correctly without allowing the unit cell dimension to respond to the applied pressure.

The temperature

and pressure dependences of the cage occupancy and the chemical potential of methane 6 ACS Paragon Plus Environment

Page 7 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


in clathrate hydrate are also discussed. The paper is organized as follows. A theoretical method is presented to treat the composition variation as well as the pressure variation in reference to the original vdWP theory and its extensions in Sec. II. A method of practical calculation is given in Sec. III. We show its application to methane clathrate hydrate to draw the phase diagram in Sec. IV. Our study is concluded in Sec. V with a few remarks.

II. THEORY A. van der Waals and Platteeuw Theory for Semi-open System Let us begin with the original vdWP theory1-2,

16

which provides a basis of the

revised theories. A clathrate hydrate is made of Nw water molecules at temperature T (or

β; the inverse of T times the Boltzmann constant, k B ) confined in a volume V. A guest molecule is accommodated in a cage of type j, where j runs over the possible types. We consider a single component of guest species and a fixed number of guest molecules accommodated in each cage. An extension to multiple-component of guest species is straightforward and is given elsewhere.2,

16, 21

According to the assumptions in the

original vdWP theory, the canonical partition function, Z , for a clathrate hydrate accommodating n j guest molecules out of its N j cages of type j is written by the product of the host and guest parts as Nj  Z = exp( − β Ac0 ) × ∏   exp( − β n j f j )   j n j 

(1)

where Ac0 stands for the Helmholtz free energy of the hypothetical empty clathrate hydrate and f j indicates the free energy of the cage occupation by the guest species



Page 8 of 42

(note that CS-II empty clathrate hydrate turned out to be a metastable but real entity in 2014).31 In deriving eq 1, we assume that occupation of one cage is independent of occupation of other cages, which is conducted from neglect of the guest-guest interactions and random distribution of guests over the available cages, each corresponding to energetic and entropic contribution. The canonical partition function is converted to the grandcanonical one, Ξ , by multiplying exp( β n j µ g ) and carrying out the summation over the possible occupations. This leads to N j   exp( − β n j f j ) Ξ = exp( − βA ) × ∏ ∑ exp( βn j µ g ) n  j n j =0  j  , Nj 0 = exp( − β Ac ) × ∏ (exp[ β ( µ g − f j )] + 1)

(2)

where µg is the chemical potential of the guest species in the guest-rich phase.

The

Nj

0 c

j

above formulation implies that the clathrate hydrate thus described is an open system with respect to the guest species. Therefore, the mean occupancy is derived from the partition function as

< x >=

exp[ β ( µ g − f j )] ∂ ln Ξ = ∑α j = ∑α j x j N w ∂ ( βµ g ) 1 + exp[β ( µ g − f j )] j j

(3)

where αj =Nj / Nw is the ratio of the number of j-type cages to that of water and x j is its occupancy (the filling ratio) by the guest. The chemical potential of water is calculated from the thermodynamic potential of the clathrate hydrate, Φ(= −k B T ln Ξ) and is expressed by the two contributions; (1) the chemical potential of water in the empty clathrate hydrate, µ c0 =

∂Ac0 , which is regarded as a reference state, and (2) the term ∂N w

associated with the cage occupation. It is written as


Page 9 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


µc =

∂Φ = µ c0 − k BT ∑ α j ln{1 + exp[ β ( µ g − f j )]} ∂N w j

 exp[ β ( µ g − f j )]   = µ c0 + k BT ∑ α j ln 1 −   1 + exp[ ( − f )] β µ j g j   0 = µ c + k BT ∑ α j ln(1 − x j )

.

(4)

j

B. Stability of Closed Systems with Fixed Number of Guest Molecules In order to calculate a complete phase diagram of a clathrate hydrate, we must describe the phase behaviors in pressure-temperature-composition space. In order for the composition variable to be implemented in predicting the phase behaviors, we consider three types of phase equilibria as i) Two phase coexistence of aqueous solution (ice) and clathrate hydrate, ii) Two phase coexistence of clathrate hydrate and guest fluid, iii) Single stable phase of clathrate hydrate. For simplicity, water-rich and guest-rich phases in equilibrium with a clathrate hydrate are assumed to be pure liquid water (or ice) and pure fluid methane, respectively. This is justified by the fact that the mutual solubilities are extremely low.3 We examine thermodynamic stability of clathrate hydrates with the fixed number of guest molecules, which enables to determine a region of a clathrate hydrate being the only stable phase. We choose an appropriate thermodynamic potential at a given temperature and a pressure, focusing particularly on the single stable phase. It is evident that the stability under the above circumstances is estimated by the difference in the Gibbs free energy between the product (clathrate hydrate) and the reactant components (water and guest species) denoted by ∆G , i.e., the Gibbs free energy relative to the pure substances having the chemical potential of water, µw, and that of guest, µg, assuming



Page 10 of 42

that the clathrate hydrate is in equilibrium with either pure water or pure guest species as mentioned above. Under the assumption in deriving eq 1, the Gibbs free energy difference is given by Nj   ∆G = N w ( µ c0 − µ w ) + ∑ n j ( f j − µ g ) − k B T ∑ ln  n  j j j   0 ≈ N w (µ c − µ w ) + ∑ n j ( f j − µ g )

(5)

j

− k B T ∑ [ N j ln N j − n j ln n j − ( N j − n j ) ln( N j − n j )] j

where Stirling’s approximation is applied. Since we can specify only the total number of guest molecules, the number of guest molecules in each cage type is subject to the constraint as n = ∑ n j . The occupancy of j-type cage is given by j

xj =

exp[ β ( µ − f j )]

(6)

exp[ β ( µ − f j )] + 1

with the Lagrange multiplier, exp(βµ ) , where µ corresponds to the chemical potential of the guest species in the clathrate hydrate (µ should be distinguished from the chemical potential of the guest species in the guest-rich phase, µg). The Gibbs free energy difference,

[

]

∆ G / N w = ( µ c0 − µ w ) + ∑ α j k B T x j ln x j + (1 − x j ) ln(1 − x j ) + β x j ( f j − µ g ) ,

(7)

j

should satisfy the following inequalities relative to the pure components for the clathrate hydrate to be the most stable phase. Since only Nj is proportional to Nw, we have

∂∆G = ∆µ w = ( µ c0 − µ w ) + k BT ∑ α j ln(1 − x j ) < 0 . ∂N w j

(8)

Since nj is irrelevant to Nw (its maximum value is given a priori), the inequality holds for any j as


Page 11 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


∂∆G = ∆µ j = ( f j − µ g ) + k BT [ln x j − ln(1 − x j )] = ( µ − µ g ) < 0 . ∂n j

(9)

C. Conversion of Independent Variable from Volume to Pressure In order to convert the independent variable of solids from the volume to the pressure, the following standard procedure is adopted. That is, for a single component solid such as ice, the Helmholtz free energy, A, is a function of T and V and the mean volume, , is calculated by minimizing the function,

Ψ (T , p;V ) = A(T ,V ) + pV ,

(10)

with respect to the volume at a given pressure, p. The Gibbs free energy is calculated as G(T , p) = A(T , V ) + p V

(11)

Similarly, the Gibbs free energy of the clathrate hydrate, Gc (T , p, y ) , is obtained by minimizing the following function

[

]

Ψ (T , p, y;V ) = Ac0 + pV + k B TN w ∑ α j β x j f j + x j ln x j + (1 − x j ) ln(1 − x j ) .

(12)

j

where y is the mole fraction of the guest species. It is rewritten as Ψ (T , p, y;V ) = Ac0 + pV + nµ + k B TN w ∑ α j ln(1 − x j ) .

(13)

j

The relation between the mole fraction, y, and the occupancy of each cage type, xj, is established as

y = n /(n + Nw ) = ∑n j /(∑ n j + Nw )

(14)

n j = α j x j Nw

(15)

j

j

For a usual clathrate hydrate containing two types of cages, large and small, denoted by subscripts l and s, the chemical potential of guest species is calculated with

C = exp(βµ)

(16)



Page 12 of 42

as

[(α l + α s )(1 − y ) − y ]exp[ − β ( f l + f s )]C 2 + [(1 − y ){α l exp( − β f l ) + α s exp( − β f s )} − {exp( − β f l ) + exp( − β f s )} y ]C − y = 0

.

(17)

From the above thermodynamic potential, we obtain the chemical potentials of water and guest as

∂Gc = µc = µc0 + k BT ∑α j ln(1 − x j ) , ∂N w j ∂Gc ∂µ ∂µ = µ+n − k BTN w ∑ α j x j β =µ ∂n ∂n ∂n j

(18) .

(19)

D. Stability Conditions and Phase Boundaries Hereafter, the quantities such as A, µ, and fj denote those values at the mean volume except for the minimizing process of the corresponding function. The free energy difference is given by

[

]

∆G = Ac0 + pV − N w µ w + ∑ α j N w x j ( µ − µ g ) − k BT ln[1 + exp{β ( µ − f j )}] .

(20)

j

We consider three cases below to draw a pressure-temperature-composition diagram as mentioned in subsection B. i) For two-phase equilibrium between liquid water (or ice) and clathrate hydrate, the chemical potential of water in the pure liquid water (ice) phase should be the same value as  ( Ac0 + pV )   − k B T ∑ α j ln[1 + exp{ β ( µ − f j )}] Nw j   .

µ w (T , p ) = 

= µ − k B T ∑ α j ln[1 + exp{β ( µ − f j )}]

(21)

0 c

j

Note that the right hand side of eq 21, the chemical potential of water in clathrate hydrate, can also be evaluated from a different method, semi-grandcanonical MC


Page 13 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


simulations.19,28-29,

32

The equilibrium volume at a given pressure is calculated by

minimizing the corresponding function as ∂Ac0 ∂ + p+ ∂V ∂V

∑α

j

[

]

N w x j µ − k BT ln[1 + exp{β ( µ − f j )}] = 0 .

(22)

j

ii) Two phase equilibrium between clathrate hydrate and guest species occurs under the following condition,

x j = exp[β ( µ g − f j )] /{exp[β ( µ g − f j )] + 1} .

(23)

Since the chemical potential of the guest species, µg, is independent of the volume, the volume at a given pressure is obtained by ∂Ac0 ∂ + p− ∂V ∂V

∑α

j

N w k B T ln[1 + exp{β ( µ g − f j )}] = 0

(24)

j

iii) The following inequalities should hold for a single stable phase of clathrate hydrate relative to water and guest species,  Ac0 + pV   Nw

  − k BT ∑ α j ln[1 + exp{β ( µ − f j )}] < µ w and µ < µ g j 

(25)

Each of the above conditions is also derived from an equilibrium of coexisting two phases.33 The equilibrium volume of clathrate hydrate with variable composition is obtained in the same way as in the case of the clathrate hydrate coexisting with water, eq 22. The chemical potential of the guest species, µ, in the clathrate hydrate is connected smoothly to that of pure fluid, µg . It is also noted that the volume of the clathrate hydrate of the single stable phase calculated according to eq 22 is connected to that at the coexistence of the two phases described by eq 24 considering the fact that µg is pertinent solely to the guest fluid and does not depend on V.



III. COMPUTATIONAL METHOD A. Intermolecular Interaction and Generation of Ice and Clathrate Hydrate Structures Clathrate hydrates examined in the present study are composed of water and guest molecules.

The water-water interaction is described by the TIP4P/ice,34 which is

advantageous to reproduce the phase diagram of water, especially at low pressures. Methane is chosen as a guest species, which is approximated by a spherical molecule represented by a single Lennard-Jones (LJ) interaction site. The size and energy parameters are set to the values either from a classical textbook35 or the OPLS model36 and as tabulated in Table I. The Lorentz-Berthelot rule is applied to any unlike interaction. For calculation of the free energies, crystalline structures of clathrate hydrates CS-I and CS-II, and ice Ih are required. We generate 100 hydrogen disordered configurations for each structure containing 384 water molecules for CS-I, 1,088 for CS-II, and 1,120 for ice Ih using the GenIce tool.37 Each generated configuration has vanishing net polarization and satisfies the ice rule where an individual water molecule has four hydrogen bonded neighbors.

Table I Lennard-Jones size (σ) and energy (ε) parameters for methane. Type

σ / nm

ε / kJ mol-1

A1)

0.3882

1.139

B2)

0.3817

1.232

1) reference 35, 2) reference 36.


Page 14 of 42

Page 15 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


B. Free Energy and Chemical Potential of Host Water The Helmholtz free energy of ice Ih and empty clathrate hydrate at temperature T is assumed to be the sum of the three terms; the interaction energy at T = 0 K, i.e., in the inherent structure,38 the free energy of the intermolecular vibrations, and the residual entropy as

A(T , V ) = Uq (V ) + Fv (T ,V ) − TSr .

(26)

The first term Uq is obtained by averaging the interaction energies over the 100 hydrogen-disordered configurations at a given volume V. The free energy of the intermolecular vibrational motions, Fv, is divided into the harmonic and anharmonic contributions as, Fv (T ,V ) =F h(T ,V ) + Fa (T ,V ) .

(27)

The classical harmonic vibrational free energy, Fhc , is expressed as  hν Fhc (T , V ) = k BT ∑ ln j j  k BT

  

(28)

and quantum one, Fhq , as

  hν j   , Fhq (T ,V ) = kBT ∑ ln2 sinh j  2k BT  

(29)

where h stands for the Planck constant and νi is the frequency of i-th normal mode of the lattice.39 These lattice dynamics equations are combined with the minimization of the function given as eq 10 to obtain the Gibbs free energy and the equilibrium volume. This is called a quasi-harmonic approximation and was applied to ice Ih and cubic (Ic) ice.40 The residual entropy in eq 26, Sr, is given by kB Nw ln(3/2) but this entropy is common to ice and clathrate hydrates considered in the present study.41 The chemical potential of ice, µ i (T , p ) , is decomposed into two terms, the chemical



Page 16 of 42

potential under the quasi-harmonic approximation of the intermolecular vibrations denoted as µ ih (T , p ) and the remaining term as ∆µ ia (T , p ) . The former property is calculated according to eq 26 by ignoring the anharmonic contribution in eq 27 as,

µ ih (T , p ) = [U q (< V > ) + Fh (T , < V > ) − TS r + p < V > ] / N w ,

(30)

and the latter is associated with the anharmonic vibrational motions. Similarly, the chemical potential of the empty clathrate hydrate µ c0 (T , p ) is divided into µ c0 h (T , p ) and ∆µ ca (T , p ) . An equilibrium of clathrate hydrate with ice described by eq 21 satisfies the following condition,

µ ih (T , p ) + ∆µ ia (T , p ) = µ c0 h (T , p ) + ∆µ ca (T , p ) − k B T ∑ α j ln[1 + exp{β ( µ − f j )}] (31) j

The anharmonic vibrational free energies may not be negligible. However, we can expect that the magnitudes of these properties are hardly different from each other, ∆µia (T , p ) ≈ ∆µ ca (T , p ) . This enables to subtract the anharmonic contributions from the

both sides of eq 31 in exploring the stability of clathrate hydrate relative to ice, i.e.,

µih (T , p ) = µ c0 h (T , p ) − k BT ∑ α j ln[1 + exp{β ( µ − f j )}]

(32)

j

For temperatures higher than the melting point of ice, Tm, the equilibrium of clathrate hydrate with liquid water is described by

µ l (T , p ) = µ lh (T , p ) + ∆µ ia (Tm , p ) = µ c0 h (T , p ) + ∆µ ca (T , p ) − k B T ∑ α j ln[1 + exp{ β ( µ − f j )}]

(33)

j

where µ l (T , p ) is the chemical potential of liquid water and µ lh (T , p ) does not mean the harmonic contribution in liquid water but is rather defined by

µ lh (T , p ) = µ l (T , p ) − ∆µ ia (Tm , p ) .

(34)


Page 17 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The melting point of the TIP4P/ice model at a given pressure is calculated from the reported values of Tm and dp/dT at 1 bar.34

The chemical potential of liquid water,

µl(T), at temperature T=Tm+∆T can be estimated roughly as

µ l (T , p ) ≈ µ i (Tm , p ) − ∆s m ∆T − ∆C p

∆T 2 2

(35)

from a set of the thermodynamic quantities, the chemical potential of ice, µi(Tm,p), the entropy difference between ice and liquid water, ∆ s m , and the heat capacity difference at constant pressure, ∆Cp. The heat of fusion, ∆H m , leading to ∆sm ( = ∆H m / Tm ) and the heat capacity difference are simply calculated from the enthalpies of both ice Ih and liquid water. Using eq 34, we have

µ lh (T , p) ≈ µ ih (Tm , p) − ∆s m ∆T − ∆C p

∆T 2 . 2

(36)

Assuming ∆µ ca (T , p ) ≈ ∆µ ia (Tm , p ) , we can remove again the anharmonic contributions appearing in eq 33. Then, the equilibrium is attained by the condition as

µ lh (T , p ) = µ c0 h (T , p ) − k BT ∑ α j ln[1 + exp{ β ( µ − f j )}] ,

(37)

j

where both sides can be calculated with the harmonic vibrations of crystalline solids. This may lead to a little less accurate phase behavior as being away from Tm compared with that between ice and clathrate hydrate but we still anticipate that the temperature

 ∂∆µ ca dependence of the contribution from anharmonic vibration,   ∂T

  ∆T , is not so  T =Tm

large compared with the terms in eq 37 considering the fact that large amplitude motions of host water are suppressed by the presence of the guest molecules.

C. Free Energy of Cage Occupation



Page 18 of 42

The free energy of cage occupation for a spherical molecule, f, is calculated according to classical mechanics as f c = − k B T ln[λ −3 ∫ exp{− βψ (r )}dr ]

(38)

where ψ(r) stands for the potential of the guest due to surrounding water molecules at a position r and λ is the de Broglie thermal wave length.

The guest-guest interaction is

included by filling all the cages for simplicity and thus a little overestimated. But this contribution to the total free energy is much less significant. The thermodynamic properties of the clathrate hydrate at a low temperature should be calculated based on quantum mechanics.

The quantum mechanical free energy is given by

f q = − k BT ln ∑ exp( − E j / k BT )

(39)

j

where Ei is i-th energy level of the guest molecule in the cage. The energy levels are obtained by solving numerically the radial Schrödinger equation assuming that the potential of the guest in a cage is spherically isotropic. The detailed method is given elsewhere.42

D. Chemical Potential of Methane Fluid The chemical potential of the guest plays a crucial role in determining the occupancies of individual cages. The ideal gas approximation is not accurate enough to evaluate this quantity. The critical temperature of methane, 190.56 K is just below the lowest bound temperature we are interested in ( T ≥ 200 K) and pure methane is a fluid phase under this condition.43 We calculate the chemical potential of guest methane with Redlich-Kwong equation of state44 which is advantageous to lead directly to the free energy of the fluid phase from the critical point of methane, Tc = 190.56 K and pc =


Page 19 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


4.599 MPa.43

IV. RESULTS AND DISCUSSION A. Free Energies of Ice and Empty and “Effective Reference” Clathrate Hydrates The Helmholtz free energy of ice Ih or empty clathrate hydrate is calculated according to eq 26. The Gibbs free energy of ice Ih, eq 11, with the harmonic vibrational free energy calculated from classical and quantum equations, eq 28 and 29, at 0.1 MPa is plotted in Figure 2 against temperature. Although the free energy value according to classical mechanics is fairly different from that according to quantum mechanics, the difference arises mainly from the zero-point energy. More importantly, the slope against temperature is negative at temperatures above 200 K, which is always satisfied in the quantum free energy. The molar volumes, Vm, of ice Ih by both the quantum and classical quasi-harmonic approximations at 0.1 MPa are also plotted in Figure 2 against temperature. The molar volumes at 0 K from the quasi-harmonic approximation with the classical and the quantum vibrational free energy are different from each other but the thermal expansivities are similar except for the low temperature region where ice Ih shrinks upon heating.40, 42, 45 The molar volume calculated from Ac+pV at 250 K and 0.1 MPa of 19.94 cm3 mol-1, is close to that obtained from a

classical MD simulation, 19.80 cm3 mol-1,34 and agrees fairly with the experimental value, 19.55 cm3 mol-1.43, 46



Figure 2. Chemical potential (black, left) and molar volume (red, right) of ice Ih against temperature at 0.1 MPa with the classical (solid line) or quantum (dotted line) harmonic free energy.

The free energies of the empty CS-I and CS-II clathrate hydrates are calculated in the same way as that of ice Ih. Although the Gibbs free energy of the empty clathrate hydrate at a fixed pressure alone is of less significance in the practical equilibrium, it contributes to the most part of the whole free energy of clathrate hydrate. The free energies for the empty CS-I and CS-II clathrate hydrates at 0.1 MPa are plotted against temperature in Figure 3. Either of them is higher than that for ice Ih by about 1 kJ mol-1. The free energy of the empty CS-II clathrate hydrate is close to that of the CS-I but is always lower by 0.1 to 0.2 kJ mol-1.


Page 20 of 42

Page 21 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 3. Chemical potentials of water in ice Ih (black), CS-I clathrate hydrate (red and purple), and CS-II clathrate hydrate (light blue and blue) against temperature at pressure

p = 0.1 MPa. See text for the distinction between “empty” and “effective reference” clathrate hydrates.

So far, we assume that a guest molecule moves in the force field due to all the surrounding water molecules without affecting the host lattice sites. However, this assumption is not necessarily correct because the vibrational frequencies of the host lattice are influenced by the existence of a large guest molecule. This contribution is calculated by adding the second derivatives of the host-guest interaction potential ψ

∂ 2ψ appearing in eq. 38, ( χ i and ζ i are either Cartesian coordinate or Euler ∂χ i ∂ ζ i angle of a water molecule labelled i) to the Hessian matrix of the size of 6 N w × 6 N w 17,21

(for a water model with the rigid framework such as the TIP4P family34 ). Hence,

the reference state, which is originally the empty clathrate hydrate, should be replaced



by the clathrate hydrate whose lattice vibrations are altered by the guests. Incorporation of the effect by guest molecules usually exhibits blue shift of the vibrational frequencies. This new reference state is referred to as the “effective reference” clathrate hydrate and it does not contain the free energy of cage occupation calculated from either eq 38 or eq 39. This shift is more significant for a large guest and has to be taken into consideration in evaluating the stability of clathrate hydrates.17-19 The Gibbs free energies of the effective reference CS-I and CS-II clathrate hydrates at 0.1 MPa are also shown in Figure 3. In effect, the existence of the guest molecules increases the Gibbs free energy for the CS-I clathrate hydrate by about 0.1 kJ mol-1 and the increase is somewhat enhanced for the CS-II clathrate hydrate. This is because the relative number of the small cages in the CS-II clathrate hydrate is more than that in the CS-I clathrate hydrate and the frequency shift is obviously more pronounced by the small cages. In the following calculation, we replace the chemical potential of the empty clathrate hydrate by that of the effective reference clathrate hydrate in eq 21.

B. Free Energies of Cage Occupation The quantum free energy of cage occupation is obtained by a set of the energy levels of a guest with the radial potential V(r), which is an approximate form of the interaction with the surrounding water molecules fixed to the lattice sites as described previously.42 The radial potential energy curves and the energy levels for CS-I clathrate hydrate are shown in Figure 4. Alternately, the classical free energy of cage occupation is calculated simply by integration of the Boltzmann factor of the interaction of a guest with all the host water molecules. The classical free energy, eq 38, is compared with the quantum one according to eq 39 in Figure 5. Unlike the free energy of ice, the classical free


Page 22 of 42

Page 23 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


energy value is almost the same as the quantum one and the slope for the classical free energy is negative except for the very low temperatures. Although the minimum value of V(r) in a large cage (14-hedral) and a small cage (12-hedral) is similar to each other as shown in Figure 4, the free energy of cage occupation for the large cage is lower than that for the small cage. The large cage is more likely to be occupied than the small cage.

Figure 4. Potentials (black curves) against radial distance r and energy levels (horizontal lines) of a methane molecule in (a) a large cage and (b) a small cage of CS-I clathrate hydrate with unit cell dimension of 1.203 nm.



Figure 5. Free energies of cage occupation by a methane molecule in the large cage (black) and the small cage (red) of CS-I clathrate hydrate. The free energies are calculated according to classical (solid line) and quantum (dotted line) mechanics.

C. Dissociation Pressure of Methane Clathrate Hydrate A stringent test of the methodology is to examine reproducibility of the experimental dissociation pressure as a function of temperature. The chemical potential values of the TIP4P/ice model in liquid water and ice Ih are required to explore the phase behavior of clathrate hydrate in a wide range of temperature. Those quantities associated with the phase transition of water are available only on the basis of the classical mechanical free energies of liquid and solid water.34 Therefore, we rely mostly on the calculations based on the classical free energy and comparison is made with that on the quantum one. The phase boundary between liquid water and ice Ih is estimated approximately by the melting temperature at 1 bar and the slope of the boundary near the melting point at 1 bar. The chemical potential of water is calculated using the chemical potential of ice at


Page 24 of 42

Page 25 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the melting point as described in Sec. III. In Figure 6, the dissociation pressure is compared with the experimental ones for the two sorts of methane models listed in Table I.47-51 Either LJ parameter set recovers the experimental temperature dependence fairly well considering simplicity of the functional form of the potential for guest methane. The experimental measurements are reproduced more satisfactorily with the larger LJ size parameter of 0.3822 nm.35 In general, a larger ε enhances the stability of clathrate hydrate and a larger σ diminishes it. For a further agreement with the experiment, we may make a slight adjustment of the LJ parameters in future. Another serious test is whether methane forms not CS-II but CS-I. This requires more favorable accommodation of methane in the large cage (14-hedron in CS-I, 16-hedron in CS-II) than in the small (12-hedral) cage for both CS-I and CS-II and/or similar stability of the empty CS-I to the empty CS-II. Although the empty CS-II is known to be slightly more stable than the empty CS-I,2 the chemical potential difference between the empty CS-I and CS-II must not be large for CS-I encaging methane to be preferentially formed. In fact, the difference is small as shown in Figure 3 and it is even smaller for the “effective reference” clathrate hydrate in which the influence by the guest molecules on the host vibrational motions is taken into account. Therefore, preferential formation of CS-I clathrate hydrate encaging guest methane is realized as seen in Figure 6: The dissociation pressure of CS-I clathrate hydrate is always lower than that of CS-II, though the difference is small. We may calculate both the quantum free energy of the host lattice and the cage occupation, although the phase behavior of pure water for the TIP4P/ice model has been estimated only by classical MD simulations. The quantum free energies lead to the



similar dissociation pressure of methane clathrate hydrate but it is a little more stable as shown in Figure 6.

Figure 6. Dissociation pressures of methane clathrate hydrate against temperature for two kinds of LJ methane potentials together with experimental measurements. Solid and dotted lines denote CS-I and CS-II, respectively. C and Q indicate the free energies are calculated according to the classical and quantum mechanics.

D. Phase Diagram of Methane Clathrate Hydrate in Temperature-Composition Space The region where CS-I methane clathrate hydrate is the only stable phase is surrounded by the hydrate/water and hydrate/guest boundaries. These boundaries are estimated according to either eq 21 or eq 23. The intersection of them at a given pressure corresponds to the dissociation temperature, which is plotted in Figure 6. Figure 7 illustrates the stable regions of methane clathrate hydrate at pressures, 1, 2.5,


Page 26 of 42

Page 27 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


10, 100, and 200 MPa in the temperature-composition plane. A more detailed illustration is given in Figure S1. The solid lines indicate the phase boundaries between water (ice) and methane clathrate hydrate at various pressures. The dotted lines are those between methane clathrate hydrate and methane fluid. Methane clathrate hydrate is the only stable phase in the region surrounded by the left solid and the right dotted lines at each pressure. The left boundary at a low pressure such as 1 MPa corresponds to equilibrium between clathrate hydrate and ice Ih, where the dissociation temperature is not higher than the melting point of ice. The shape of the enclosed region is similar to what was obtained previously.28 The left boundary at a certainly high pressure such as 10 MPa has a sharp kink caused by the phase transition of water from ice to liquid. Depicted in Figure 8 are the chemical potentials of either liquid water or crystalline ice Ih in equilibrium with the methane clathrate hydrate and the “effective reference” clathrate hydrate, the latter value of which must be decreased by filling methane molecules to meet the former one at the equilibrium as conducted from eq 21.

The

chemical potential of liquid water decreases more harshly against temperature than ice and the difference from the effective reference clathrate hydrate becomes larger as shown in Figure 8. For equilibrium between water and methane clathrate hydrate, the chemical potential of water in the methane clathrate hydrate must equal to that of the pure water phase. It follows that the occupancy increases with increasing temperature above the melting point of ice. This results in the positive slope of the boundary in Figure 7. If the phase transition between ice and liquid water is ignored and the chemical potential of ice Ih is extended to a higher temperature above the melting point, the stable region of the methane clathrate hydrate would be unphysically large.



Page 28 of 42

Figure 7. Phase boundaries of clathrate hydrate with water (solid) and with methane fluid (dotted) at pressure p = 1, 2.5, 10, 100, and 200 MPa in temperature-composition diagram. Methane clathrate hydrate is the only stable phase in the region surrounded by the solid and dotted lines at each pressure.

The free energies are calculated at a fixed

pressure with variable molar volume (heavy lines) and at a fixed molar volume (thin lines).

On the other hand, the hydrate/guest boundary has a negative slope in the temperature-composition diagram for all the pressures.

The slope is calculated from

eqs 4, 14, and 15 as

 ∂y   p=  ∂T 

∑α j

j

 ∂{( µ g − f j ) / T }  x j (1 − x j )   ∂T     k B  ∑ α j x j + 1  j 

2

(40)

The negative slope is caused by the dominant interaction energy of methane in the clathrate hydrate over that in the fluid phase. That is, the sign is determined by the


Page 29 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


enthalpy difference between its fluid, hg, and the encaged state, hj, as

∂{( µ g − f j ) / T } ∂T

=−

hg − h j T2

.

(41)

We expect that hg > h j due to the low density of the fluid methane under moderate pressures. The hydrate/water and hydrate/guest boundaries have always different slopes at the dissociation temperature exhibiting a cusp-like feature and cannot constitute a region enclosed by a smooth convex curve in the composition-temperature diagram.26

Figure 8. Chemical potentials of water against temperature in methane clathrate hydrate (solid line), which are equivalent to the chemical potential of pure water (ice) at equilibria, and those of the “effective reference” clathrate hydrate (dotted line) at 10 (purple), 100 (blue), and 200 (red) MPa.

The hydrate/guest boundary has the upper bound in the methane mole fraction, 8/(46+8) and the boundary is close to the limiting value even at a moderate pressure. On the other hand, the hydrate/water boundary can afford to move toward higher mole



fraction of guest methane under high compression because the difference in the chemical potential of water between liquid water (or ice) and clathrate hydrate increases with increasing pressure, which is caused in part by the larger contribution from the pV term. As a result, the area enclosed by the two boundaries shrinks with increasing pressure above 10 MPa as shown in Figure 7. It is of interest to compare the area of the single stable phase of clathrate hydrate at a constant pressure with that at a constant volume. For the latter calculation, the molar volume is fixed to a standard value, Vm = 22.79 cm3, which corresponds to the unit cell dimension of 1.203 nm. The area with the fixed volume apparently becomes smaller with increasing pressure as shown in Figure 7, suggesting that calculation under constant pressure is essential in drawing a composition-temperature diagram, especially at low temperatures. The cage occupancies, x, at pressures 2.5 MPa and 100 MPa are plotted in Figure 9a. Methane clathrate hydrate is in equilibrium with either water or methane fluid. This is distinguished by the color of curves in Figure 9a. The occupancies for the large and small cages are also shown in Figures 9b and 9c at pressures ranging from 1.6 MPa to 204.8 MPa, although the coexistence phase is not indicated in these panels. The occupancy for the large cage is almost unity at a high pressure but it is a little lower in the small cage in equilibrium with water. A more detailed temperature-pressure dependence of the occupancies in equilibrium with water is given in Figures S2 and S3. The occupancy value increases with increasing the pressure both for the small and large cages as is evident from the pressure dependence of the chemical potential of guest methane in eq 6. The occupancies decrease with increasing temperature when the hydrate is in equilibrium with ice but increase when it is in equilibrium with liquid


Page 30 of 42

Page 31 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


water. This is elucidated by the temperature dependence of the chemical potential difference of water as discussed above: The difference between liquid water and empty clathrate hydrate increases with increasing temperature while the difference between ice and empty clathrate hydrate decreases slightly as shown in Figure 8.



Page 32 of 42

Figure 9. Occupancies of cages against temperature, which are drawn only in the temperature range where methane clathrate hydrate is formed.

(a) large cage (solid

line) and small cage (dotted line) against temperature in equilibrium with water at pressure p = 2.5 (blue) and 100 MPa (red) and those in equilibrium with fluid methane at p = 2.5 (light blue) and 100 MPa (purple). Green lines indicate those occupancies at p = 100 MPa with fixed molar volume. (b) Occupancies of the large cage at p = 1.6, 3.2, 6.4, 12.8, 25.8, 51.2, 102.4, 204.8 MPa. (c) Occupancies of the small cage at the same pressures as in (b). Note that the scales of the ordinate axes in panel b and c are different.

Another new feature obtained in the present calculation is the chemical potential of methane in methane clathrate hydrate. Those in the stable clathrate hydrate at several temperatures are depicted in Figure 10 as a function of composition. The chemical potential increases with increasing the methane ratio and each right endpoint is


Page 33 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


equivalent to the chemical potential of methane fluid in equilibrium with the clathrate hydrate as given in eq 23. The highest occupancy (each right endpoint in Figure 10) increases with cooling and approaches the limiting value of 8/(8+46), which is realized by a certain high chemical potential of methane. However, rapid decrease from the right endpoint is observed and there exists a fairly large margin to the left endpoint where the clathrate hydrate is in equilibrium with water.

Figure 10. Chemical potentials of methane inside methane clathrate hydrate at pressure p = 10 MPa and temperature T = 220 (purple), 240 (blue), 260 (green), 280 K (red). Each right endpoint corresponds to the chemical potential of methane fluid.

An equilibrium between clathrate hydrate and liquid water is of particular interest since clathrate hydrates of natural gases usually exist in (sea) water. This boundary is obtained by our calculation as depicted in Figure 7. Some studies on the equilibrium have been carried out using empirical methods such as the equation of state and the so called Langmuir constants which are associated with the free energies of cage



occupation.52-53 They anticipated that the methane concentration of the aqueous phase in the presence of the methane clathrate hydrate is lower than that in the artificial equilibrium with methane fluid without the intervening clathrate hydrate and the solubility of methane decreases with decreasing temperature in the presence of the clathrate hydrate. This temperature dependence of the solubility was confirmed experimentally54, which is opposite to the conventional trend in hydrophobic hydration of gaseous solute. In our calculation at 280 K, the chemical potential of methane in clathrate hydrate at the left endpoint (in equilibrium with liquid water) is fairly lower than that at the right endpoint (in equilibrium with methane fluid) as shown in Figure 10. This seems to provide a firm ground for a future research on the observed temperature dependence of the solubility.54

V. CONCLUDING REMARKS A theoretical method is developed based on the vdWP theory to calculate the phase diagram of clathrate hydrate in temperature-pressure-composition space. Emphasis is placed on estimating a region where a clathrate hydrate is the only stable phase. Instead of the original vdWP theory with a semi-open system, the thermodynamic stability of the closed system containing fixed numbers of host water and guest species is examined and the conditions for phase equilibrium is estimated for methane clathrate hydrate using currently available force fields. The chemical potentials of crystalline solids (ice and empty clathrate hydrates) are calculated by the lattice dynamics and the subsequent quasi-harmonic approximation. A slight manipulation makes it possible to estimate the chemical potential of liquid water from only the thermodynamic properties on the


Page 34 of 42

Page 35 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


transition of ice to liquid water without calculating the anharmonic vibrational free energy. This simple but efficient way allows us to extend the range of temperature in the phase diagram. The present method is applied to methane clathrate hydrate, which results in an excellent agreement in the dissociation pressure with the experimental observations. We find that the phase boundary between water and clathrate hydrate is significantly interfered by the phase transition of the coexisting host water from ice to liquid, which puts a limit on the area of the stable clathrate hydrate in temperature-composition space at pressures higher than 2.5 MPa. We also reveal that it is essential to choose the pressure as an independent variable so as to recover the accurate phase equilibria in the composition space. The present method establishes a temperature-composition relation when a clathrate hydrate containing any guest species is the only stable phase at a given pressure. Our calculation may be improved by considering the anharmonic contribution to the free energy of the intermolecular vibrations. At this moment, we expect optimistically cancelation of this term between ice and clathrate hydrate.

In order to check the

validity of the approximation roughly, we estimate the individual anharmonic vibrational free energy by the thermodynamic integration method at a fixed volume with the reference of the harmonic system. Those for ice Ih and empty CS-I are calculated to be -1.28 and -1.28 kJ mol-1 at 250 K and to be -1.73 and -1.66 kJ mol-1 at 300 K, respectively.

The difference is smaller than 0.1 kJ mol-1, which can alter the

dissociation pressure with minor degree even at the higher temperature of 300 K. Our series of theories17-21 have eliminated various restrictions in the original vdWP theory and extended the applicability to broader targets with higher precisions. The



present research incorporates further extensions and brings the theory to completion. Our method presents quantitative features of stability of clathrate hydrates with high reproducibility of various thermodynamic and structural properties in their pure states. Since the pressure is one of the independent properties and the volume (cell dimension) can vary in the present method, it is applicable to a clathrate hydrate (or filled ice) under very high pressure such as 1 GPa with the knowledge only on the stable states of aqueous (ice) and guest-rich phases.14, 55-57

Another merit is that our method does not

rely on MD or MC simulations which require heavy computational cost and suffer from possible statistical errors; Our method can handle many hydrogen-disordered structures of crystalline states.

Furthermore, the evaluation of stability via calculation of the free

energy of cage occupation makes it easier to understand the origin of the stability of clathrate hydrates in a wide range of thermodynamic condition by fine control of temperature and pressure.

Supporting Information More detailed phase behaviors and occupancies are shown in the Supporting Information. The phase boundaries of methane clathrate hydrate with water and fluid methane are drawn in temperature-pressure-composition space in Figure S1. The occupancies of guest methane in large and small cages of methane clathrate hydrate in equilibrium with water are depicted in Figures S2 and S3.


Page 36 of 42

Page 37 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Author Information Corresponding Author Phone: +81-(0)86-251-7769 E-mail: [email protected]

Acknowledgments The present work was supported by JSPS KAKENHI Grant Number 17K19106 and MEXT as “Priority Issue on Post-Kcomputer” (Development of new fundamental technologies for high-efficiency energy creation, conversion/storage and use). The authors are grateful to Prof. M. Taniguchi, Prof. K. Koga, and T. Nakayama for valuable discussions.

References 1.

Davidson, D. W., Water - a Comprehensive Treatise; Plenum: New York, 1973; Vol. 2.

2.

Sloan, E. D.; Koh, C. A., Clathrate Hydrates of Natural Gases; CPC Press: Boca Raton, 2008.

3.

Franks, F., Water - a Comprehensive Treatise; Plenum: New York, 1972; Vol. 2.

4.

Guillot, B.; Guissani, Y.; Bratos, S., A Computer ‐ Simulation Study of Hydrophobic Hydration of Rare Gases and of Methane. I. Thermodynamic and

5.

Structural Properties. J. Chem. Phys. 1991, 95, 3643-3648. Tanaka, H.; Nakanishi, K., Hydrophobic Hydration of Inert Gases:

6.

Thermodynamic Properties, Inherent Structures, and Normal‐Mode Analysis. J. Chem. Phys. 1991, 95, 3719-3727. Tse, J. S.; Klein, M. L., Dynamic Properties of the Structure Ii Clathrate

7.

Hydrate of Krypton. J. Phys. Chem. 1987, 91, 5789-5791. Walsh, M. R.; Koh, C. A.; Sloan, E. D.; Sum, A. K.; Wu, D. T., Microsecond Simulations of Spontaneous Methane Hydrate Nucleation and Growth. Science

2009, 326, 1095-1098. 37 ACS Paragon Plus Environment


8.

Walsh, M. R.; Beckham, G. T.; Koh, C. A.; Sloan, E. D.; Wu, D. T.; Sum, A. K., Methane Hydrate Nucleation Rates from Molecular Dynamics Simulations: Effects of Aqueous Methane Concentration, Interfacial Curvature, and System

9.

Size. J. Phys. Chem. C 2011, 115, 21241-21248. Yagasaki, T.; Matsumoto, M.; Andoh, Y.; Okazaki, S.; Tanaka, H., Effect of Bubble Formation on the Dissociation of Methane Hydrate in Water: A

10.

Molecular Dynamics Study. J. Phys. Chem. B 2014, 118, 1900-1906. Conde, M. M.; Vega, C., Determining the Three-Phase Coexistence Line in Methane Hydrates Using Computer Simulations. J. Chem. Phys. 2010, 133, 064507.

11.

Katsumasa, K.; Koga, K.; Tanaka, H., On the Thermodynamic Stability of

12.

Hydrogen Clathrate Hydrates. J. Chem. Phys. 2007, 127, 044509. Papadimitriou, N. I.; Tsimpanogiannis, I. N.; Papaioannou, A. T.; Stubos, A. K., Evaluation of the Hydrogen-Storage Capacity of Pure H2 and Binary H2-Thf Hydrates with Monte Carlo Simulations. J. Phys. Chem. C 2008, 112, 10294-10302.

13.

Papadimitriou, N. I.; Tsimpanogiannis, I. N.; Peters, C. J.; Papaioannou, A. T.; Stubos, A. K., Hydrogen Storage in Sh Hydrates: A Monte Carlo Study. J. Phys.

14.

Chem. B 2008, 112, 14206-14211. Hakim, L.; Koga, K.; Tanaka, H., Phase Behavior of Different Forms of Ice

15.

Filled with Hydrogen Molecules. Phys. Rev. Lett. 2010, 104, 115701. Nakayama, T.; Koga, K.; Tanaka, H., Augmented Stability of Hydrogen Clathrate Hydrates by Weakly Polar Molecules. J. Chem. Phys. 2009, 131, 214506.

16.

van der Waals, J. H.; Platteeuw, J. C., Clathrate Solutions. Adv. Chem. Phys.

17.

1959, 2, 1. Tanaka, H.; Kiyohara, K., On the Thermodynamic Stability of Clathrate

18.

Hydrate .I. J. Chem. Phys. 1993, 98, 4098-4109. Tanaka, H.; Kiyohara, K., The Thermodynamic Stability of Clathrate Hydrate. II. Simultaneous Occupation of Larger and Smaller Cages. J. Chem. Phys.

19.

1993, 98, 8110-8118. Tanaka, H., The Thermodynamic Stability of Clathrate Hydrate. III. Accommodation of Nonspherical Propane and Ethane Molecules. J. Chem.

20.

Phys. 1994, 101, 10833-10842. Tanaka, H.; Nakatsuka, T.; Koga, K., On the Thermodynamic Stability of Clathrate Hydrates IV: Double Occupancy of Cages. J. Chem. Phys. 2004, 121,


Page 38 of 42

Page 39 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


5488-5493. 21.

Tanaka, H.; Matsumoto, M., On the Thermodynamic Stability of Clathrate Hydrates V: Phase Behaviors Accommodating Large Guest Molecules with

22.

New Reference States. J. Phys. Chem. B 2011, 115, 14256-14262. Belosludov, V. R.; Lavrentiev, M. Y.; Dyadin, Y. A., Theory of Clathrates. J.

23.

Inclusion Phenom. Mol. Recognit. Chem. 1991, 10, 399-422. Zubkus, V. E.; Tornau, E. E.; Belosludov, V. R., Theoretic Physicochemical Problems of Clathrate Compounds. In Adv. Chem. Phys., John Wiley & Sons, Inc.: 1991; Vol. 81, pp 269-359.

24.

Kobayashi,

R.;

Katz,

D.,

Vapor-Liquid

Equilibria

for

Binary

25.

Hydrocarbon-Water Systems. Ind. Eng. Chem. 1953, 45, 440-446. Dendy Sloan, E., Clathrate Hydrate Measurements: Microscopic, Mesoscopic,

26.

and Macroscopic. J. Chem. Thermodyn. 2003, 35, 41-53. Huo, Z.; Hester, K.; Sloan, E. D.; Miller, K. T., Methane Hydrate

27.

Nonstoichiometry and Phase Diagram. AlChE J. 2003, 49, 1300-1306. Nezbeda, I.; Smith, W. R.; Kolafa, J., Molecular Theory of Phase Equilibria in Model Associated Mixtures. I. Binary Mixtures of Water and a Simple Fluid. J.

28.

Chem. Phys. 1994, 100, 2191-2201. Wierzchowski, S. J.; Monson, P. A., Calculating the Phase Behavior of Gas-Hydrate-Forming Systems from Molecular Models. Ind. Eng. Chem. Res.

29.

2006, 45, 424-431. Wierzchowski, S. J.; Monson, P. A., Calculation of Free Energies and Chemical Potentials for Gas Hydrates Using Monte Carlo Simulations. J. Phys. Chem. B

30.

2007, 111, 7274-7282. Petrenko, V. F.; Whitworth, R. W., Physics of Ice; Oxford University Press: Oxford, 1999.

31.

Falenty, A.; Hansen, T. C.; Kuhs, W. F., Formation and Properties of Ice Xvi Obtained by Emptying a Type Sii Clathrate Hydrate. Nature 2014, 516, 231-233.

32.

Waage, M. H.; Vlugt, T. J. H.; Kjelstrup, S., Phase Diagram of Methane and Carbon Dioxide Hydrates Computed by Monte Carlo Simulations. J. Phys.

33.

Chem. B 2017, 121, 7336-7350. Rowlinson, J. S.; Swinton, F. L., Liquids and Liquid Mixtures, 3 ed.; Butterworth: London, 1982.

34.

Abascal, J. L. F.; Sanz, E.; García Fernández, R.; Vega, C., A Potential Model for the Study of Ices and Amorphous Water: Tip4p/Ice. J. Chem. Phys. 2005,



122, 234511. 35.

Hirschfelder, J. O.; Curtiss, C. F.; Bird, R. B., Molecular Theory of Gases and Liquids; Wiley: New York, 1954.

36.

Jorgensen, W. L.; Madura, J. D.; Swenson, C. J., Optimized Intermolecular Potential Functions for Liquid Hydrocarbons. J. Am. Chem. Soc. 1984, 106, 6638-6646.

37.

Matsumoto, M.; Yagasaki, T.; Tanaka, H., Genice: Hydrogen-Disordered Ice Generator. J. Comput. Chem. 2018, 39, 61-64.

38.

Stillinger, F. H.; Weber, T. A., Hidden Structure in Liquids. Phys. Rev. A 1982, 25, 978-989.

39.

Pohorille, A.; Pratt, L. R.; LaViolette, R. A.; Wilson, M. A.; MacElroy, R. D., Comparison of the Structure of Harmonic Aqueous Glasses and Liquid Water. J.

40.

Chem. Phys. 1987, 87, 6070-6077. Tanaka, H., Thermodynamic Stability and Negative Thermal Expansion of

41.

Hexagonal and Cubic Ices. J. Chem. Phys. 1998, 108, 4887-4893. Pauling, L., The Structure and Entropy of Ice and of Other Crystals with Some

42.

Randomness of Atomic Arrangement. J. Am. Chem. Soc. 1935, 57, 2680-2684. Yagasaki, T.; Matsumoto, M.; Tanaka, H., Anomalous Thermodynamic Properties of Ice XVI and Metastable Hydrates. Physical Review B 2016, 93, 054118.

43.

Crc Handbook of Chemistry and Physics, 89 ed., 2008.

44.

Redlich, O.; Kwong, J. N. S., On the Thermodynamics of Solutions. V. An Equation of State. Fugacities of Gaseous Solutions. Chem. Rev. 1949, 44, 233-244.

45.

Dantl, G., Wärmeausdehnung Von H2O- Und D2O-Einkristallen. Z. Phys. 1962, 166, 115-118.

46.

Eisenberg, D.; Kauzmann, W., The Structures and Properties of Water Oxford: London, 1969.

47.

Deaton, M. W.; Frost, E. M., Gas Hydrates and Their Relation to the Operation

48.

of Natural-Gas Pipe Lines. U.S. Bureau of Mines Monograph 1946, 8, 1-101. De Roo, J. L.; Peters, C. J.; Lichtenthaler, R. N.; Diepen, G. A. M., Occurrence of Methane Hydrate in Saturated and Unsaturated Solutions of Sodium Chloride and Water in Dependence of Temperature and Pressure. AlChE J.

49.

1983, 29, 651-657. Nakano, S.; Moritoki, M.; Ohgaki, K., High-Pressure Phase Equilibrium and Raman Microprobe Spectroscopic Studies on the Methane Hydrate System. J.


Page 40 of 42

Page 41 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


50.

Chem. Eng. Data 1999, 44, 254-257. Nakamura, T.; Makino, T.; Sugahara, T.; Ohgaki, K., Stability Boundaries of Gas

Hydrates

Helped

by

Methane—Structure-H

Hydrates

of

Methylcyclohexane and Cis-1,2-Dimethylcyclohexane. Chem. Eng. Sci. 2003, 58, 269-273. 51.

Mohammadi, A. H.; Anderson, R.; Tohidi, B., Carbon Monoxide Clathrate Hydrates: Equilibrium Data and Thermodynamic Modeling. AlChE J. 2005, 51, 2825-2833.

52.

Handa, Y. P., Effect of Hydrostatic Pressure and Salinity on the Stability of Gas

53.

Hydrates. J. Phys. Chem. 1990, 94, 2652-2657. Zatsepina, O. Y.; Buffett, B. A., Phase Equilibrium of Gas Hydrate: Implications for the Formation of Hydrate in the Deep Sea Floor. Geophys. Res.

54.

Lett. 1997, 24, 1567-1570. Servio, P.; Englezos, P., Measurement of Dissolved Methane in Water in

55.

Equilibrium with Its Hydrate. J. Chem. Eng. Data 2002, 47, 87-90. Vos, W.; Finger, L.; Hemley, R.; Mao, H.-k., Novel H2-H2o Clathrates at High

56.

Pressures. Phys. Rev. Lett. 1993, 71, 3150-3153. Chou, I.-M.; Sharma, A.; Burruss, R. C.; Shu, J.; Mao, H.-k.; Hemley, R. J.; Goncharov, A. F.; Stern, L. A.; Kirby, S. H., Transformations in Methane

57.

Hydrates. Proc. Natl. Acad. Sci. U.S.A. 2000, 97, 13484-13487. Hirai, H.; Uchihara, Y.; Fujihisa, H.; Sakashita, M.; Katoh, E.; Aoki, K.; Nagashima, K.; Yamamoto, Y.; Yagi, T., High-Pressure Structures of Methane Hydrate Observed up to 8 Gpa at Room Temperature. J. Chem. Phys. 2001, 115, 7066-7070.



Table of contents image


Page 42 of 42

On the Thermodynamic Stability of Clathrate Hydrates VI: Complete

Recommend Documents