Theory and Kinetic Monte Carlo Simulation of Guest Molecule

Subscriber access provided by UNIV AUTONOMA DE COAHUILA UADEC

C: Physical Processes in Nanomaterials and Nanostructures

Theory and Kinetic Monte Carlo Simulation of Guest Molecule Transport in sI Clathrate Hydrates Based on Cage Hopping Lee-Shin Chu, David T. Wu, and Shiang-Tai Lin J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.9b01109 • Publication Date (Web): 08 Apr 2019 Downloaded from http://pubs.acs.org on April 8, 2019

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 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 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.

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 28 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

Theory and Kinetic Monte Carlo Simulation of Guest Molecule Transport in sI Clathrate Hydrates Based on Cage Hopping Lee-Shin Chua, David T. Wub* and Shiang-Tai Lina* aDepartment

of Chemical Engineering, National Taiwan University, Taipei 10617, Taiwan bDepartment of Chemical and Biological Engineering and Department of Chemistry Colorado School of Mines, Golden, CO 80401, USA *Corresponding authors: [email protected], [email protected]

Abstract Guest migration in clathrate hydrates is a slow but important process for reaching thermodynamic equilibrium. The transport of guest molecules in a hydrate lattice is considered as a series of hopping events from an occupied cage to an empty neighboring cage facilitated by water vacancies and without significant lattice restructuring in the bulk. In this work, we developed an analytical model for determining the equilibrium distribution and the diffusivity of gas molecules in the cages of sI clathrate hydrate based on their hopping rate. Furthermore, kinetic Monte Carlo (KMC) simulations were performed to verify the analytical model. The equilibrium occupancies, transport (Fickian) and jump (Maxwell-Stefan) diffusion coefficients and the thermodynamic correction factors determined from the analytical model are in excellent agreement with the simulation results. Using the hopping rate constants determined based on transition path sampling calculations, we obtain methane transport diffusion coefficient at 275 K to be 5.06×10-14 m2/s, which is in good agreement with the recent experimental measurements (4.00×10-14 m2/s at 275 K). In addition, the analytical model is useful for identifying the predominant cage hopping pathways for equilibrium and transport properties.

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 28

1. Introduction Solid state diffusion of methane in clathrate hydrates is the rate limiting step of several important processes, including the formation of methane hydrate from ice powders and the gas replacement of methane hydrates with CO2. Kuhs et al. conducted neutron diffraction and gas-consumption experiments to study the kinetics of methane hydrate formation from hydrogenated and deuterated ice powder samples in the temperature range of 245 ~ 270 K.1 According to the proposed model, the hydrate shell growth consists of two steps. The first step is an interfacial clathration reaction and the second step is the gas and water diffusion through the hydrate layer surrounding the shrinking ice cores. The diffusion process in the second step is the rate limiting step at temperatures above 263 K. The permeation coefficients (or jump diffusivity) of methane were reported to be in the range of 4.17 × 10−17 ~ 3.61 × 10−16 m2/s at the temperature range from 245 to 270 K. The replacement of methane in natural gas hydrates with CO2 is an intriguing approach that can simultaneously extract energy resources and sequester greenhouse gases.2-3 Several experimental works have tested the replacement process in the laboratory scale.4-14 In the Iġnik Sikumi gas hydrate field experiment, mixed CO2 and N2 gas (77.5% N2/22.5% CO2) was injected into the methane reservoirs.15-16 It was reported that along with the methane production, about 70% of the injected N2, but only 40% of the injected CO2, was recovered during the flowback phase. This observation of favored retention of CO2 implies that the bulk exchange (net change in reservoir pore fill regardless of the process by which it occurred) of CO2 for CH4 did occur. Salamatin et al. simulated the gas exchange process by developing a mathematical model for the nonequilibrium binary diffusion of guest molecules.17-19 The whole process of gas exchange is conceived as a two-stage switching, starting with a rapid formation of a mixed hydrate layer on the hydrate surface followed by a much slower diffusion-controlled process. They performed gas replacement experiments with CO2 in CH4-hydrate at a higher temperature range of 274 ~ 280 K.19 Based on their model, the transport diffusivities of CH4 and CO2 were calculated to be 4.00 × 10−14 and 8.58 × 10−14 m2/s, respectively, at 275 K, where the diffusion coefficient of CO2 is approximately two times higher than the diffusion coefficient of CH4. The transport of guest molecules is considered as a series of hopping events between occupied and empty neighboring cages without significant lattice restructuring. Several studies have shown that the energy barrier for guest molecules hopping through cages is significantly reduced if a water vacancy exists in the water ring at the interface shared by the donor and acceptor cages.20-22 With evidence from several simulation studies ACS Paragon Plus Environment

Page 3 of 28 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


that the mobility of water vacancies in the hydrate cages is much larger than that of the guest molecules20, 23, it is reasonable to assume that each hopping event of a guest molecule is facilitated by the presence of water vacancies in the hydrate lattice. Several simulations studies were devoted to determine the transport properties of gas in sI hydrates. Demurov et al. used umbrella sampling and Monte Carlo (MC) simulations to calculate free energy barriers of carbon dioxide hopping through different pathways, including L5L, L6L, L5S, and S5L, with and without water vacancies present in the system. The pathway nomenclature L5L indicates a hopping from a large (L) cage through a pentagonal face to another large cage, and similarly for the other pathways. Transition state theory was used to calculate the hopping rates from the free energy barriers.20 Peters et al. provided a detailed analysis of the water vacancy assisted diffusion in sI methane hydrates.21 They used path sampling methods to calculate the free energy barrier and further reactive flux correlation function corrections to determine hopping rate constants. Kinetic Monte Carlo simulation was used to calculate the diffusivity of methane in 8 × 8 × 8 unit hydrates with vacancies. The self-diffusivity of methane was estimated to be 7 × 10-15 X m2/s at 250 K, where X is the fraction of unoccupied cages. The more recent simulation work by Liang et al. used molecular dynamics (MD) simulation to study the hopping mechanisms of CO2 hydrates at relatively high temperature, close to the dissociation temperature (310 − 320 K at 100 MPa).24 They observed a water vacancy-interstitial defect pair temporarily appearing at the interface of the two cages while a hopping event occurs. The short-lived defect pair is recovered after the hopping event completes. Lo et al. also explored the cage hopping of guest molecules by molecular dynamics (MD) simulation on single phase sI hydrates with water vacancies introduced in the system.22 Without water vacancies present, no guest hopping events were observed. With a water vacancy, they observed the microscopic mechanisms of CO2 cage hopping in stable hydrate within 100 ns simulations. However, there was no methane guest hopping within 100 ns simulations even with water vacancies present. Waage et al. studied the diffusion of gas mixtures in sI hydrate by umbrella sampling and MC simulations.25 The free energy barriers of CH4, CO2, N2 and H2 passing through different pathways, with and without water vacancies present were calculated. The computed free energy barriers for CH4 agree with those from Peters et al.,21 while the results for CO2 are significantly higher compared to the results from Demurov et al.20 The disparity might arise from the different molecular models of water used in their simulations (Waage et al. used TIP4P/Ice and Demurov et al. used SPC/E). ACS Paragon Plus Environment


