Mechanisms and Nucleation Rate of Methane Hydrate by Dynamical

Subscriber access provided by UNIV NEW ORLEANS

Article

Mechanisms and Nucleation Rate of Methane-Hydrate by Dynamical Nonequilibrium Molecular Dynamics Marco Lauricella, Giovanni Ciccotti, Niall Joseph English, Baron Peters, and Simone Meloni J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.7b05754 • Publication Date (Web): 06 Oct 2017 Downloaded from http://pubs.acs.org on October 9, 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 C 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 55

The Journal of Physical Chemistry 1

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

Mechanisms and Nucleation Rate of Methane-Hydrate by Dynamical Nonequilibrium Molecular Dynamics. Marco Lauricella a,b), Giovanni Ciccotti a,b,c), Niall J. English b),, Baron Peters e), Simone Meloni f,b) a) Istituto per le Applicazioni del Calcolo, Consiglio Nazionale delle Ricerche, Rome, Italy b) School of Physics, University College Dublin, Belfield, Dublin 4, Ireland c) Dipartimento di Fisica, Università La Sapienza, Rome, Italy d) School of Chemical and Bioprocess Engineering, University College Dublin, Belfield, Dublin 4, Ireland. e) Dept. of Chemical Engineering, University of California, Santa Barbara, California 93106, United States f) Dipartimento di Ingegneria Meccanica e Aerospaziale, Università La Sapienza, Rome, Italy PACS: 64.60.qe

ABSTRACT: We investigate the effects of high solvated-methane concentration on methane-hydrate nucleation at 250 K and 500 atm. We consider solutions at four levels of methane molar fraction in the initial H2O-CH4 solution, = 0.038, 0.044,

0.052 and 0.058, which are higher than (metastable) bulk supersaturation. is controlled independently of the temperature and pressure thanks to the use of special simulation techniques [Phys. Chem. Chem. Phys. 13, 13177 (2011)].

These

conditions mimic a possible increase of local methane concentration beyond supersaturation induced, for example, by freeze-concentration or thermal fluctuations. The nucleation mechanism and kinetics are investigated using the dynamical approach to non-equilibrium molecular dynamics. We demonstrate a hydrate- forming/ordering process of solvated methane and water molecules in a manner consistent with both the ‘blob’ hypothesis and ‘cage adsorption hypothesis’ hypotheses: the system initially forms an amorphous nucleus at high methane concentration, which then gets ordered forming the clathrate crystallite. We evaluate nucleation rates using both the methods of the mean first-passage time, i.e. the curve of the average time the system takes to reach a crystalline nucleus of given size, and survival-probability, i.e. probability that up to a given time the system has not nucleated yet. We found a dependence of the nucleation rate on initial methane concentration of a form consistent with the dependence of classical nucleation theory rate on supersaturation, and determined the

ACS Paragon Plus Environment

The Journal of Physical Chemistry

Page 2 of 55 2

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

relevant parameters of this relation. We found a very rapid increase of nucleation rate with solvated-methane concentration, proving that methane molar fraction, even beyond bulk supersaturation, is key at triggering the homogenous nucleation of clathrate. We derive a kinetic equation that allows for estimation of the nucleation rate over a wide range of concentration conditions.


Page 3 of 55


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

1. Introduction Clathrate hydrates are non-stoichiometric crystalline inclusion compounds in which a water host lattice encages small guest atoms or molecules in cages; the empty hydrogen bonded lattice is thermodynamically unstable, and its existence is due to stabilisation resulting from the enclathration of the trapped solutes in its cages.1,2 There are three known common hydrate structures: sI, sII and sH. In type I hydrate, the unit cell is formed from two small 512 (the base indicates the type of polygon forming the face of the polyhedron – 5 stands for pentagon – and the exponent the number of such faces) pentagonal dodecahedral cages and six slightly larger tetrakaidecahedral 51262 cages, with 46 water molecules.1,2 In type II hydrate, each unit cell contains 136 water molecules and 24 cages, 16 of which are small 512 cages and eight of which are medium-sized hexadecahedral 51264, while sH unit cells are composed of 34 water molecules with three small 512 cages, two small 435663 cages and one large 51268 cage.1,2 Methane hydrates are the most widespread type of clathrate hydrates, and are thought to exist in nature primarily as type I in the permafrost and deep ocean regions. With methane compositions of around 85-95 % of the maximum theoretical occupation, ~ 0.15 molar fraction of methane, methane clathrate hydrates constitute a possible significant energy resource.3,4 The large quantity of available methane hydrates has driven significant interest in their use as a possible resource of methanegas fuel,5-7 with the relative success of the several production trial (JOGMEC, Malik, SUGAR - no drilling required8) Another potential application of methane clathrate is as medium for methane-fuelstorage and transport, which is possible at 273.15 K and 26 atm9,10 while conventional Liquefied Natural Gas transport methods require temperature lower than 110 K at 20 atm with a very high energy cost of storage. ACS Paragon Plus Environment


Page 4 of 55 4

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

Another aspect of fossil fuels industry concerning clathrates is their formation in pipelines,1,2,11 which results in the plugging transmission lines and expensive production stoppages of the order of several months.12 An extreme case of damage resulting from undesired clathrate formation happened during the attempt to collect oil leaking into the Gulf of Mexico with giant dome after the Deepwater Horizon disaster, which failed because it became clogged with methane hydrates. Despite the important application- and engineering-oriented relevance of clathrates,58

the fundamental mechanisms of hydrate formation, which are key in the applications

mentioned above, are understood rather poorly, especially at the molecular level. For example, theoretical predictions of rates of nucleation often differ from experimentally measured rates by at least an order of magnitude, and often more.2 Improving our understanding of the underlying mechanisms of hydrate kinetics is desirable in terms of enhancing the viability of large-scale technological applications of methane hydrates.3-8 For example, an improved understanding of hydratenucleation mechanisms would facilitate the design process of new promoters for hydrate formation in gas-transport and storage applications, or inhibitors for its formation in gas pipelines. 13 In the literature, four principal mechanisms have been proposed for hydrate nucleation. The earliest is the ‘labile cluster hypothesis’ (LCH),14 which proposes that unstable entities composed of a guest molecule surrounded by a certain number of water molecules in their first coordination shells ordered forming cages of the kind discussed above. These entities diffuse in the liquid, and when two or more of them encounter can form larger stable clusters which grow into the hydrate crystal. However, Guo et al.15,16 later proved that labile clusters survive for only picoseconds before breaking apart. As consequence, the probability of two entities merging into a larger cluster is negligible and the proposed mechanism is not likely. ACS Paragon Plus Environment

Page 5 of 55


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

Secondly, Radhakrishnan and Trout have proposed the ‘local structuring’ hypothesis (LSH),17 wherein a thermal fluctuation causes a set of methane molecules to be arranged in a clathate-like configuration. If the number of guest molecules in the local ordered arrangement exceeds that in a critical nucleus the system grows into a clathrate crystallite. Thirdly, Guo et al proposed the ‘cage adsorption hypothesis’ (CAH).18-20 Solvated methane molecules are surrounded by water molecules that, usually, do not have an ordered structure, like cages in sI or sII clathrtes. However, ordered cages can be randomly formed; in this case, there is a minimum of the free energy for a CH4 molecule at a distance of 6.2 Å, consistent with the distance between nearest neighbor guest molecules in clathrates, and another at 10.2 Å separated by a potential barrier at 8.8 Å. The strong attraction between a water cage and methane molecules may act as the driving force of hydrate nucleation. Fourthly, Jacobson et al proposed a refined version of the local structuring hypothesis, dubbed the ‘blob hypothesis’ (BH).21 According to this, stochastic local density fluctuations of guest particles in water favour the formation of a cluster, also called blob. Blobs are amorphous clusters of solvent-separated guests at high concentration of solvated methane close to the expected value in clathrate phase. The blobs are characterised by the absence of ordered structures for both host and guest molecules, in contrast to the LCH where the guest molecules are instead in configuration consistent with the final crystal structure. If the blob size is larger than a given critical value, the blob does not vanish and it can undergo an ordering process providing the final crystal structure. English et al. studied the size of the watermethane interface at realistic solvated-methane concentration, finding similar earlystage nucleation mechanisms to the blob hypothesis.22