Page 4 of 28

Guest hopping between cages in hydrates is a rare event even when water defects exist in the system. In a traditional molecular dynamics (MD) simulation, few hopping events would be observed (for example, approximately 11 hopping events will happen within 1 s and a box size of 1000 nm3, estimated based on the hopping rate constants at 250 K21), which might result in poor sampling for calculating the diffusivity of a guest molecule in the solid state. Kinetic Monte Carlo (KMC) simulations may circumvent such problems as long as the rate constants of the possible events are known. The thermodynamic and transport properties of guest molecules in a hydrate can then be calculated by properly chosen parameter settings of the simulation systems. In order to gain physical insights into the thermodynamic and transport properties of the cage hopping model or guests studied in such KMC simulations, we developed an analytical model according to the Maxwell-Stefan approach to mass transfer,26 which links the transport-diffusivity of the guest molecules to the hopping rate constants and the total occupancy in sI hydrates. Similar approaches have been used to model guest diffusion in porous crystalline materials such as zeolites, carbon nanotubes and metal-organic frameworks.27-28 For example, Coppens et al. derived analytical expressions for occupancy and self-diffusivity in zeolites consisting of weak and strong adsorption sites based on mean field theory.29-30 In this work, proper descriptions of Maxwell-Stefan diffusivity and thermodynamic correction factors for sI hydrate are developed. Similar to the zeolites studied by Coppens et al., sI hydrate contains two types of cages; however, unlike those zeolites, the pathways between the cages are different and our focus here includes transport-diffusivity of gas molecules in addition to self-diffusivity. To the best of our knowledge, this theoretical approach for solid-state diffusion has not previously been applied to hydrate systems. We show that the analytical model allows for the direct connection between the hopping rate constants and the derived thermodynamic and transport properties, providing useful insights to the understanding of macroscopic properties based on molecular hopping events.

2. Theory Molecular model. The unit cell of sI hydrate contains 6 large (L) cages and 2 small (S) cages. Guest molecules residing in the cages can hop to one of its unoccupied neighboring sites. Thus the number of possible hopping events of a guest molecule depends on the type of cage it resides in and the occupancy of its ACS Paragon Plus Environment

Page 5 of 28 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


neighboring sites. Four kinds of pathways exist in the sI hydrate system, according to the type of the donor and acceptor cages (small 512 and large 51262) and the number of water molecules in the face they share (a 5 or 6 membered water ring), and are correspondingly named L6L, L5L, L5S and S5L. A large cage (51262) in the hydrate lattice has 14 neighboring cages: ten large cages and four small cages. Two of the ten neighboring large cages are connected via a six-membered water ring, and the other eight large cages are connected via a five-membered water ring. A small cage (512) in the hydrate lattice only has 12 neighboring large cages: all connected via a five-membered water ring. The hopping rate constants depend on the temperature of the system and the types of pathways a molecule hops through. Table 1 summarizes the guest hopping pathways and the rate constants of sI methane hydrate. Aside from not being able to simultaneously occupy the same cage, the mutual interaction of the guest molecules inside the cages is neglected in this work, hence the local concentration of empty sites will not affect the hopping rate constants of the guest molecules. The hopping events are thus considered to be statistically independent.

Table 1. Summary of the parameters used in describing the sI methane hydrate system.

Cage type

Type of neighboring cage

Number of neighboring cages

Interface water geometry

Type of pathway (i)

Number of pathways (ni)

Hopping distance (di) (nm)

512

S

0

-

-

0

-

225 K

250 K

(S)

L

12

5-ring

S5L

12

0.6652

5.3

35

S

4

5-ring

L5S

4

0.6652

0.14

2.6

5-ring

L5L

8

0.7299

0.07

1.5

L

10 6-ring

L6L

2

0.6003

31.8

101

51262 (L)

Hopping rate constant21 (ki) (1/ms)

Equilibrium guest distribution in sI hydrate lattice. Consider a perfect sI clathrate hydrate consisting of NS small and NL large cages (NS:NL=1:3) that is filled with N methane molecules. Supposing OS and OL methane molecules reside in the small and large cages, respectively, (with one methane per cage, OS+OL=N), the fractional occupancies of large and small cages are OL

OS

XL = NL, XS = NS.

[1] ACS Paragon Plus Environment


Page 6 of 28

and the total occupancy is Xt =

OL + OS NL + NS

NL

NS

3

1

[2]

= NL + NSXL + NL + NSXS = 4XL + 4XS,

The hopping rate constants of the four different pathways are denoted as kL6L, kL5L, kL5S and kS5L. The number of neighboring pathways of a cage (large or small) are described as nL6L, nL5L, nL5S and nS5L. For sI hydrate, nL6L = 2, nL5L = 8, nL5S = 4 and nS5L =12.

We start with the rate of change of the number of guest molecules in the small cages, dOS dt

= kL5S × OL × nL5S × (1 ― XS) ― kS5L × OS × nS5L × (1 ― XL) ,

where we have made a mean-field assumption by using the average occupancies in the rate expression.

[3] As

we are also assuming interactions between guests in neighboring cages are negligible, this should be a reasonable approximation when local equilibrium holds, e.g. when spatial concentration gradients are not too extreme. Substituting eq. 1 into 3 and given the fact that NLnL5S = NSnS5L, we have dOS dt

= [ kL5S × XL × (1 ― XS) ― kS5L × XS × (1 ― XL) ] NLnL5S.

[4]

Since a guest molecule is either inside a large cage or a small cage, the rate of change of the number of guest molecules in the large cages has the relation: dOL dt

=―

dOS dt

[5]

.

At equilibrium,

dOL dt

=

dOS dt

= 0, and eq. 4 becomes

kL5S × XL × (1 ― XS) ― kS5L × XS × (1 ― XL) = 0.

[6]

Eqs. 2 and 6 can be used to determine the fractional occupations in terms of the total occupancy (see Supporting Information for details) XL = XS =

((4Xt ― 1) ― 𝐾(4Xt + 3)) + ((4Xt ― 1) ― 𝐾(4Xt + 3))2 + 48𝐾(1 ― 𝐾)Xt 6(1 ― 𝐾)

((4Xt + 1) ― 𝐾(4Xt ― 3)) ― ((4Xt + 1) ― 𝐾(4Xt ― 3))2 ― 16(1 ― 𝐾)Xt 2(1 ― 𝐾)

;

.

[7] [8]

kS5L

where 𝐾 = kL5S. The hopping rate constants kL6L and kL5L do not appear in the above equations because the hopping along these two pathways do not change the fractional occupancies XL and XS. Therefore, the equilibrium fractional occupancies are independent of these two hopping rate constants. ACS Paragon Plus Environment

Page 7 of 28 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


According to the van der Waals–Platteeuw model,31 the fugacity, f, of encaged gas solutes in thermodynamic equilibrium with the ambient gas atmosphere can be expressed in terms of fractional occupancies as 1

XL

XS

1

[9]

f = CL1 ― XL = CS1 ― XS