Page 6 of 55 6

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 the BH, it is proposed that the local high density of solvated methanes in the blob triggers nucleation rather than the order of water molecules and methane molecules, like in the LCH and LSH. Hawtin et al.,23 Walsh et al.,24 and Jacobson et al.25 have carried out MD simulations of water-methane interfaces to study methane-hydrate nucleation. Indeed, these studies have shown that hydrate nucleates initially into a phases not fully consistent with the common bulk crystal structures, but containing a variety of types of cages23,24 as well as ‘amorphous’ or dis-ordered units,24,25 somewhat consistent with the BH. Knott et al.68 predicted that clathrate homogeneous nucleation rate at equilibrium

methane concentration in solution is extremely low. This finding suggests that nucleation of methane hydrates under this condition cannot occur by a homogeneous mechanism. However, since methane solubility in water is several orders of magnitude lower than the methane concentration in clathrate phase at normal conditions, intuition suggests that high CH4 concentrations promote methane clathrate nucleation, and this step is crucial in the crystallization process. Indeed, consistent with this hypothesis, Walsh et al.24 in their brute force MD investigation made the key observation that nucleation occurs when stochastic fluctuations increase the concentration of methane in solution well beyond the saturation value. High local methane concentration is also expected to enhance nucleation within CAH as it increases the number of molecules at a distance lower than the barrier from the guest methane of the relevant hydrated unit, 8.8 Å. In particular, Guo and Rodger predicted that within the CAH a ~ 0.05 critical local mole fraction exist beyond which nucleation is, essentially, barrierless.26 Poon and Peters40 have developed a stochastic model of freezing-assisted nucleation, in which a boundary layer of giant supersaturation develops ahead of the moving ice


Page 7 of 55


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

front

enhancing

clathrate

nucleation

and

possibly

allowing

homogeneous

crystallization on the experimental timescale. Consistent with the conclusions of Walsh et al.24 and Poon and Peters,40 on the basis of their metadynamics simulations Lauricella et al.27 observed that the nucleation mechanism consists of four steps: firstly the solution is enriched in methane beyond the equilibrium value; secondly, solvated methane molecules meet to form an unstructured cluster; thirdly, the excess water in the cluster is expelled, and an amorphous nucleus of clathrate composition is formed. Finally, methane and water molecules of the amorphous nucleus undergo a disorder-order transition to form the clathrate crystallite. This extended CAH or BH mechanism is characterised by a free-energy barrier of several hundreds of kBT, the thermal energy available to the system. Other biased MD methods have been applied recently to methane-hydrate nucleation. Surupria and Debenedetti have generated an ensemble of initial homogeneous conditions at composition higher than the equilibrium concentration starting from a melt of CO2 water clathrate, which they then let to evolve by standard MD. 28 They observed the formation of subcritical water-methane clusters that over a time scale of several hundred nanoseconds aggregate to form critical hydrate nuclei. These clusters, which undergo continuous structural rearrangements, contain cages typical of sI clathrates, such as 512 and 51262 cages, as well as other uncommon 51263 and 51264 cages, but do not possess long-range order. Małolepsza and Keyes applied an enhanced version of the generalised replica exchange method (gREM)29 to methane-hydrate nucleation. They found that initially methane is trapped in the normal, hexagonal, or cubic ice.30 Then, once a small amount of clathrate hydrate is formed, water molecules rearrange to form empty cages; this transforms the remainder of the system into metastable β ice, which is hypothesised to be an intermediate in hydrate formation for hosting more guest ACS Paragon Plus Environment


Page 8 of 55 8

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

molecules when available. This mechanism contrasts with the BH, which has received a growing consensus over the last few years. Bi and Li31 applied the Forward-Flux Sampling (FFS) method32 to methane-hydrate nucleation. FFS allows for computing the transition rate and the ensemble of reactive trajectories. Bi and Li have shown that hydrate formation is initiated by nucleation events occurring in the vicinity of water/methane interface, and proceeds via a gradual transition from amorphous to crystalline structure, supporting the BH mechanism. In a following study, Bi et al.33 have shown that, although most of the reactive trajectories are consistent with the two-step mechanism, there exist also nucleation events in which hydrate crystallites are formed directly, without going through the amorphous stage. This short survey of recent articles based on molecular simulation underscores that nucleation-mechanism results are not completely consistent. Considering the difficulty in the direct comparison with experiments, and the very limited statistics which can be obtained with brute-force MD, accelerated simulations with minimal biasing, i.e., employing methods which do not enforce a specific nucleation mechanism, might help to progress the field of elucidating nucleation. In this sense, it is worth stressing that simulations show that nucleation is favoured by supersaturation and/or fluctuations of methane concentration above the (stable or metastable) bulk value. For example, Moon et al. find that it takes place at the water/methane interface, where methane concentration is expected to be higher than that given by equilibrium solubility (see also Fig. 3 of the present work).40 We refer to supersaturation and concentration fluctuations which are two different quantities, the first refers to the equilibrium methane concentration at the pressure and temperature of the initial metastable phase, the second to the present non-equilibrium concentration which can be the result of thermal fluctuations or other processes, e.g. crystallization of the solution. ACS Paragon Plus Environment

Page 9 of 55


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 general, the dependence (increase) of nucleation rate on solute concentration, in particular on supersaturation, is predicted by several models based on extensions of the standard single component classical nucleation theory (see, e.g., Kashchiev et al.35 and Peters),36 and is confirmed experimentally.37-41 This ‘backdrop’ provides a twofold motivation for the present work: i) to identify the nucleation mechanism and its dependence on (non-equilibrium) concentrations using a simulation technique with minimal biasing, and ii) to establish a relationship between rate constant and methane concentration for methane-hydrate nucleation. Both objectives are difficult to achieve experimentally. The determination of a rate vs methane concentration law will allow, in turn, to identify the local methane molar fraction triggering nucleation on the experimental timescale. We have achieved this aim by using the dynamical approach to non-equilibrium molecular dynamics (D-NEMD) with macroscopic initial conditions corresponding to a controlled degree of supersaturation. Nucleation rates are quantified using the mean first-passage time42 and survival-probability approaches:43 in the first case, one measures the average time the system takes to produce a clathrate nucleus of prescribed size, from which one can obtain the nucleation rate; in the second case, one measures the probability density that no nucleation event is occurred within a given time. Also in this case, from the distribution one can obtain the nucleation rate. Yuhara et al. have used similar approaches in their non-equilibrium simulation of nucleation, but without independently controlling the methane concentration in solution.44 Their results are in fair agreement with free-energy calculations at the same thermodynamic conditions.45 This article is organised as follows: in section 2, we describe the computational methods used in this work. In section 3, we report the simulation details. In section 4 we present and discuss the results of our simulations. Finally, in section 5, we draw conclusions. ACS Paragon Plus Environment


Page 10 of 55 10

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

2. Methods To study the relaxation to the equilibrium of a supersaturated solution of methane in water we use the dynamical approach to non-equilibrium molecular dynamics (DNEMD) of Orlandini, Meloni and Ciccotti.46 In this extension of a previous method introduced by Ciccotti and co-workers,47,48 the initial state is determined by a macroscopic condition, being the concentration of well solvated solute molecules in an homogeneous methane/water solution in the present case. Let us denote (Γ) as a generic observable which is a function of the point in the phase space Γ = (, ). In a

non-equilibrium system the expectation value of this observable, () , is the

ensemble average over the probability density function (PDF) (Γ, ): () = (Γ) (Γ, t) dΓ (Γ, t)

satisfies

the

Liouville

(1) equation