where CL and CS are the temperature dependent Langmuir constants of the gas in the large and small cages, respectively. The ratio of the Langmuir constants is related to the hopping rate constant as CL CS

XL(1 ― XS)

kS5L

[10]

= XS(1 ― XL) = kL5S = 𝐾.

where the first equality is derived from eq. 9 and the second equality from eq. 6. This also highlights that the parameter, K, is effectively an equilibrium constant between guests and vacancies in the large and small cages. Note that eq. 9 allows for the calculation of fractional occupancies from the Langmuir constants CLf

C Sf

[11]

XL = 1 + CLf and XS = 1 + CSf.

Combined with eq. 2, the total occupancy can be determined from the gas fugacity in the gas or liquid phase as 3 CLf

1 C Sf

[12]

Xt = 41 + CLf + 41 + CSf.

The transport properties of guest molecules in sI hydrate lattice. For a system not in equilibrium, the molecular flux, expressed in number of particles per square meters per second, is given by Fick’s law32 [13]

q = ―ρDtran∇Xt,

where ρ is the water cage number density of sI hydrate, expressed in number of hydrate cages per m3, Dtran is the transport or Fickian diffusion coefficient. A more general approach using chemical potential as the driving force to describe the diffusion process is given by32 q = ―ρDjump

(

Xt

),

[14]

RT∇T,Pμ

where Djump is referred to as the jump or Maxwell-Stefan diffusion coefficient, R is the gas constant and T is the temperature. The chemical potential μ of the encaged guest molecules is related to its fugacity as μ = RTlnf.

[15]

Substituting eq. 15 into eq. 14 and comparing with eq. 13, we obtain the relation between transport- and jumpdiffusivities as Dtran = ΓDjump,


[16]


Page 8 of 28

where Γ is the thermodynamic correction factor: Γ ≡

∇T,Plnf ∇lnXt

dlnf

[17]

= dlnXt.

The thermodynamic correction factor Γ can be determined as a function of the total occupancy Xt and the ratio 𝐾 of hopping rate constants between large and small cages using eq. 12 (see Supporting Information for details), i.e., 1

Γ = 1 ― Xt +

4Xt(𝐾 ― 1)(1 ― 3𝐾 + 4Xt) + α(1 + 3𝐾 ― α) α(α ― 3𝐾 + 4Xt(1 + 𝐾) ― 1)

[18]

.

where α = 16(𝐾 ― 1)2X2t ― 8(3𝐾 ― 1)(𝐾 ― 1)Xt + (3𝐾 + 1)2.

[19]

1

It is noteworthy that, when 𝐾=1, eq. 18 reduces to Γ = 1 ― Xt, which is the thermodynamic correction factor for a one cage type relation.

The next step is to relate the jump-diffusion coefficient with the total occupancy and the hopping rate constants. The jump-diffusion coefficient Djump is linearly dependent on the total occupancy Xt in a “one cage type” system,32 given as Djump = Djump(0) × (1 ― Xt).

[20]

However, since the sI hydrate system has two types of cages (small and large), eq. 20 is not applicable. Instead, we consider the jump diffusivity as the sum of the contributions to the diffusion from each kind of pathway 1

Djump = 6∑𝑖Pi × kidi2 (with i = L6L, L5L, L5S, S5L).

[21]

where ki and di are the rate constant and hopping distance of a specific pathway (as summarized in Table 1). 1

The term 6kidi2 is analogous to the zero-loading diffusion coefficient Djump(0), and Pi is the probability of a guest molecule occupying a particular cage type hopping to an adjacent empty cage of a particular type through a specific kind of pathway, which depends on the occupancies of neighboring cages and the number of connecting pathways, as shown in Table 2. Note that both eq. 20 and the values of Pi given in Table 2 use mean-field values for the concentration of empty sites. As mentioned earlier, the fractional occupancies of the large and small cages (XL and XS) are functions of the total occupancy Xt and the hopping rate constants. As a consequence, the jump-diffusion coefficient is also a function of Xt and the hopping rate constants. Eq. 21 is a generalization of eq. 20, which takes into account the difference of cage types and the connectivity of the ACS Paragon Plus Environment

Page 9 of 28 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


structure, and will be shown later to agree with the KMC simulation results better than eq. 20 for the jumpdiffusion coefficient in the sI hydrate system.

Finally, the transport-diffusion coefficient as a function of the hopping rate constants and the total occupancy can be obtained from combining eqs. 16 and 21 Dtran =

(

1 1 ― Xt

4Xt(𝐾 ― 1)(1 ― 3𝐾 + 4Xt) + α(1 + 3𝐾 ― α)

+

α(α ― 3𝐾 + 4Xt(1 + 𝐾) ― 1)

)×

1 4 ∑ 6 i = 1Pi

× kidi2.

[22]

Table 2. Probability of a guest molecule occupying a particular cage type hopping to an adjacent empty cage of a particular type, for the four different pathways in sI hydrate. Pathway L6L L5L L5S S5L

𝐏𝐢 NL NL + NS NL NL + NS NL NL + NS NS NL + NS

× nL6L × (1 ― XL) × nL5L × (1 ― XL) × nL5S × (1 ― XS) × nS5L × (1 ― XL)

Thus far, the diffusivities (Djump, Dtran) and the thermodynamic correction factor Γ are determined as a function of the total occupancy Xt and the hopping rate constants (kL6L, kL5L, kL5S and kS5L). The KMC simulation results of the same properties based on the same set of total occupancy and hopping rate constants will be calculated to compare with the results obtained from the analytical model.

To describe the self-diffusivity, the inter-relation between self- and jump-diffusivity for a single component system according to the Maxwell-Stefan formulation is derived by Krishna and Paschek32 as Dself =

1 (D

1

jump

Xt

[23]

,

+D ) i,i

where Djump can be determined according to eq. 21. However, the self-exchange coefficient Di,i, which accounts for the correlation effects due to intermolecular interactions and guest-host topology,33 does not have an exact expression. Krishna and Paschek suggested a practical solution by assuming the self-exchange ACS Paragon Plus Environment


Page 10 of 28

coefficient to be equal to the jump-diffusivity, i.e. Di,i = Djump(Xt), which works well for a few cases.32 Another approach to determine the self-exchange coefficient Di,i is by fitting the KMC simulation results.33 Since the self- and jump- diffusivities can be determined independently by the simulations, the self-exchange coefficient Di,i can be back-calculated. In this work, we applied the approach suggested by Krishna and Paschek32 to estimate the self-diffusivity by assuming that the self-exchange coefficient is equal to the jumpdiffusivity.

Kinetic Monte Carlo method. The kinetic Monte Carlo (KMC) method we adopted in this work is based on the Bortz−Kalos−Lebowitz (BKL) algorithm.34-36 In this scheme, we scan through the simulation lattice cell and create a list including all possible hopping events, i, and their associated rate constants κi at the current system configuration. Random numbers (γ1,γ2) ϵ (0,1] are generated, and the hopping event j is selected satisfying the following condition j

j+1

[24]

∑i = 1κi ≤ γ1κtot ≤ ∑i = 1κi

where the overall rate constant κtot represents the sum of rate constants of all possible hopping events at the current configuration. After the system configuration is adjusted according to the selected event, the time is advanced stochastically to t→t + ∆t where ∆t is given as ∆t = ―

ln (γ2) κtot

[25]

.

From the updated configuration, a new list of possible events is regenerated and the whole procedure starts again. Following this scheme, the simulation system can be propagated step by step until the desired step number is achieved.

Self-diffusivity. Using the random walk theory, the tracer diffusion coefficient, which is a useful quantity to describe mobility of guest molecules, can be extracted from KMC simulations in terms of the mean square displacement,32 1

N

Dself = lim 6Nt∑i = 1〈[ri(t) ― ri(0)]2〉

[26]

t→∞


Page 11 of 28 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


where the vector ri(t) ― ri(0) indicates the change in position of a particle i after passing a time t. N is the number of guest molecules in the system. The angular brackets represent an expectation value which is obtained by averaging over a time interval.

Jump-diffusivity. The jump-diffusivity is calculated by the KMC simulations using the following relation37-38: 1

〈[

N

]〉

Djump = lim 6Nt ∑i = 1ri(t) ― ri(0) t→∞

2

[27]

which represents the mean square displacement of the center of mass of the N guest molecules in the system.

Thermodynamic correction factor. The thermodynamic correction factor Γ is calculated by relating it to the fluctuation in the number of particles in the system, given as37-39 〈N〉

[28]

Γ = 〈N2〉 ― 〈N〉2 where N is the number of guest molecules in the system.

3. Computational Details



Page 12 of 28

Figure 1. (a) Two-dimensional projection of three-dimensional (3D) simulation lattice for sI methane hydrates. Water molecules are visualized in blue, guest molecules are in red, and hydrogen bonds are in green. The model consists of 8 × 8 × 8 hydrate unit cells, and possesses 4096 available cages for guest molecules. (b) The 3D simulation lattice for a closed system, with periodic boundary conditions in all three dimensions. (c) The 3D simulation lattice for an open system. Molecules are allowed to exchange with the outside in one of the dimensions, and periodic boundary conditions are applied in the other two dimensions.

Systems. A system of 8 × 8 × 8 unit cells of sI hydrate is used in our simulations, as illustrated in Figure 1. Such a system contains 4096 cages in total, including 3072 large and 1024 small cages. In the closed system runs, periodic boundary conditions are applied in all three dimensions, while in the open system runs, periodic boundary conditions are applied in two of the dimensions. The box length is 9.5 nm in each dimension. Models are constructed at different levels of total occupancy. 106 ~ 108 KMC simulation steps were conducted for determining the equilibrium fractional occupancies, the thermodynamic correction factors, the self-, jumpand transport-diffusivities.


Page 13 of 28 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


Rate constants. The defect-driven rate constants of guest molecules hopping through different pathways, including L6L, L5L, L5S, and S5L, are the pre-determined inputs of this study. Since the water vacancies diffuse much faster than the guest molecules20, 23, it is assumed that each hopping event is facilitated by a water vacancy present in the interface. In this work, the hopping rate constants of methane were determined based on the calculations of Peters et. al21, which were calculated at 225 K and 250 K, as summarized in Table 1.

Simulation details. Two types of KMC simulations were conducted. One was run in a closed system: the total number of particles inside the system is fixed (Fig. 1b). The other was run in an open system: particles inside the system are allowed to cross the left and right boundaries of the simulation lattice. The boundary contains one layer of a hydrate lattice with width equal to the side length of a unit cell, which is an 8 × 8 × 1 hydrate layer (Fig. 1c). Before each KMC step, the fractional occupancies of the cages on the left and right boundaries are refreshed. A random number between 0 and 1 is generated. If the random number is smaller than the specified fractional occupancy and the cage (small or large) is originally empty, a new particle will be placed in the cage. On the other hand, if the random number is larger than the specified fractional occupancy and the cage is originally occupied, the particle inside the cage will be removed and the cage will become empty. The same procedure is applied to all the cages located in the boundaries. As a result, as the simulation proceeds, the average fractional occupancies within the boundary layers will approach the specified fractional occupancies. The simulation is thus similar to a Grand Canonical Monte Carlo simulation, in which the chemical potential of particles is kept fixed. In the analysis of mean square displacement, it is convenient to record the positions of particles using a constant time increment. However, the time intervals between events in KMC varies with each step. Therefore, a time interval t for sampling positions was chosen, and calculated as the accumulated KMC simulation time divided by the number of KMC steps. The positions measured at these constant time interval increments were then determined based on the variable time interval trajectory from the KMC simulation, and used in the analysis for the mean square displacement.

The self- and jump- diffusivities were determined according to eqs. 26 and 27 respectively, by tracking the trajectories of all the particles in the closed system over 106 KMC steps. The results were obtained by ACS Paragon Plus Environment


Page 14 of 28

averaging over three independent simulation runs with different random numbers. The transport-diffusivity was calculated according to eq. 13, by measuring the molecular flux and the gradient of occupancy across the open system over 107 KMC steps. The open system KMC simulations were performed by setting different occupancies in the left and right boundary layers. The system will reach a steady state with molecular flux flowing toward the low occupancy region. Three independent simulations were conducted with different random numbers to obtain the average results of the transport-diffusivity. Both types of simulation systems (closed and open) were used to calculate the thermodynamic correction factor and to check for consistency. According to eq. 28, the thermodynamic correction factor in the closed system is measured by monitoring the fluctuation in the number of particles contained in a small probe volume randomly placed in the system. (We neglect a small finite-size correction by assuming the system outside the probe volume acts as an infinite reservoir.) For the open system, the particles can be inserted into or deleted from the simulation lattice, and we measured the fluctuation in the total number of particles within the entire simulation box (including the final layer at the boundary). 108 KMC steps in both of the systems were performed to calculate the thermodynamic correction factor. The results of closed system run were obtained by averaging over each small probe volume, and the results of open system run were calculated by averaging over every 107 KMC steps.

4. Results and Discussion 4.1 Equilibrium guest distribution in sI hydrate lattice. The equilibrium fractional occupancies for large and small cages are calculated by averaging over 106 KMC simulation steps, weighted by the length of the time increment of each step, as shown in the following equations: t

⟨XL⟩ =

∫0XL(t)dτ t

∫0dτ

t

and ⟨XS⟩ =

∫0XS(t)dτ

[29]

t

∫0dτ

Figure 2 compares the fractional occupancies XL and XS of methane hydrate determined from KMC and the analytical model (eqs. 7 and 8) at 225 and 250 K. The analytical model can perfectly describe the equilibrium guest distribution obtained from the KMC simulations under different total occupancies. This is as expected, since the analytical model gives an exact solution for the KMC model being simulated, aside from statistical noise, and is equivalent to the van der Waals-Platteeuw model and the model derived by Coppens et al.29 The kS5L

analytical model is useful for understanding that the equilibrium occupancy is only controlled by 𝐾 = kL5S. ACS Paragon Plus Environment