(Γ, t)/ = − (Γ, t), ! (Γ, t)" =

−#$()(Γ, t), with {·,·} denoting the Poisson bracket and L(t) the Liouville operator. Generally speaking, the direct (exact or numerical) evaluation of f (Γ, t) is not possible apart for simple cases. However, one can resort to the Onsager-Kubo relation () = (Γ) % ∗ ()(Γ, t ' ) dΓ = %() (Γ) (Γ, t ' ) dΓ

(2)

where S (t) is the time-evolution operator of the phases space and S*(t), its adjoint, is the time-evolution operator of the PDF. f (Γ, t0) is the PDF at the initial time t0. Eq. (2) suggests a possible alternative method to compute O ( t ) : i) run trajectories starting from a set of phases space points Γ'( "()*,+ extracted from (Γ, ' ) ; ii) along these trajectories compute ,Γ(, Γ' )- , where Γ(, Γ' ) denotes the time-evolution of at Γ' time t; iii) average over the trajectories started from the various initial phase space (

points () = 1/. ∑+ ()* ( Γ (, Γ' )) .


Page 11 of 55


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

We assume (Γ, ' ) to be the PDF of a system at constant number of particles, N, temperature, T, and pressure, P, with the additional condition that the solution is at a prescribed value of the molar fraction χ of solvated methane molecules in solution: (Γ, t ' ) = (Γ| χ3() = χ∗ ) . We stress that 4 ∗ may differ from the (stable or metastable) equilibrium value at the present thermodynamic conditions; this is not a violation of any fundamental law because χ is only the most probable value of the methane molar fraction, not the only one that the system can take at the given thermodynamic condition. Summarising, to implement D-NEMD, we are left with the problem of extracting the set of initial conditions Γ'( "()*,+ from the conditional PDF (Γ| χ3() = χ∗ ) . This can be done, for example, by using the restrained MD method (RMD) of Maragliano and Vanden-Eijnden. 49 Assuming the system is governed by a suitable Hamiltonian

(Γ) , with standard

techniques, e.g., via a Nosé-Hoover chain thermostat,50 or incorporation of a Melchionna et al. barostat,51 one can sample the respective NVT or NPT ensembles. Maragliano and Vanden-Eijnden49 have shown that using the modified Hamiltonian 5 (Γ; 4

∗)

5

= (Γ) + 8 (4̂ () − 4 ∗ )8

(3)

in conjunction with the standard MD techniques mentioned above, and for sufficiently large values of : , one samples the conditional PDF (Γ| χ3() = χ∗ ). Here 4̂ (Γ) is the microscopic observable associated with the macroscopic property to be controlled, and : is a tunable parameter controlling how tightly the condition 4̂ (Γ) = 4 ∗ has to be met. The key point to appreciate in this above MD sampling vis-à-vis the relevant PDF is that the term :/2 (4̂ (Γ) − 4 ∗ )8 acts as an additional potential and, thus, the configurational part of the distribution is



Page 12 of 55 12

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

5 (;|4 ∗ ) =

F ()@E∗ )G H ?@AB()CD@A5/8(E 5 F ()@E∗ )G H I ?@AB()CD@A /8(E

(4)

Taking the limit for infinite k, one has the Gaussian J2K/L: MNOP−L :/2 (4̂ () −

4 ∗ )8 Q → S (4̂ () − 4 ∗ ), and thus lim5→W 5 (;|4 ∗ ) =

F ()@E∗ ) ?@AB()CX(E F ()@E∗ ) I ?@AB()CX(E

=

Y(Z,E∗ ) Y(E∗ )

= (;|4 ∗ ) . Summarising, it is

possible to sample the relevant conditional PDF by performing MD driven by the Hamiltonian

5 (Γ; 4

∗)

with a large value of k, which, in the following, is referenced

with the term ‘restrained MD’ (ReMD). To investigate hydrate nucleation at given concentration of methane, we need to sample the conditional PDF of a homogeneous liquid (a liquid containing neither methane bubbles nor clathrate nucleus) at the prescribed methane molar fraction. To

\ ]^ () , ` a () and b3cd () , defined do this, we use three collective variables, [. _ _ below. The first imposes that the liquid does not contain any methane bubble; the second prevents water and methane from forming nuclei of clathrate before the constraints are released, especially at high methane concentrations; the third controls the number of methane molecules within a prescribed volume of the computational sample, i.e., methane concentration. To impose the homogeneity of the liquid, i.e. the absence of methane bubbles in the solution, we control the arithmetic mean over CH4 molecules in solution of the methane-methane coordination number:

\ ]^ () = 1/ [. _ .

]^_

∑()*,+ij ∑gh(D1 − Θ,f;( − ;g f − ;]^_ -H _

(5)

where .]^_ is the number of methane molecules in the liquid domain of the sample, ;( and ;g are the positions of two methane molecules and Θ(⋅) is the Heaviside step function, which is equal to 0 when its argument is negative, and 1 otherwise. ;]^_ =


Page 13 of 55


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

5.5 Å is the radius of the first methane-methane coordination shell in methane bubbles ( () \ ]^ in water. The sum [. = ∑gh(D1 − Θ,f;( − ;g f − ;]^_ -H counts the number of _

methane molecules, j, coordinated to the molecule i. For ;]^_ = 5.5 Å the methanemethane coordination number is negligible for solvated CH4 but large in the case of a

\ ]^ () to be close to zero one enforces a methane bubble (Fig. 1). Thus, restraining [. _ homogeneous water-methane solution. To prevent any impulsive forces on the atoms arising from the biasing potential when they cross the threshold ;]^_ , the Heaviside step function is replaced by a smooth approximation, the Fermi function: *

Θ,f;( − ;g f − ;]^_ - = *nD@o,fZ @Z p

(6)

q f@Zij_ -H

where λ is a parameter controlling the smoothness of the approximation.52,53 In order to prevent crystallisation of water we control the ` a () observable,54 which is the mean of the cosine of three times the dihedral angle ξ(,g,5,s among all the fourterm chains of hydrogen-bonded water molecules in the simulation box: * ` a () = + ∑()*,+jG{ ∑g)*,yp ∑5)*,yq ∑s)*,yz cos,3ξ(g5s -

(7)

t

where .| is the number of four-terms water chains. In order to prevent impulsive forces when a water molecule enters or exits from the coordination shell, like in the

\ ]^ we replace Θ with a Fermi function: case of [. _

*

*nD@o,fZp @Zq f@;

2

. The ξ(,g,5,s

-H

dihedral angles in liquid water take on all possible values, with a 25 % difference between the most (60°, 180°, and 300°) and least (120° and 240°) probable angles.55 Because of this almost uniform distribution of ξ(,g,5,s the value of F4 of liquid water is zero.55 In clathrates, where water molecules form (essentially) planar five- and sixterm rings, most of the dihedral angles take values close to 0°. In addition, four-terms chains can be formed by water molecules belonging to different faces of a cage, or



Page 14 of 55 14

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

different cages, whose dihedral angles take values different from 0° (the value depends on the types of faces and cages). Thus, in bulk methane clathrate we observe a value of F4 ≈ 0.8.56 In practice, for each water molecule one searches for water molecules falling within its first coordination cell, i.e., within an oxygen-oxygen distance of 3.5 Å, corresponding to the first minimum of the O-O pair correlation function. The procedure is iterated for all neighbors and for neighbors of neighbors (Fig. 2). Given the position of the oxygen atoms of the four water molecules forming the chain, one can compute the dihedral angle ξ(g5s between the planes identified by the first and last three molecules of the chain. Concerning the molar fraction, 4̂ cd_ () is defined as the number of moles of solute molecules over the total number of moles of solvent plus solute. At the low concentrations of methane/water solution, in which the number of moles of solute is negligible with respect to the number of moles of solvent, and at given water density, the methane molar fraction is proportional to the number of methane molecules in the reference volume, b3cd_ (). Thus, we use this simpler collective variable, b3]^_ , in our calculations:

b3]^_ () = ∑()*,+ij_ ΔB (;( )

(8)

where ΔB (∙) is the characteristic function of the reference volume, which is equal to 1 if ;( lies within the reference volume and zero otherwise. As discussed in the following, the reference volume (Fig. 3) consists of a slice of thickness Δz that spans the complete box in the x-y directions, parallel to the interface between the reservoir of methane and the solution, and a portion of the box between the lower and upper limit s and , respectively, in the direction. Thus, the characteristic function can be written with only reference to the component of the atomic position: ACS Paragon Plus Environment

Page 15 of 55


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

Δ(;( ) = Δ(( ) = Θ(( − s ),1 − Θ(( − )-

(9)

Like in the cases of the other collective variables, the Heaviside step function is replaced by a Fermi function.

\ ]^ () and ` a () to zero to Using the ReMD approach, we restrain the values of [. _ enforce the homogeneity of the solution, and b3cd_ () to a suitable value to set the methane mole fraction to the desired concentration.

3. Simulation Details In all of our simulations, we use the ‘mW’ force field of Molinero and coworkers,57,58 a Stillinger-Weber-like interaction potential59 in which methane and water molecules are represented as point particles. This model consists of a sum of pairwise and three-body contributions, the latter encouraging the formation of tetrahedrally-coordinated molecules to reproduce “hydrogen-bonded” structures. In particular the three-body term adds an energy cost to configurations with nontetrahedral angles. Methane-methane and methane-water interactions are described by the two-body term of a Stillinger-Weber type potential. Using this force-field, one obtains a significant reduction in computational costs with respect to other models for two main reasons: firstly, using the united-atoms representation of water and methane molecules significantly reduces the number of particles composing the system. Secondly, the united-atom model, with the elimination of fast vibrating O-H bonds, allows to use a much longer timestep. Thirdly, the absence of electrostatics allows to avoid the computationally expensive method of Ewald summation or equivalent techniques. The mW model has been critiqued vis-à-vis other more popular models for methane hydrates60 and ice/water crystallisation/melting.61 In general, the mW family is very good for prediction of thermodynamic properties, but has somewhat



Page 16 of 55 16

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

artificially fast kinetics and dynamic properties. Still, bearing particularly in mind the low computational cost, mW constitutes a rather excellent approach for obtaining qualitative (indeed, semi-quantitative) insights into hydrate nucleation and formation. The computational samples used in our simulations consist of an orthorhombic box elongated along the z-axis with three-periodic boundary conditions of dimension

$= ~$ ~50 and $ ~110 Å (Fig. 3, the symbol “ ~ ” is used because in constant pressure simulations the size of the simulation box fluctuates). The box contains two phases; the first is an aqueous solution of solvated methane at the given concentration and the second is a methane liquid phase. The interfaces between the two phases are flat and perpendicular relative to the z-axis. For all samples, the water-methane solution occupies a portion of the box long ~ 80 Å in the z direction, with flat watermethane interfaces positioned at z1 ~ −42 Å and z2 ~ 42 Å. Within the solution domain we define the sub cell of ~ 60 Å in the z direction (Fig. 3) over which we compute

\ ]^ (), ` a () and b3cd (), and impose the conditions for the collective variables, [. _ _ the sampling of the initial probability density function. We considered four samples at methane mole fraction 4cd_ , 0.038, 0.044, 0.052 and 0.058, labelled A, B, C and D, respectively. The methane and solution domains of the simulations boxes were randomly filled up with the methane and water molecules, with numbers specified in Table 1. For each of the initial configurations, we performed a NPT ReMD at 250 K and 500 atm with a 1 fs time step for 20 ns to equilibrate the system. We remark that the timestep adopted here is shorter than the typical value used with the mW, 5 fs, because in ReMD the restraints forces are stiffer than the physical ones. We used the modified form of the Hoover barostat proposed by Melchionna et al.51All the


Page 17 of 55


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 were performed using our in-house modified version of the DL POLY Classic package.62 Given the fluctuations of the simulation box in NPT MD, the definition of the borders of sub cell has to be specified at each time step. Assuming that the bulk modulus is the same in the different regions of Fig. 3, the position of the boundaries of the box , and , obey the law

̅() =

()

(')

̅(0)

(10)

and the same for . In Eq. 10 $(0) and $() is the length of the box along the direction at the time 0 and t, respectively.

( ()

)

The conditional PDF f Γ χ r = χ * at each methane molar fraction is sampled taking a number of configurations (see Table 1) along the corresponding ReMD trajectory at a distance of 10 ps from each other after the initial 10 ns of thermalisation. Following the D-NEMD algorithm described in section 2, from each of these configurations we started unrestrained NPT simulations at 250 K and 500 atm. In this second set of calculations, we found it possible to extend the time step to 2 fs, which is still shorter than the typical 5 fs value. We made this conservative choice to avoid problems with initial out-of-equilibrium configurations characterized by high atomic forces. All of the relaxation trajectories were evolved in time until a nucleation event was observed. Finally, as a reference we computed the equilibrium methane molar through (standard) MD simulation over 100 ns, obtaining a value 4cd_ ~ 0.0025.

4. Results and Discussion In a previous work,27 we computed the minimum free-energy path (MFEP) connecting the local and global minima, i.e. the most probable nucleation path, using ACS Paragon Plus Environment


Page 18 of 55 18

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 string method63 on a free energy landscape obtained from the well-tempered metadynamics.64 We also computed the free-energy nucleation barrier, 173 kBT, defined as the difference in free energy between the maximum along the nucleation path (saddle point of the free-energy landscape) and local minimum corresponding to the two-phase initial condition, consisting of the methane-water solution plus a drop of methane. However, the nucleation path and kinetics can be affected by some of the hypotheses underpinning the rare-event techniques used in Ref. 27, namely that the sampling is performed in quasi-static (unrealistic) conditions and that methanemethane and methane-water coordination numbers are sufficiently good observables to describe the whole process. Here, using D-NEMD, we want to go beyond the inherent limitations of the simulation techniques used in this previous work. As discussed in Section 2, we study nucleation starting from a high concentrated methane-water solution. The reason of this choice is that the MFEP identified in Ref. 27 consisted of three steps, summarized in the following: the first step represents an increase of concentration of methane in solution and simultaneous formation of a disordered cluster of water and solvated methane molecules, i.e., the amorphous ‘blob’21 of BH or a cluster of ‘face saturated incomplete cages’ of CAH {Ref. 2 referee} shown in Fig. 4a. Since BH and CAH are both compatible with our results, hereafter we will use the term ‘blob’ with reference to both the blob of BH and cluster of ‘face saturated incomplete cages’ of CAH. In the second step, the blob is sufficiently large to be stable, i.e. it can resist to rapid regression to the initial metastable solution, and an ordering of its methane and water molecules takes place bringing to the formation of methane-filled [512]-, [51262]-, or [51264]-water cages (Fig. 4b). These cages are not yet arranged in a configuration consistent with any of the clathrate crystal structure. In the third step, the crystalline nucleus grows and the water cages undergo an ordering process producing a nucleus with a structure II clathrate ACS Paragon Plus Environment

Page 19 of 55


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

hydrate (Fig. 4c).65,66 Along the process the methane present in the bubble is consumed to maintain the concentration of methane in solution approximately constant at the high level of the first step, while the methane in solution is absorbed by the blob and used to grow the clathrate crystallite. This three-step mechanism suggests that a local fluctuation of methane concentration is a necessary precondition for homogeneous nucleation. The initial state of the system, as identified from the configurations sampled during the ReMD simulations, consists of a metastable methane-water blob involving 8090% of the CH4 molecules of the solution domain at a molar fraction, 4 = 0.06 −

0.07 , higher than the nominal value of the solution, 4 = 0.036 − 0.058 . Solvated methane molecules belonging to the blob are at a distance shorter than 9 Å each other, i.e. within the typical maximum distance of CH4 molecules belonging to sI or sII clathrate structures at the relevant temperature. Along D-NEMD trajectories the system first fluctuates around the initial composition, then overcome the nucleation barrier and its composition suddenly increases (Fig. 5). Indeed, the growing nuclei do not reach stoichiometric clathrate composition because of the mismatch between the simulation box and the randomly oriented clathrate nucleus. This prevents the complete formation of a single clathrate crystallite, resulting in the growth of a large crystalline domain surrounded by an amorphous region at lower methane concentration. Nucleation kinetics is analysed through the first-passage time (FPT), Y , and its probability density (FPTD), O(Y ). FPT is defined as the first time at which the largest cluster of connected cages contains a given number of cages, b : Y (b) =

#b > 0| b| b". FPTD is the distribution of this (stochastic) time, i.e., the probability density that the threshold b is reached for the first time at time Y . The



Page 20 of 55 20

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

number of connected cages, n, is computed with a code developed by Jacobson et al.65 and based on the graph theory, where water molecules are the vertices of the graph and neighboring molecules, within a distance of 3.5 Å from each other, form its edges. Analysing the connectivity and ordering of the edges it is possible to identify complete clathrate cages.27,66 Wedekind et al.37 have shown that for activated processes Y (b) follows the law: ∗

Y (b) = ?1 + M;((b − b∗ ))C, b ≥ 0 8

(11)

where M;(N) = 2/ ' M @ with −∞ < N < ∞ , and b∗ is the inflection point of √K =

G

Y (b), corresponding to the transition state which is the size of the crystalline nucleus

expressed in number of cages. determines how suddenly the Y changes with b. ∗ , the stationary (asymptotic) nucleation time, is related to the (homogeneous, volume independent) nucleation rate by y = 1/( ∗ ), where is the volume of the solution domain. We remark that recently Jungblut and Dellago have shown that Eq. (11) overestimates the nucleation time because of the non-Markovian character of the dynamics caused by the projection to a poorly chosen reaction coordinate.67 In this article, the error on the nucleation time of a Lennard-Jones crystal is estimated to be of a factor two. However, here we are more interested in the dependence of ∗ on methane concentration than its absolute value, and, as we show below, ∗ changes of 1-2 orders of magnitude with χ in the range considered. Thus, we expect that this error does not affect significantly our qualitative and quantitative conclusions. The time needed to produce a (post-critical) crystallite of size b is well expressed by Eq. (11) only if the (nucleation) time ∗ is much larger than that needed to go from b∗ to b , with b large enough to be experimentally detectable. In other words, it is assumed that growth is much faster than nucleation. We notice that at the lowest


Page 21 of 55


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

concentration the MFTP curve reaches a plateau (Fig. 6), indicating that the growth is much faster than nucleation. At higher concentrations, however, no plateau is observed, indicating that nucleation and growth occur on similar timescales. In this case, Y (b) is not well described by Eq. (11) (inset of Fig. 6). Nevertheless, Y (b) still present two well distinguished regimes, characterized by different slopes of the curve. This indicates that at all concentrations the process can still be described by the nucleation and growth steps. Chkonia et al.36 have shown that in this case the fitting of the early stage of the curve gives a reasonable estimate of the nucleation time. Here, however, we take a different approach that consists in adding to Eq. (11) a term expressing the contribution coming from the growth step to Y (b).38 Thanks to this approach we will be able to estimate the growth rate, which being relatively slow, with a rate that depends on the solution composition, represents an important element of the crystallization kinetics. This will help us understanding the possible effect of concentration on the overall crystallization process. Assuming that the grow rate is constant, b/ = = 1/ , we have ∗

Y (b) = 8 ?1 + M;((b − b∗ ))C + (b − b∗ )Θ(b − b∗ )

(12)

where the Heaviside step function, Θ(b − b∗ ), acts to turn on the growth term only after nucleation. Figs. 6 shows that Eq. 12 fits the present results very well, confirming the hypotheses of the assumed model. From the fitting of Y (b), one obtains the nucleation time at each given methane molar fraction, from which one gets the dependence of ∗ on 4cd_ . In particular, we found values for ∗ equal to 35.2, 4.5, 0.9 and 0.4 ns for 4cd_ = 0.038, 0.044, 0.052 and 0.058, respectively, with associated rates, y , 0.16, 1.3, 6.3 and 12.6 × 1033 m−3s−1. The increase of the rate with methane molar fraction suggests a corresponding decrease of the barrier with the concentration of methane in solution.



Page 22 of 55 22

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 output of the fitting, in particular b∗ , ∗ , , allows us to garner some insight into the crystallisation mechanism. The size of the critical nucleus b∗ is rather small at all the methane concentrations considered in this work. We measured a size between 4 and 2 for a saturation ratio % =

4cd_ 4 between ~15 and 23, at variance with cd_

b∗ ~300 for a saturation ratio of %~6 (Knott et al.68) and b∗ ~100 for a saturation ratio of %~8 (Lauricella et al.27) for the same force field model. This indicates a very rapid decrease of the critical size with saturation ratio. The combination of present results, focusing on nucleation in a highly concentrated solution, and literature data,24 can be used to suggest a clathrate crystallization mechanism. Walsh et al., in their investigation of nucleation of clathrates based on extensive standard MD, observed an increase of the methane mole fraction to the levels investigated in this work before observing nucleation events. This confirms a multi-step mechanism, consistent with the CAH and BH: through increase of the local methane concentration, the system is brought in conditions in which the critical nucleus consists of few clathrate cages, which can easily be formed by spontaneous fluctuations. In addition, the fit of Y (b) tells us that higher levels of saturation ratio increase the grow rate, which passes from ~17 to 25, 43, 63 cages/ns for 4cd_ going from 0.038 to 0.044, 0.052, 0.058. This brings to a quick formation of large clathrate crystallite, which can better resist to dissolution in presence of a reduction of the local concentration of methane. Complementary to the previous analysis, we can show the dependence of the nucleation barrier on saturation ratio by using the cumulative distribution function (CDF) of the nucleation time, f FPT (tf), i.e. the probability that a nucleation event


Page 23 of 55


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 > b∗ ) has taken place within time tf (Fig. 7). Higher saturation ratios affect the CDF making it steeper, i.e., the corresponding distribution of nucleation time is sharper. This indicates that at higher saturation ratios nucleation becomes more deterministic, which is typically due to a decrease of the nucleation barrier. This result is consistent with the reduction of the average nucleation time ∗ with % discussed above. The nucleation time can be obtained with different techniques. In particular,69 the rate y is computed from the survival probability (), i.e., the probability that no nucleation event is occurred up to time t, given by

() = exp?−y C = exp ?− / ∗ C

(13)

where is the volume of the methane-water solution. Eq. (13) is the solution of the first order master equation for the formation of a supercritical methane clathrate nucleus under the hypothesis that thermodynamic conditions do not change during nucleation, which implies that y is constant. The advantage of this approach is that, in experiments, it is simpler to measure the survival probability. In our simulations,

( ) = 1 − (). Thus, from the linear fitting of log ,Y - vs Y (Eq. (13)), one obtains ∗ = 1/( ) (Fig. 8). The nucleation times so obtained were 29.3, 4.3, 0.5 y and 0.2 ns (corresponding to y = 0.2 , 1.4, 11.5 and 33.5 × 1033 m−3s−1) for 4cd_ = 0.038, 0.044, 0.052 and 0.058, respectively, consistent with the corresponding results obtained from the mean first-passage time. Let us now discuss a quantitative relation between the nucleation rate (or time) and saturation ratio . Nucleation theory suggests that the frequency of nucleation events is ¢

y = ¡ exp P− A£ Q

(14)

¤



Page 24 of 55 24

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 ¡ is a weak function of supersaturation and depends on nucleus shape, interfacial energy, and temperature.70,71 In the exponent, the factors of L¥¦ , with

L = 1/:§ ¨ , are often approximated by the square of the natural logarithm of the supersaturation, which is here replaced with the saturation ration % , i.e. L¥¦ =

ln8 (%). Here we assume that ¡ and do not do not depend on %. Although the values of ¡ and predicted by CNT (Eq. 14) are inaccurate, and S is the saturation ratio rather than supersaturation, the general form of this equation can be used to fit the data quite well.70,72-75 We report in Fig. 9 the nucleation rate log y versus 1/ln8 (%) for the four values of saturation ratio and the associated linear-regression line; this gives a = 1.69×1040 m−3s−1 and b = -135.89 for the mean first passage time and a = 1.95×1041 m−3s−1 and b = -153.71 for the survival probability, respectively. Extrapolating from Eq. 14 the nucleation time at lower values of local saturation ratio on a nanoscopic volume (53 nm3, Fig. 10) one obtains that at 4cd_ = 0.03 (corresponding to % =12) the ∗ is in the microsecond timescale and at 4cd_ = 0.01 (corresponding to % =4) the value ∗ is of the order of months. This confirms the hypothesis that forming a high concentrated solution (by spontaneous local concentration fluctuations, water freezing,40 etc) is a precondition of homogeneous nucleation and that this fluctuation on a nanoscopic portion of solution is sufficient to induce nucleation. In particular, nucleation occurs on the relevant experimental timescale, ∗ ~1 − 60ª, for 4cd_ = 0.0175-0.015 (corresponding to %=7-6).

5. Conclusions In this work, we have investigated the mechanism of (homogeneous) nucleation of methane hydrate as a function of methane concentration in water using the D-NEMD method. Nucleation is a very complex phenomenon and the role of methane


Page 25 of 55


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

concentration in controlling the process has still to be clarified. In particular, here we computed the nucleation rate at four different methane molar fraction conditions, and obtained the theoretical curve rate vs saturation ratio. This allows for estimating the nucleation rate at thus-far unexplored methane concentrations. We observed that slight variations in methane saturation ratio leads to huge effects on the nucleation rates, confirming that the local increase of methane mole fraction is the rate-determining step in methane-hydrate nucleation. An important remark is that the time scale necessary to observe nucleation for solvated-methane concentration of 0.0125, is of the order of months, suggesting that, most probably, the formation of methane hydrates in nature occurs either via heterogeneous nucleation or, whenever an homogeneous nucleation takes place, the homogenous process requires a relevant fluctuation in solvated methane concentration, such as that induced by the freezingassisted mechanism. This calls for the investigation in the future of heterogeneous nucleation of methane clathrates at suitable surfaces, where, again, advanced MD techniques will be necessary.

Acknowledgements The authors thank Valeria Molinero for providing her code for cages identification and interesting discussions. GC, SM and ML acknowledge financial support from SFI Grant No. 08–IN.1–I1869. GC and SM also acknowledge financial support from the Istituto Italiano di Tecnologia under the SEED project grant No. 259 SIMBEDD – Advanced Computational Methods for Biophysics, Drug Design and Energy Research. S.M. acknowledges financial support from the MIUR-FIRB Grant No. RBFR10ZUUK. The authors thank the Irish Centre for High-End Computing (ICHEC), and the Science Foundation Ireland (SFI) for the provision of computational resources.



Page 26 of 55 26

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

REFERENCES 1

Makogon, I. F. Hydrates of Hydrocarbons. Pennwell Books: Houston, U.S., 1997.

2

Sloan Jr, E. D.; Koh, C. Clathrate Hydrates of Natural Gases. CRC press: Boca Raton, U.S., 2007.

3

MacDonald, G. J. The Future of Methane as an Energy Resource. Annu. Rev. Energy 1990, 15, 53-83.

4

Kvenvolden, K. A. Methane Hydrate - a Major Reservoir of Carbon in the Shallow Geosphere? Chem. Geol. 1998, 71, 41-51.

5

Max, M. D. Natural Gas Hydrate in Oceanic and Permafrost Environments. Vol. 5. Springer Science & Business Media: Berlin, Germany, 2003.

6

Max, M. D.; Johnson A. H.; Dillon, W. P. Economic Geology of Natural Gas Hydrate. Vol. 9. Springer Science & Business Media: Berlin, Germany, 2005.

7

Max, M. D.; Johnson, A. H. Exploration and Production of Oceanic Natural Gas Hydrate. Springer International Publishing: Berlin, Germany, 2016.

8

Spagnoli, G.; Finkenzeller, S.; Freudenthal, T.; Hoekstra, T.; Woollard, M.; Storteboom, O.; Weixler, L. First Deployment of the Underwater Drill Rig MeBo200 in the North Sea and its Applications for the Geotechnical Exploration. SPE Offshore Europe Conference and Exhibition. Society of Petroleum Engineers: London, U.K., 2015.

9

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. AIChE J. 1983, 29, 651-657.

10

Feng, Y.; Xu, Y; Yu, Y.; Xie, Z.; Lin, X. Mechanisms of Biochar Decreasing Methane Emission from Chinese Paddy Soils. Soil Biol. Biochem. 2012, 46, 80-88.

11

Hammerschmidt, E. G. Formation of Gas Hydrates in Natural Gas Transmission Lines. Ind. Eng. Chem. 1934, 26, 851-855.

12

Sloan, E. D. Fundamental Principles and Applications of Natural Gas Hydrates. Nature 2003, 426, 353-363.


Page 27 of 55


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

13

Chatti, I.; Delahaye, A.; Fournaison, L.; Petitet, J. P. Benefits and Drawbacks of Clathrate Hydrates: a Review of their Areas of Interest. Energy Convers. Manage

2005, 46, 1333-1343. 14

Sloan, E. D.; Fleyfel. F. A Molecular Mechanism for Gas Hydrate Nucleation from Ice. AIChE J. 1991, 37, 1281-1292.

15

Guo, G. J.; Zhang, Y. G.; Zhao, Y. J.; Refson, K.; Shan, G. H. Lifetimes of Cagelike Water Clusters Immersed in Bulk Liquid Water: A Molecular Dynamics Study on Gas Hydrate Nucleation Mechanisms. J. Chem. Phys. 2004, 121, 15421547.

16

Guo, G. J.; Zhang, Y. G.; Li, M.; Wu, C. H. Can the Dodecahedral Water Cluster Naturally Form in Methane Aqueous Solutions? A Molecular Dynamics Study on the Hydrate Nucleation Mechanisms. J. Chem. Phys. 2008, 128, No. 194504.

17

Radhakrishnan, R.; Trout, B. L. A New Approach for Studying Nucleation Phenomena Using Molecular Simulations: Application to CO2 Hydrate Clathrates. J. Chem. Phys. 2002, 117, 1786-1796.

18

Guo, G. J.; Li, M.; Zhang, Y. G.; Wu, C. H. Why Can Water Cages Adsorb Aqueous Methane? A Potential of Mean Force Calculation on Hydrate Nucleation Mechanisms. Phys. Chem. Chem. Phys. 2009, 11, 10427-10437.

19

Guo, G. J.; Zhang, Y. G.; Liu, H. Effect of Methane Adsorption on the Lifetime of a Dodecahedral Water Cluster Immersed in Liquid Water: a Molecular Dynamics Study on the Hydrate Nucleation Mechanisms. J. Phys. Chem. C 2007, 111, 25952606.

20

Guo, G. J.; Zhang, Y. G.; Liu, C. J.; Li, K. H. Using the Face-Saturated Incomplete Cage Analysis to Quantify the Cage Compositions and Cage Linking Structures of Amorphous Phase Hydrates. Phys. Chem. Chem. Phys. 2011, 13, 12048-12057.

21

Jacobson, L. C.; Hujo, W.; Molinero, V. Amorphous Precursors in the Nucleation of Clathrate Hydrates, J. Am. Chem. Soc. 2010, 132, 11806-11811.

22

English, N. J.; Lauricella, M.; Meloni, S. Massively Parallel Molecular Dynamics Simulation of Formation of Clathrate-Hydrate Precursors at Planar Water-Methane Interfaces: Insights into Heterogeneous Nucleation. J. Chem. Phys. 2014, 140, No. 204714.



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

23

Hawtin, R. W.; Quigley, D.; Rodger, P. M. Gas Hydrate Nucleation and Cage Formation at a Water/Methane Interface. Phys. Chem. Chem. Phys. 2008, 10, 4853-4864.

24

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

25

Jacobson, L. C.; Molinero, V. Can Amorphous Nuclei Grow Crystalline Clathrates? The Size and Crystallinity of Critical Clathrate Nuclei. J. Am. Chem. Soc. 2011, 133, 6458-6463.

26

Guo, G. J.; Rodger, P. M. Solubility of Aqueous Methane under Metastable Conditions: Implications for Gas Hydrate Nucleation. J. Phys. Chem. B 2013, 117, 6498-6504.

27

Lauricella, M.; Meloni, S.; English, N. J.; Peters, B.; Ciccotti, G. Methane Clathrate Hydrate Nucleation Mechanism by Advanced Molecular Simulations. J. Phys. Chem. C 2014, 118, 22847-22857.

28

Sarupria, S.; Debenedetti, P.G. Homogeneous Nucleation of Methane Hydrate in Microsecond Molecular Dynamics Simulations. J. Phys. Chem. Lett. 2012, 3, 2942-2947.

29

Małolepsza, E.; Keyes, T. Pathways through Equilibrated States with Coexisting Phases for Gas Hydrate Formation. J. Phys. Chem. B 2015, 119, 15857-15865.

30

Petrenko, V. F.; Whitworth, R. W. Physics of Ice. Oxford University Press: Oxford, U.K., 1999.

31

Bi, Y.; Li, T. Probing Methane Hydrate Nucleation through the Forward Flux Sampling Method. J. Phys. Chem. B 2014, 118, 13324-13332.

32

Allen, R. J.; Frenkel, D.; ten Wolde, P. R. Simulating Rare Events in Equilibrium or Nonequilibrium Stochastic Systems. J. Chem. Phys. 2006, 124, No. 024102.

33

Bi, Y.; Porras, A.; Li, T. Free Energy Landscape and Molecular Pathways of Gas Hydrate Nucleation. J. Chem. Phys. 2016, 145, No. 211909.

34

Moon, C.; Taylor, P. C.; Rodger, P. M. Molecular Dynamics Study of Gas Hydrate Formation. J. Am. Chem. Soc. 2003, 125, 4706-4707.


Page 29 of 55


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

35

Kashchiev, D.; Firoozabadi. A. Kinetics of the Initial Stage of Isothermal Gas Phase Formation. J. Chem. Phys. 1993, 98, 4690-4699.

36

Peters, B. Supersaturation Rates and Cchedules: Nucleation Kinetics from Isothermal Metastable Zone Widths. J. Cryst. Growth 2011, 317, 79-83.

37 38

Nývlt, J. Kinetics of Nucleation in Solutions. J. Cryst. Growth 1968, 3, 377-383. Kubota, N. A New Interpretation of Metastable Zone Widths Measured for Unseeded Solutions. J. Cryst. Growth 2008, 310, 629-634.

39

Sangwal, K. On the Effect of Impurities on the Metastable Zone Width of Phosphoric Acid. J. Cryst. Growth 2010, 312, 3316-3325.

40

Poon, G. G.; Peters, B. A Stochastic Model for Nucleation in the Boundary Layer During Solvent Freeze-Concentration. Cryst. Growth Des. 2013, 13, 4642-4647.

41

Kashchiev, D.; Borissova, A.; Hammond, R. B.; Roberts, K. J. Dependence of the Critical Undercooling for Crystallization on the Cooling Rate. J. Phys. Chem. B 2010, 114, 5441-5446.

42

Chkonia, G.; Wölk, J.; Strey, R.; Wedekind, J.; Reguera, D. Evaluating Nucleation Rates in Direct Simulations. J. Chem. Phys. 2009, 130, No. 064505.

43

Wedekind, J.; Strey, R.; Reguera, D. New Method to Analyze Simulations of Activated Processes. J. Chem. Phys. 2007, 126, No. 134103.

44

Yuhara, D.; Barnes, B. C.; Suh, D.; Knott, B. C.; Beckham, G. T.; Yasuoka, K.; Wubd, D. T.; Sum, A. K., Nucleation Rate Analysis of Methane Hydrate from Molecular Dynamics Simulations. Faraday Discuss. 2015, 179, 463-474.

45

Barnes, B. C.; Knott, B. C.; Beckham, G. T.; Wu, D. T.; Sum, A. K. Reaction Coordinate of Incipient Methane Clathrate Hydrate Nucleation. J. Phys. Chem. B 2014, 118, 13236-13243.

46

Orlandini, S.; Meloni, S.; Ciccotti, G. Hydrodynamics from Statistical Mechanics: Combined Dynamical-NEMD and Conditional Sampling to Relax an Interface between Two Immiscible Liquids. Phys. Chem. Chem. Phys. 2011, 13, 1317713181.



Page 30 of 55 30

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

47

Ciccotti, G.; Jacucci. G. Direct Computation of Dynamical Response by Molecular Dynamics: The Mobility of a Charged Lennard-Jones Particle. Phys. Rev. Lett.

1975, 35, 789-792. 48

Ciccotti, G.; Jacucci, G.; McDonald. I. R. “Thought-Experiments” by Molecular Dynamics. J. Stat. Phys. 1979, 21, 1-22.

49

Maragliano, L.; Vanden-Eijnden, E. A Temperature Accelerated Method for Sampling Free Energy and Determining Reaction Pathways in Rare Events Simulations. Chem. Phys. Lett. 2006, 426, 168-175.

50

Martyna, G. J.; Klein, M. L.; Tuckerman, M. Nosé–Hoover Chains: the Canonical Ensemble via Continuous Dynamics. J. Chem. Phys. 1992, 97, 2635-2643.

51

Melchionna, S.; Ciccotti, G.; Lee Holian, B. Hoover NPT Dynamics for Systems Varying in Shape and Size. Mol. Phys. 1993, 78, 533-544.

52

Orlandini, S.; Meloni, S.; Colombo, L. Order-Disorder Phase Change in Embedded Si Nanoparticles. Phys. Rev. B 2011, 83, No. 235303.

53

Orlandini, S.; Meloni, S.; Ciccotti, G. Combining Rare Events Techniques: Phase Change in Si Nanoparticles. J. Stat. Phys. 2011, 145, 812-830.

54

Rodger, P. M.; Forester, T. R.; Smith, W. Simulations of the Methane Hydrate/Methane Gas Interface near Hydrate Forming Conditions. Fluid Phase Equilib. 1996, 116, 326-332.

55

Mezei, M.; Speedy, R. J. Simulation Studies of the Dihedral Angle in Water. J. Phys. Chem. 1984, 88, 3180-3182.

56

Rodger, P. Methane Hydrate: Melting and Memory. Ann. N. Y. Acad. Sci. 2000, 912, 474-482.

57

Molinero, V.; Moore, E. B. Water Modeled as an Intermediate Element between Carbon and Silicon. J. Phys. Chem. B 2008, 113, 4008-4016.

58

Jacobson, L. C.; Molinero, V. A Methane-Water Model for Coarse-Grained Simulations of Solutions and Clathrate Hydrates. J. Phys. Chem. B 2008, 114, 7302-7311.

59

Stillinger, F. H.; Weber, T. A. Computer Simulation of Local Order in Condensed Phases of Silicon. Phys. Rev. B 1985, 31, 5262-5271.


Page 31 of 55


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

60

Waldron, C. J.; Lauricella, M.; English, N. J. Structural and Dynamical Properties of Methane Clathrate Hydrates from Molecular Dynamics: Comparison of Atomistic and More Coarse-Grained Potential Models. Fluid Phase Equilib. 2016, 413, 235-241.

61

English, N. J. Massively Parallel Molecular-Dynamics Simulation of Ice Crystallisation and Melting: The Roles of System Size, Ensemble, and Electrostatics. J. Chem. Phys. 2014, 141, No. 234501.

62

Smith, W.; Forester, T. R. DL_POLY_2. 0: A General-Purpose Parallel Molecular Dynamics Simulation Package. J. Mol. Graphics 1996, 14, 136-141.

63

Maragliano, L.; Fischer, A.; Vanden-Eijnden, E.; Ciccotti, G. String Method in Collective Variables: Minimum Free Energy Paths and Isocommittor Surfaces. J. Chem. Phys. 2006, 125, No. 024106.

64

Barducci, A.; Bussi, G.; Parrinello, M. Well-Tempered Metadynamics: a Smoothly Converging and Tunable Free-Energy Method. Phys. Rev. Lett. 2008, 100, No. 020603.

65

Jacobson, L. C.; Hujo, W.; Molinero, V. Thermodynamic Stability and Growth of Guest-Free Clathrate Hydrates: a Low-Density Crystal Phase of Water. J. Phys. Chem. B 2009, 113, 10298-10307.

66

Lauricella, M.; Meloni, S.; Liang, S.; English, N. J.; Kusalik, P. G.; Ciccotti, G. Clathrate Structure-Type Recognition: Application to Hydrate Nucleation and Crystallisation. J. Chem. Phys. 2015, 142, No. 244503.

67

Jungblut, S.; Dellago, C. Caveats of Mean First-Passage Time Methods Applied to the Crystallization Transition: Effects of non-Markovianity. J. Chem. Phys. 2015, 143, No. 064103.

68

Knott, B. C.; Molinero, V.; Doherty, M. F.; Peters, B. Homogeneous Nucleation of Methane Hydrates: Unrealistic under Realistic Conditions. J. Am. Chem. Soc.

2012, 134, 19544-19547. 69

Kashchiev, D.; Verdoes, D.; Van Rosmalen, G. M. Induction Time and Metastability Limit in New Phase Formation. J. Cryst. Growth 1991, 110, 373380.

70

Kashchiev, D. Nucleation. Butterworth-Heinemann: Oxford, U.K., 2000. ACS Paragon Plus Environment


Page 32 of 55 32

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

71

Kalikmanov, V. Nucleation Theory, Springer: Dordrecht, Netherlands, 2013.

72

Davey, R. J.; Schroeder, S. L.; ter Horst, J. H. Nucleation of Organic Crystals-a Molecular Perspective. Angew. Chem., Int. Ed. 2013, 52, 2166-2179.

73

Kelton, K. F.; Greer, A. L.; Thompson, C. V. Transient Nucleation in Condensed Systems. J. Chem. Phys. 1983, 79, 6261-6276.

74

Kashchiev, D.; Van Rosmalen, G. M. Review: Nucleation in Solutions Revisited. Cryst. Res. Technol. 2003, 38, 555-574.

75

Kashchiev, D. Forms and Applications of the Nucleation Theorem. J. Chem. Phys.

2006, 125, No. 014502.

Table 1: Composition of systems simulated Label

No. of Trajectories

Methane molar fraction χCH4 in solution

No. of water molecules

Total No. of methane molecules

A B C D

64 64 64 40

0.038 0.044 0.052 0.058

7800 6847 6847 6847

1280 1153 1153 1153


No. of restrained methane molecules 240 240 276 312

Page 33 of 55


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 Captions Figure 1: Average coordination number CNCH4−CH4 in a methane bubble (dashed blue) and in an homogeneous solution at equilibrium concentration (red). Figure 2: Schematic diagram showing how the water hydrogen-bonded network is defined. One four-terms hydrogen bonded chain of mW water molecules is reported. The coordination shell of the first three molecules of the chain is also shown together with the vectors ri , j , r j , k and rk ,l joining pairs of oxygen atoms. Figure 3: a) Snapshot of the computational sample. Cyan and red spheres represent CH4 and H2O water molecules, respectively. The system consists of three domains: i) at the borders of the simulation box there is a (planar) drop of methane at the interface with ii) a methane-water solution, iii) the central sub cell, , highlighted in yellow, is the region of methanewater solution on which we enforce the initial methane concentration,

4]^_ = 0.038, 0.044, 0.052, 0.058 . b) Methane concentration profile

illustrating that in the central sub cell the concentration is uniform.

4]^_ in the region between the methane drop and the sub cell presents

a decreasing-increasing profile resulting from the tendency of the solution to reach the equilibrium value 4]^_ ~0.0025 and the formation of smooth interfaces with the neighboring domains at higher methane concentration. Figure 4: Snapshots at the system configuration during nucleation at 250 K and 500 atm starting from a two-phase system containing a drop of methane.65 Cyan and blue spheres represent CH4 in the bubble and wellhydrated methane molecules in solution (more than 16 H2O in the first coordination shell of the methane molecule). Water molecules belonging to the blob, sI and sII are colored in red, green and orange, respectively. Water molecules bridging sI and sII domains are colored in yellow. The rest of water molecules are colored in pink. Water molecules forming complete cages of types [512], [51262], or [51264] are represented as connected. Cages of this type are present also in the blob structure



Page 34 of 55 34

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

however they are not part of a sI or sII nucleus with a suitable combination of cages of different types. Nucleus structure recognition is performed by means of the algorithm discussed in Ref. 66. Figure 5: methane molar fraction in the blob/crystalline nucleus, 4¢s¬¢ , as a function of time along single D-NEMD trajectories for solutions at the four nominal composition 4]^_ = 0.038, 0.044, 0.052 and 0.058. Figure 6: Plot of the mean first passage time Y vs n at 4]^_ = 0.038 (a), 0.044 (b),

0.052 (c), 0.058 (d). Red symbols with error bars denote simulation results and the blue dashed curve the fitting to Eq 11. In the inset of

panel (b) it is reported the fitting of Y vs n with Eq. 10, i.e. without including the growth term, to illustrate the inadequacy of this model. Figure 7: Cumulative distribution function FPTD f FPT (τf), i.e., the probability that nucleation has generated a hydrate-like cluster with n > 12 Figure 8: (natural) log(PS) vs Y at 4]^_ = 0.038 (a), 0.044 (b), 0.052 (c), 0.058 (d). The symbols with error bars denotes results obtained by simulations and the red dashed line is the linear-regression line of Eq. 12. Note that PS was computed as the probability no nucleus larger than 12 cages is present in the sample. Figure 9: log (y ) vs 1/ log2 (S) at the four methane concentrations together with the linear-regression fit line. We report both the results obtained with the mean first passage time that with the survival probability. Figure 10: Nucleation time ∗ as function of solvated methane mole fraction 4]^_ for a box of volume 125x103 Å3. Together with the computational results (symbols) we reports the predictions obtained by using eqn. 13 (continuous lines), with parameters obtained from the fitting of mean first passage time (blue) and survival probability (red) simulation data.


Page 35 of 55


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

Fig. 1



Page 36 of 55 36

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

Fig. 2


Page 37 of 55


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

Fig. 3



Page 38 of 55 38

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

(a)

(b)

(c)

(d)

Fig. 4


Page 39 of 55


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

Fig. 5



Page 40 of 55 40

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

Fig. 6


Page 41 of 55


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

Fig. 7



Page 42 of 55 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

Fig. 8


Page 43 of 55


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

Fig. 9



Page 44 of 55 44

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

Fig. 10


Page 45 of 55

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



The Journal of Physical Chemistry Page 46 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15


Page 47 of 55

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




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


Page 48 of 55

Page 49 of 55

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




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 50 of 55

Page 51 of 55

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




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 52 of 55

Page 53 of 55

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




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


Page 54 of 55

Page 55 of 55

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



Mechanisms and Nucleation Rate of Methane Hydrate by Dynamical

Recommend Documents