Page 15 of 28 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 2. KMC simulations and analytical model results of equilibrium fractional occupancies versus total occupancy at 225 K (𝐾=37.86) and 250 K (𝐾=13.46). Open symbols are the results for large (circles) and small (diamonds) cage occupancies obtained from the KMC simulations. The lines represent the results of large (solid lines) and small (dash lines) cage occupancies from the analytical model (eqs. 7 and 8). XL and XS represent the fractional occupancies of large and small cages respectively. The uncertainties from KMC simulations are smaller than the data symbols.

Most of the experimental measurements and simulation estimations were conducted in the temperature range 273 – 291 K. Sum et al. and Uchida et al. measured the occupancy of methane hydrate by Raman spectroscopy.40-41 Their measurements of fractional occupancies of large and small cages were reported under several sets of temperature and pressure conditions. In order to compare the results, experimental fractional occupancies are expressed as a function of the total occupancy here. Prior simulation estimations were conducted by ab initio calculation,42-45 Monte Carlo method46-48 and thermodynamic models.49 To estimate the fractional occupancies XL and XS by eqs. 7 and 8, we calculated the value of 𝐾 at temperatures ranging from 273 to 291 K based on the Langmuir constants (eq. 10) from the empirical relation which is valid in the temperature range 260 to 300 K according to Parrish and Prausnitz.44, 50 The calculated 𝐾 values at 273 and ACS Paragon Plus Environment


Page 16 of 28

291 K are 5.49 and 5.45, respectively. Using these values of 𝐾, we found that the fractional occupancies XL and XS are not sensitive to temperature. For XL, the differences in the temperature range (273 – 291 K) are within 0.03%, and for XS, the differences are within 0.16%. Figure 3 compares the fractional occupancies XL and XS obtained from the literature and the analytical model (eqs. 7 and 8). The results of fractional occupancies from the analytical model agree well with the experimental measurements from Sum et al.40 However, the small cage occupancy is overestimated comparing with the Raman spectrometry work of Uchida et al.41 This overestimation of small cage occupancy was also reported by Brumby et al. when the TIP4P/Ice water model is used as their simulation force field.46 It is also possible that the overestimation results from the relatively large uncertainties in the experimental measurements from Uchida et al. Comparing with the ab initio calculations, the results of fractional occupancies from the analytical model are in good agreement with the calculations from Klauda and Sandler43-44, while the small cage occupancy is lower comparing with the results from Cao et al.42, who predicted higher small cage occupancies than large cage occupancies. Except for these disparities, a fairly good agreement with the prior experimental measurements and simulation estimations was found by using the analytical model to calculate the equilibrium guest distribution. This good agreement also provides a check on (the ratio of) the hopping rate constants used here.


Page 17 of 28 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. Comparison of equilibrium fractional occupancies XL (blue) or XS (orange) versus total occupancy (Xt) in sI methane hydrate at 273 – 291 K. The symbols are experimental (squares (Sum et al.)40, diamonds (Uchida et al.)41) and simulation (triangles (Cao et al.)42, stars (Klauda & Sandler)44, asterisks (Sun & Duan)45, circles (Brumby et al.)46, dashes (Waage et al.)47, hyphens (Jensen et al.)48, crosses (Bhawangirkar et al.)49) results from the literature and the solid lines represent the results from the analytical model.

4.2 The transport properties of guest in sI hydrate lattice. Figure 4 compares the thermodynamic correction factors, jump- and transport- diffusivities from KMC simulations and the analytical model at 250 K. (Results for 225 K are given in Fig. S1 of Supporting Information.) The jump- and transport- diffusivity results in the plot are normalized by the zero-loading diffusivity, which is the diffusivity when the simulated system contains only one guest molecule (infinite dilution limit). It can be seen that the analytical model (eqs. 18, 21 and 26) agrees very well with the KMC simulation results, while the one cage type thermodynamic correction factor 1/(1-Xt) and jump-diffusivity using the factor (1-Xt) do not. For systems containing only one type of cages (e.g. many zeolites), the transport-diffusivity is independent of the total occupancy since the thermodynamic correction factor will be the reciprocal of the jump-diffusivity.51 For sI hydrate, the transportdiffusivity decreases with the total occupancy. The decrease becomes significant especially when the total occupancy is higher than 0.6. For example, the transport-diffusivity for total occupancy higher than 95% is approximately half of the zero-loading diffusivity. Comparing with the transport-diffusivity, the jumpdiffusivity shows a stronger dependence on the total occupancy. When the total occupancy is higher than 95%, the jump-diffusivity is about two orders of magnitude smaller than the transport-diffusivity.



Page 18 of 28

Figure 4. Comparison of the thermodynamic correction factor, normalized jump- and transport- diffusivities versus total occupancy for KMC simulations and the analytical model at 250 K. The symbols are KMC simulation results, the solid lines represent the results from analytical model, and the dashed lines represent the one cage type relation. The uncertainties are smaller than the data symbols.

4.3 Contribution of hopping events from the four different pathways. The analytical model can be used to analyze the distribution of hopping events from the four pathways (L6L, L5L, L5S, and S5L). According to eq. 21, the fractional contribution to the jump-diffusivity through a specific kind of pathway is given by: kiPi ∑𝑖kiPi

[30]

(with i = L6L, L5L, L5S, S5L).

In the KMC simulation, the contribution can be determined by the number of hopping events along a specific pathway divided by the total number of hopping events. The results from the analytical model and the KMC simulations at 250 K are shown in Figure 5. When the occupancy is low, the hopping events are dominated by the hopping through the L6L pathway. As the occupancy becomes higher, the contributions of hopping events through the L5S and S5L pathways gradually increase, while the contribution of hopping events through the L5L pathway becomes even lower. The hopping through L6L is more dominant at a lower ACS Paragon Plus Environment

Page 19 of 28 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


temperature (See Fig. S2 in Supporting Information for results at 225 K.) Note that the contributions from the L5S and S5L pathways are equal at equilibrium to satisfy detailed balance. According to eq. 16, the transportdiffusivity is the product of the jump-diffusivity and the thermodynamic correction factor. Since the thermodynamic correction factor is a constant when the total occupancy is fixed, the contributions of the four pathways to the transport-diffusivity are identical to their contributions to the jump-diffusivity.

Figure 5. Distribution of the hopping events through the four pathways (L6L, L5L, L5S, and S5L) versus the total occupancy at 250 K. The symbols are the results from the KMC simulation and the lines represent the results from the analytical model.

4.4 Self-diffusivity. While an exact expression of the self-exchange coeffficent Di,i is not available, we follow the suggestion of Krishna and Paschek32 to assume that the self-exchange coeffficent is equal to the jump-diffusivity, i.e. Di,i = Djump(Xt). Figure 6 compares the self- and jump- diffusivities obtained from the KMC simulations and the analytical model (eqs. 21 and 23) at 250 K. (Results for 225 K are given in Fig. S3 in Supporting Information.) As expected, the self-diffusivity is lower than the jump-diffusivity except for the zero-loading case (one guest molecule in the system) due to correlation effects. The analytical model can give ACS Paragon Plus Environment


Page 20 of 28

a fairly good estimation of the self-diffusivity in comparing with the KMC simulation results.

Figure 6. KMC simulations and analytical model results of self- and jump- diffusivities versus total occupancy at 250 K. The hollow symbols are the KMC simulation results and the solid lines represent the results from the analytical model.

4.5 Prediction of methane transport diffusivity in sI hydrate under experimental conditions. According to the validations that have been made in the previous sections, the analytical model can properly describe the KMC simulation results of the thermodynamic and transport properties at low temperatures (225 K and 250 K). Here we compare the results of jump- and transport- diffusivities determined by the analytical model with existing experimental studies at higher temperatures. We use the Arrhenius law to extrapolate the jump- and transport- diffusivity to different temperatures EA

[31]

D = D0 × exp ( ― RT)

The total occupancies are set to be 99% at 225 K and 98% at 250 K according to eq. 12, with the Langmuir constants obtained according to Klauda and Sandler44, 50 and the gas fugacity equal to 5.09 MPa according to ACS Paragon Plus Environment

Page 21 of 28 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 experimental methane fugacity by Salamatin et al.19 Fitting the jump-diffusivity under the two conditions mentioned above to the Arrhenius form, we obtain that D0 is 3.41×10-4 (m2/s) and the activation energy EA is 59.5 (kJ/mol). The energy barrier is in very good agreement with prior experiment measurements of the permeation coefficients of methane during the hydrate formation from ice powder, ranging between 52 and 62 (kJ/mol).1 As also shown in Figure 7, the predictions of jump-diffusivity are in good agreement with the experimental measurements, with an average absolute relative deviation of 63%.

In their recently developed theory of the “hole-in-cage-wall” diffusion mechanism,18 Salamatin et al. predicted a non-trivial relation between the diffusion (Dg) and permeation coefficients (Lg) in the pure sI-hydrate formation experiments even for a regular case of quasi-stationary gas transport with a local equilibrium guest distribution. The experimental data of permeation coefficients1 (Kuhs et al. 2006) were converted to the diffusion coefficients19 (Salamatin et al. 2017) via the relation given from their theoretical model.18-19 The permeation coefficient and the diffusion coefficient in their model descriptions correspond to the jump- and transport- diffusivity, respectively, in this work. The predicted transport-diffusivities from eq. 31 are also reported to compare with the experiment measurements of the diffusivity of methane measured from the gas replacement process.19 Fitting the transport-diffusivity under the two conditions mentioned above to the Arrhenius form, we obtained that D0 is 3.76×10-5 (m2/s) and the activation energy EA is 46.7 (kJ/mol). The results of the analytical model agree well with the most recent experimental measurements and model descriptions by Salamatin et al. (Figure 7), with an average absolute relative deviation of 25%.

Even though we only have the hopping rate constants at 225 and 250 K,21 the success of the extrapolation of diffusivities allowed us to estimate the transport-diffusivity of methane in hydrates over a wide range of temperature and different conditions of occupancy. This application should be useful for modeling the gas replacement process, of which the system temperature and occupancy conditions may change with time. In addition, it should also be feasible to apply the analytical model to other kinds of guest molecules such as CO2 and N2 in sI hydrates, as long as the hopping rate constants of each pathway are available.

While the agreement with experiments is impressive, it is useful to reiterate the simplifications and ACS Paragon Plus Environment


Page 22 of 28

assumptions made in the analytical model. First of all, we assume no interactions between the gas molecules (other than not being able to simultaneously occupy the same cage as another gas molecule) and so the guest hopping events are assumed to be statistically uncorrelated. We also make mean-field assumptions for the surrounding vacancy concentration of neighboring cages when computing hopping probabilities. Furthermore, perfect sI hydrate structure with water vacancies is assumed, while in reality, there could be many types of defects other than the water vacancy (e.g., Bjerrum defect, grain boundary), which might result in a higher diffusion rate.52-54 Perhaps the most remarkable fact is that the ratio of rate constants, 𝐾 = kS5L/kL5S, obtained from Peters et al. significantly deviates from the experimental measurements of the ratio of Langmuir constant (CL/CS). For example, the 𝐾 value obtained from Peters et al. is 37.86 at 225 K and 13.46 at 250 K, while from experiment the value is temperature insensitive, about 5.9 (See Fig. S4 in Supporting Information for further details). To perform a sensitivity analysis on the effect of these rate constants, we varied the value of 𝐾 from 5 to 45 at 225 K by fixing either kS5L or kL5S in the analytical model. We found that the jump-diffusivity increases by at most 80% when fixing kS5L and increases by at most 56% when fixing kL5S as 𝐾 decreases from 45 to 5. At 250 K, the jump-diffusivity increases by at most 30% when fixing kS5L and increases by at most 10% when fixing kL5S as 𝐾 decreases from 15 to 5. Given these uncertainties, we conducted the temperature extrapolations of jump-diffusivities based on the maximum deviations (80% at 225 K and 30% at 250 K). To obtain the upper and lower bounds of the results, we made the extrapolations with the combination of the largest Djump at 225 K and the lowest Djump at 250 K (or the lowest Djump at 225 K and the largest Djump at 250 K). The extrapolated jump-diffusivities compare reasonably well with experimental values (Fig. S5 in Supporting Information). The results imply that the ratio 𝐾 of rate constants does not directly affect the prediction of jump-diffusivity given that the absolute value of rate constants are within a reasonable range. In addition, the thermodynamic correction factor in the high occupancy situation shows little dependence on the ratio 𝐾 of rate constants. There is only an approximately 1% difference in the calculated thermodynamic correction factor when changing the ratio 𝐾 of rate constants from 5.9 to 37.86. Therefore, this error can be ignored when predicting the transport-diffusivity.


Page 23 of 28 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 7. Comparison of the extrapolation results of jump- (orange) and transport- (blue) diffusivities with the reported experiment measurements. The symbols are the literature results1, 19 and the solid lines represent the results from the analytical model.

5. Conclusion The transport-diffusivity of methane in sI hydrate is an important property for understanding the solid state replacement process and equilibration of methane hydrate. In this work, we derived a mathematical expression for transport-diffusivity in sI hydrate according to the Maxwell-Stefan diffusion theory, which relates the transport-diffusivity to the hopping rate constants of guest molecules and the total occupancy. All of the results calculated by our analytical model are in excellent agreement with the corresponding KMC simulation results, including equilibrium occupancies, thermodynamic correction factors, and diffusivities (self-, jump- and transport-diffusivities). According to the analytical model, the thermodynamic properties (equilibrium distribution, thermodynamic correction factor) are only affected by the ratio K = kS5L/kL5S of the hopping rate constants, while the transport properties will depend on all the hopping rate constants of the possible pathways. Comparing to the most recent experiment measurements of methane diffusivity in hydrates, our results for transport-diffusivity agrees well with experimental estimations over a wide range of temperatures. The ACS Paragon Plus Environment


Page 24 of 28

developed analytical model would be useful for modeling the solid-state guest exchange process, which is of great importance in the replacement of methane with CO2 in natural gas hydrate reservoirs.

Supporting Information Details for the derivation of eqs. 7,8, 18, and 19; The results from the analytical model and the KMC simulations at 225 K; comparison of the ratio of Langmuir constants (CL/CS) from the literature; comparison of the extrapolation results (upper and lower bounds) of the jump-diffusivities with the experimental measurements

Acknowledgement This research was supported in part by the Ministry of Science and Technology of Taiwan (MOST 104-2221E-002-186-MY3, 106-3113-M-002-001, 107-3113-M-002-001) and National Taiwan University (NTU-CDP106L7827). David Wu would like to acknowledge the Ministry of Science and Technology of Taiwan (MOST 106-2811-E-002-020) for supporting his visit to NTU from 2017 to 2018. The computational resources from the National Center for High-Performance Computing of Taiwan and the Computing and Information Networking Center of the National Taiwan University are acknowledged.


Page 25 of 28 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


Reference 1. Kuhs, W. F.; Staykova, D. K.; Salamatin, A. N., Formation of Methane Hydrate from Polydisperse Ice Powders. Journal of Physical Chemistry B 2006, 110, 13283-13295. 2. Komatsu, H.; Ota, M.; Smith, R. L.; Inomata, H., Review of Co2-Ch4 Clathrate Hydrate Replacement Reaction Laboratory Studies - Properties and Kinetics. Journal of the Taiwan Institute of Chemical Engineers 2013, 44, 517-537. 3. Zhao, J. F.; Xu, K.; Song, Y. C.; Liu, W. G.; Lam, W.; Liu, Y.; Xue, K. H.; Zhu, Y. M.; Yu, X. C.; Li, Q. P., A Review on Research on Replacement of Ch4 in Natural Gas Hydrates by Use of Co2. Energies 2012, 5, 399419. 4. Baldwin, B. A.; Stevens, J.; Howard, J. J.; Graue, A.; Kvamme, B.; Aspenes, E.; Ersland, G.; Husebo, J.; Zornes, D. R., Using Magnetic Resonance Imaging to Monitor Ch4 Hydrate Formation and Spontaneous Conversion of Ch4 Hydrate to Co2 Hydrate in Porous Media. Magnetic Resonance Imaging 2009, 27, 720726. 5. Kvamme, B.; Graue, A.; Buanes, T.; Kumetsoua, T.; Ersland, G., Storage of Co2 in Natural Gas Hydrate Reservoirs and the Effect of Hydrate as an Extra Sealing in Cold Aquifers. International Journal of Greenhouse Gas Control 2007, 1, 236-246. 6. Yuan, Q.; Sun, C. Y.; Yang, X.; Ma, P. C.; Ma, Z. W.; Liu, B.; Ma, Q. L.; Yang, L. Y.; Chen, G. J., Recovery of Methane from Hydrate Reservoir with Gaseous Carbon Dioxide Using a Three-Dimensional Middle-Size Reactor. Energy 2012, 40, 47-58. 7. Jung, J. W.; Espinoza, D. N.; Santamarina, J. C., Properties and Phenomena Relevant to Ch4-Co2 Replacement in Hydrate-Bearing Sediments. Journal of Geophysical Research: Solid Earth (1978–2012) 2010, 115. 8. Ors, O.; Sinayuc, C., An Experimental Study on the Co2-Ch4 Swap Process between Gaseous Co2 and Ch4 Hydrate in Porous Media. Journal of Petroleum Science and Engineering 2014, 119, 156-162. 9. Zhou, X. B.; Liang, D. Q.; Liang, S.; Yi, L. Z.; Lin, F. H., Recovering Ch4 from Natural Gas Hydrates with the Injection of Co2-N-2 Gas Mixtures. Energy & Fuels 2015, 29, 1099-1106. 10. Lee, Y.; Kim, Y.; Lee, J.; Lee, H.; Seo, Y., Ch4 Recovery and Co2 Sequestration Using Flue Gas in Natural Gas Hydrates as Revealed by a Micro-Differential Scanning Calorimeter. Applied Energy 2015, 150, 120-127. 11. Ota, M.; Abe, Y.; Watanabe, M.; Smith, R. L.; Inomata, H., Methane Recovery from Methane Hydrate Using Pressurized Co2. Fluid Phase Equilib. 2005, 228-229, 553-559. 12. Ota, M.; Morohashi, K.; Abe, Y.; Watanabe, M.; Smith, J. R. L.; Inomata, H., Replacement of Ch4 in the Hydrate by Use of Liquid Co2. Energy Conversion and Management 2005, 46, 1680-1691. 13. Zhou, X.; Fan, S.; Liang, D.; Du, J., Replacement of Methane from Quartz Sand-Bearing Hydrate with Carbon Dioxide-in-Water Emulsion. Energy Fuels 2008, 22, 1759-1764. 14. Jung, J. W.; Espinoza, D. N.; Santamarina, J. C., Properties and Phenomena Relevant to Ch4-Co2 Replacement in Hydrate-Bearing Sediments. Journal of Geophysical Research: Solid Earth 2010, 115. 15. Boswell, R.; Schoderbek, D.; Collett, T. S.; Ohtsuki, S.; White, M.; Anderson, B. J., The Inik Sikumi Field Experiment, Alaska North Slope: Design, Operations, and Implications for Co2-Ch4 Exchange in Gas Hydrate Reservoirs. Energy & Fuels 2017, 31, 140-153. 16. Schoderbek, D.; Boswell, R., Ignik Sikumi #1, Gas Hydrate Test Well, Successfully Installed on the Alaska North Slope. Fire ice 2011, 1-5. 17. Falenty, A.; Salamatin, A. N.; Kuhs, W. F., Kinetics of Co2-Hydrate Formation from Ice Powders: Data ACS Paragon Plus Environment


Page 26 of 28

Summary and Modeling Extended to Low Temperatures. Journal of Physical Chemistry C 2013, 117, 84438457. 18. Salamatin, A. N.; Falenty, A.; Hansen, T. C.; Kuhs, W. F., Guest Migration Revealed in Co2 Clathrate Hydrates. Energy & Fuels 2015, 29, 5681-5691. 19. Salamatin, A. N.; Falenty, A.; Kuhs, W. F., Diffusion Model for Gas Replacement in an Isostructural Ch4Co2 Hydrate System. Journal of Physical Chemistry C 2017, 121, 17603-17616. 20. Demurov, A.; Radhakrishnan, R.; Trout, B. L., Computations of Diffusivities in Ice and Co2 Clathrate Hydrates Via Molecular Dynamics and Monte Carlo Simulations. Journal of Chemical Physics 2002, 116, 702709. 21. Peters, B.; Zimmermann, N. E. R.; Beckham, G. T.; Tester, J. W.; Trout, B. L., Path Sampling Calculation of Methane Diffusivity in Natural Gas Hydrates from a Water-Vacancy Assisted Mechanism. Journal of the American Chemical Society 2008, 130, 17342-17350. 22. Lo, H.; Lee, M. T.; Lin, S. T., Water Vacancy Driven Diffusion in Clathrate Hydrates: Molecular Dynamics Simulation Study. Journal of Physical Chemistry C 2017, 121, 8280-8289. 23. Liang, S.; Kusalik, P. G., The Mobility of Water Molecules through Gas Hydrates. Journal of the American Chemical Society 2011, 133, 1870-1876. 24. Liang, S.; Liang, D. Q.; Wu, N. Y.; Yi, L. Z.; Hu, G. W., Molecular Mechanisms of Gas Diffusion in Co2 Hydrates. Journal of Physical Chemistry C 2016, 120, 16298-16304. 25. Waage, M. H.; Trinh, T. T.; Erp, T. S. v., Diffusion of Gas Mixtures in the Si Hydrate Structure. The Journal of Chemical Physics 2018, 148, 214701. 26. Krishna, R.; Wesselingh, J. A., Review Article Number 50 - the Maxwell-Stefan Approach to Mass Transfer. Chem. Eng. Sci. 1997, 52, 861-911. 27. Krishna, R., Describing the Diffusion of Guest Molecules inside Porous Structures. Journal of Physical Chemistry C 2009, 113, 19756-19781. 28. Karger, J.; Ruthven, D. M., Diffusion in Nanoporous Materials: Fundamental Principles, Insights and Challenges. New Journal of Chemistry 2016, 40, 4027-4048. 29. Coppens, M.-O.; Bell, A. T.; Chakraborty, A. K., Dynamic Monte-Carlo and Mean-Field Study of the Effect of Strong Adsorption Sites on Self-Diffusion in Zeolites. Chemical Engineering Science 1999, 54, 3455-3463. 30. Coppens, M. O.; Iyengar, V., Testing the Consistency of the Maxwell-Stefan Formulation When Predicting Self-Diffusion in Zeolites with Strong Adsorption Sites. Nanotechnology 2005, 16, S442-S448. 31. Waals, J. H. v. d.; Platteeuw, J. C., Clathrate Solutions. In Advances in Chemical Physics, 1958. 32. Paschek, D.; Krishna, R., Inter-Relation between Self- and Jump-Diffusivities in Zeolites. Chemical Physics Letters 2001, 333, 278-284. 33. Skoulidas, A. I.; Sholl, D. S.; Krishna, R., Correlation Effects in Diffusion of Ch4/Cf4 Mixtures in Mfi Zeolite. A Study Linking Md Simulations with the Maxwell-Stefan Formulation. Langmuir 2003, 19, 7977-7988. 34. Bortz, A. B.; Kalos, M. H.; Lebowitz, J. L., New Algorithm for Monte-Carlo Simulation of Ising Spin Systems. Journal of Computational Physics 1975, 17, 10-18. 35. Voter, A. F., Introduction to the Kinetic Monte Carlo Method; Springer, 2007; Vol. 235. 36. Battaile, C. C., The Kinetic Monte Carlo Method: Foundation, Implementation, and Application. Computer Methods in Applied Mechanics and Engineering 2008, 197, 3386-3398. 37. Reed, D. A.; Ehrlich, G., Surface-Diffusion, Atomic Jump Rates and Thermodynamics. Surface Science 1981, 102, 588-609. ACS Paragon Plus Environment

Page 27 of 28 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


38. Uebing, C.; Pereyra, V.; Zgrablich, G., Diffusion of Interacting Lattice Gases on Heterogeneous Surfaces with Simple Topographies. Surface Science 1996, 366, 185-192. 39. Viljoen, E. C.; Uebing, C., Effects of Heterogeneity on the Surface Diffusion of Interacting Lattice Gases: The Bivariate Trap Model. Surface Science 1996, 352, 1007-1011. 40. Sum, A. K.; Burruss, R. C.; Sloan, E. D., Measurement of Clathrate Hydrates Via Raman Spectroscopy. Journal of Physical Chemistry B 1997, 101, 7371-7377. 41. Uchida, T.; Hirano, T.; Ebinuma, T.; Narita, H.; Gohara, K.; Mae, S.; Matsumoto, R., Raman Spectroscopic Determination of Hydration Number of Methane Hydrates. Aiche Journal 1999, 45, 2641-2645. 42. Cao, Z.; Tester, J. W.; Sparks, K. A.; Trout, B. L., Molecular Computations Using Robust Hydrocarbon−Water Potentials for Predicting Gas Hydrate Phase Equilibria. The Journal of Physical Chemistry B 2001, 105, 10950-10960. 43. Klauda, J. B.; Sandler, S. I., Ab Initio Intermolecular Potentials for Gas Hydrates and Their Predictions. The Journal of Physical Chemistry B 2002, 106, 5722-5732. 44. Klauda, J. B.; Sandler, S. I., Phase Behavior of Clathrate Hydrates: A Model for Single and Multiple Gas Component Hydrates. Chem. Eng. Sci. 2003, 58, 27-41. 45. Sun, R.; Duan, Z., Prediction of Ch4 and Co2 Hydrate Phase Equilibrium and Cage Occupancy from Ab Initio Intermolecular Potentials. Geochimica et Cosmochimica Acta 2005, 69, 4411-4424. 46. Brumby, P. E.; Yuhara, D.; Wu, D. T.; Sum, A. K.; Yasuoka, K., Cage Occupancy of Methane Hydrates from Gibbs Ensemble Monte Carlo Simulations. Fluid Phase Equilibria 2016, 413, 242-248. 47. Waage, M. H.; Vlugt, T. J. H.; Kjelstrup, S., Phase Diagram of Methane and Carbon Dioxide Hydrates Computed by Monte Carlo Simulations. Journal of Physical Chemistry B 2017, 121, 7336-7350. 48. Jensen, L.; Thomsen, K.; von Solms, N.; Wierzchowski, S.; Walsh, M. R.; Koh, C. A.; Sloan, E. D.; Wu, D. T.; Sum, A. K., Calculation of Liquid Water−Hydrate−Methane Vapor Phase Equilibria from Molecular Simulations. The Journal of Physical Chemistry B 2010, 114, 5775-5782. 49. Bhawangirkar, D. R.; Adhikari, J.; Sangwai, J. S., Thermodynamic Modeling of Phase Equilibria of Clathrate Hydrates Formed from Ch4, Co2, C2h6, N-2 and C3h8, with Different Equations of State. Journal of Chemical Thermodynamics 2018, 117, 180-192. 50. Parrish, W. R.; Prausnitz, J. M., Dissociation Pressures of Gas Hydrates Formed by Gas Mixtures. Industrial & Engineering Chemistry Process Design and Development 1972, 11, 26-35. 51. Krishna, R.; Paschek, D.; Baur, R., Modeling the Occupancy Dependence of Diffusivities in Zeolites. Microporous and Mesoporous Materials 2004, 76, 233-246. 52. Bjerrum, N., Structure and Properties of Ice. Science 1952, 115, 385-390. 53. Kuhs, W. F.; Klapproth, A.; Gotthardt, F.; Techmer, K.; Heinrichs, T., The Formation of Meso- and Macroporous Gas Hydrates. Geophysical Research Letters 2000, 27, 2929-2932. 54. Stern, L. A.; Circone, S.; Kirby, S. H.; Durham, W. B., Anomalous Preservation of Pure Methane Hydrate at 1 Atm. The Journal of Physical Chemistry B 2001, 105, 1756-1762.



TOC Graphics


Page 28 of 28

Theory and Kinetic Monte Carlo Simulation of Guest Molecule

Recommend Documents