Dongliang Jin

Benoit Coasne

Reduced phase stability and faster formation/dissociation kinetics in confined methane hydrate

The mechanisms involved in the formation/dissociation of methane hydrate confined at the nanometer scale are unraveled using advanced molecular modeling techniques combined with a mesoscale thermodynamic approach. By means of atom-scale simulations probing coexistence upon confinement and free energy calculations, phase stability of confined methane hydrate is shown to be restricted to a narrower temperature and pressure domain than its bulk counterpart. The melting point depression at a given pressure, which is consistent with available experimental data, is shown to be quantitatively described using the Gibbs-Thomson formalism if used with accurate estimates for the pore/liquid and pore/hydrate interfacial tensions. The metastability barrier upon hydrate formation and dissociation is found to decrease upon confinement, therefore providing a molecular scale picture for the faster kinetics observed in experiments on confined gas hydrates. By considering different formation mechanisms -bulk homogeneous nucleation, external surface nucleation, and confined nucleation within the porosity -we identify a crossover in the nucleation process; the critical nucleus formed in the pore corresponds either to a hemispherical cap or a bridge nucleus depending on temperature, contact angle, and pore size. Using the classical nucleation theory, for both mechanisms, the typical induction time is shown to scale with the pore volume to surface ratio and, hence, the pore size. These findings for the critical nucleus and nucleation rate associated to such complex transitions provide a mean to rationalize and predict methane hydrate formation in any porous media from simple thermodynamic data.

Gas hydrate | Confinement and porous media | Formation/dissociation kinetics | Thermodynamics and molecular modeling | M ethane hydrate is a non-stoichiometric crystalline phase in which water molecules form hydrogen-bonded cages that entrap methane molecules [START_REF] Sloan | Fundamental Principles and Applications of Natural Gas Hydrates[END_REF][START_REF] Koh | State of the art: Natural gas hydrates as a natural resource[END_REF]. Under typical terrestrial and marine conditions, methane hydrate forms according to a structure known as sI where 46 water molecules form two small pentagonal dodecahedral cages and six tetracaidecahedral cages so that 8 methane molecules can be encapsulated at most [START_REF] Sloan | Clathrate Hydrates of Natural Gases[END_REF]. While the large methane hydrate resources available on Earth are still regarded as an important fossil energy source [START_REF] Macdonald | The Future of Methane as an Energy Resource[END_REF], their abundant presence in locations such as seafloor or permafrost is also a threat to the environment as methane is one of the worst greenhouse gases [START_REF] Serov | Postglacial response of Arctic Ocean gas hydrates to climatic amelioration. 810[END_REF][START_REF] Phrampus | Recent Changes to the Gulf Stream Causing Widespread Gas 812 Hydrate Destabilization[END_REF]. For instance, in the context of increasing concerns about climate change, even weak perturbations such as those induced by a small temperature raise could trigger the dissociation of methane hydrate and release of large amounts of methane into the atmosphere [START_REF] Wood | Decreased Stability of Methane Hydrates in Marine Sediments Owing to 814 Phase-boundary Roughness[END_REF].

Gas hydrates including methane hydrate are also thought to be a role player in the geochemistry of planets, comets, etc. as typical temperature T and pressure P in many of these systems should promote their formation whenever water and 20 gases are present [START_REF] Formisano | Detection of Methane in the Atmosphere of Mars[END_REF]. Finally, methane hydrate is also of 21 particular relevance to energy and environmental science with 22 applications for energy storage and carbon capture [START_REF] Lee | Tuning clathrate hydrates for hydrogen storage[END_REF][START_REF] Patt | A grand canonical Monte Carlo study of the N 2 , CO, and mixed N 2 -CO clathrate 819 hydrates[END_REF][START_REF] Petuya | Selective trapping of CO 2 gas and cage occupancy in CO 2 -N 2 and CO 2 -CO 821 mixed gas hydrates[END_REF][START_REF] English | Mechanisms for thermal conduction in methane hydrate[END_REF][START_REF] Pohlman | Enhanced CO 2 uptake at a shallow Arctic Ocean seep field over-825 whelms the positive warming potential of emitted methane[END_REF].

23

The phase diagram of bulk methane hydrate has been [START_REF] Seshadri | Measurements of Equilibrium Pressures and Tempera-850 tures for Propane Hydrate in Silica Gels with Different Pore-Size Distributions[END_REF] the subject of intense research to identify P/T coexistence 25 conditions as well as to determine formation/dissociation mech-26 anisms [START_REF] Hall | Evidence from mixted hydrate nucleation for a funnel 828 model of crystallization[END_REF][START_REF] Duan | The influence of temperature, pressure, salinity and capillary force on the 830 formation of methane hydrate[END_REF]. In contrast, the case of methane hydrate in 27 porous media remains unclear by many aspects with important 28 questions left regarding the role of confinement and surface 29 forces [START_REF] Casco | Methane hydrate formation in confined nanospace can surpass nature[END_REF][START_REF] Wang | Direction Controlled Methane Hydrate Growth[END_REF][START_REF] Ghosh | Clathrate hydrates in interstellar environment[END_REF]. Yet, in nature, abundant methane hydrate re-30 sources are found in the porosity of rocks and minerals such as 31 in marine sediments, silica sands, permafrost, etc. so that un-32 derstanding confinement or surface-induced shifts in the phase 33 stability, formation/dissociation kinetics, and composition of 34 gas hydrates is of utmost importance (Fig. 1). Many exper-35 imental observations have shown that confinement in pores 36 shifts the liquid-hydrate-vapor (L-H-V) phase boundaries to 37 higher P and/or lower T as a result of capillary forces (19-38 23). For both organic (e.g. carbonaceous) and mineral (e.g. 39 siliceous) porous environments, such experiments have shown behavior is found to be reminiscent even at the nanometer 45

Significance Statement

Large amounts of methane hydrate -which consists of water forming molecular ice cages around methane molecules -are trapped in the porosity of rocks and seafloors. The impact of this so-called "burning ice" on our environment remains unclear (methane is a potent greenhouse gas), especially the parameters driving its phase stability and kinetics. Here, using molecular modeling, we identify the mechanisms involved in the formation/dissociation of confined methane hydrate. The phase stability of confined methane hydrate is shown to be restricted to a narrower temperature/pressure range with faster underlying formation/dissociation kinetics. Such facilitated freezing, which helps rationalize experimental observations, is described using a mesoscale thermodynamic approach of the nucleation mechanism. In fact, while the question of a contact angle θ between solid and liquid phases at a pore surface remains to be addressed, many experimental and theoretical works assume that there is a wetting parameter ω = γ LH cos θ that controls the shift in the melting point ∆T f ∼ ω/Dp (with the usual assumption cos θ ∼ -1). Beyond confinement-induced shifts in hydrate stability, the formation and dissociation kinetics of these complex compounds remain to be investigated in detail. There are abundant literature data indicating that confinement leads to faster kinetics compared to their bulk counterpart [START_REF] Casco | Methane hydrate formation in confined nanospace can surpass nature[END_REF][START_REF] Borchardt | Methane Hydrate in Confined Spaces: An 840 Alternative Storage System[END_REF][START_REF] Linga | Enhanced rate of gas hydrate formation in a fixed bed column filled with sand compared to a stirred vessel[END_REF] but the underlying microscopic mechanisms remain to be unraveled. Theoretically, much work has been done to elucidate the molecular routes that lead to the nucleation of bulk hydrates [START_REF] Walsh | Microsecond simulations of spontaneous methane hydrate nucleation and growth[END_REF][START_REF] Yuhara | Nucleation rate analysis of methane hydrate from molecular dynamics simulations[END_REF][START_REF] Knott | Homogeneous Nucleation of Methane Hydrates: Unrealistic under Realistic Conditions[END_REF][START_REF] Arjun | Unbiased atomistic insight in the competing nucleation mechanisms of methane hydrates[END_REF][START_REF] Liang | Nucleation and dissociation of methane clathrate embryo at the gas-water interface[END_REF] but the case of confined hydrates has not been considered. Several authors have used molecular modeling to investigate the growth of an already formed methane or carbon dioxide hydrate in the vicinity of a surface [START_REF] Kang | Kinetics of Methane and Carbon Dioxide Hydrate Formation in Silica Gel Pores[END_REF][START_REF] Bai | Microsecond Molecular Dynamics Simulations of the Kinetic Pathways of Gas Hydrate Formation from Solid Surfaces[END_REF][START_REF] Yan | Molecular Dynamics Simulation of the Crystal Nucleation and Growth Behavior of Methane Hydrate in the Presence of the Surface and Nanopores of Porous Sediment[END_REF][START_REF] Liang | Crystal growth simulations of methane hydrates in the presence of silica surfaces[END_REF][START_REF] Bagherzadeh | Influence of Hydrated Silica Surfaces on Interfacial Water in the Presence of Clathrate Hydrate Forming Gases[END_REF] but microscopic nucleation remains largely unexplored.

D R A F T

Here, we present a theoretical study of confined methane hydrate to elucidate the mechanisms responsible for their reduced phase stability and faster formation/dissociation kinetics. By relying on statistical mechanics, such molecular simulations -using either the direct coexistence method (DCM) or free energy calculations -capture phase coexistence between confined liquid water and methane hydrate. In particular, our approach does not assume any thermodynamic modeling framework or kinetic pathway. We first extend the direct coexistence method to the Grand Canonical ensemble (open ensemble)

to account for three phase coexistence (L-H-V) involved in gas hydrate formation/dissociation. Using such atom-scale simulations, we show that the melting point depression for confined hydrate is consistent with the Gibbs-Thomson equation if the pore/liquid and pore/hydrate interfacial tensions are used (here, values were calculated using independent molecular simulations). This macroscopic expression is shown to be a simplified thermodynamic approach that neglects the effect of pressure/temperature and methane chemical potential. In agreement with experimental data, the formation and dissociation of confined methane hydrate as described in our S1) are used to mimic stronger or weaker surface interactions.

As shown in Fig. S4, T f remains constant as is varied. This result can be explained using the Gibbs-Thomson equation [Eq. ( 3)] which will be discussed in detail below. The shift

∆T f /T f,0 is proportional to the interfacial tension difference ∆γ = γ WL -γ WH .
At constant T and P , a Taylor expansion for ∆γ leads to ∆γ( ) ∼ ∆γ( )+∂∆γ( )/∂ ×( -). Surface interactions are found to amount for < 5% of the energy with negligible impact on ∆γ. As a result, ∆γ( ) ∼ ∆γ( ) so that wettability considered here leads to the same melting shift. 199

Gibbs

Let us consider two phases Φ, methane hydrate (Φ = H) 200 and liquid water with solubilized methane (Φ = L), in a slit 201 pore of size Dp, surface area A, and pore volume V = DpA 202 [Fig. 3(a)]. These two phases are in equilibrium with an infinite 203 bulk reservoir which imposes the water and methane chemical 204 potentials µw and µm, and temperature T . For hydrate/liquid 205 equilibrium in confinement, considering the grand potential 206 Ω Φ = -P Φ V + 2γ WΦ A for each confined phase (Φ = L, H), 207 the grand potential equality leads to the Laplace equation:

208 P L -P H = 2 (γ WL -γ WH ) /Dp [1] 209
where the factor 2 accounts for the two surfaces in the slit 210 geometry. γ WH and γ WL are the pore/hydrate and pore/liquid 211 interfacial tensions at the confined melting point (µw, µm, T f ). 212

In confinement, as shown in SI Text, the pressures P Φ at the 213 melting point (µw, µm, T f ) can be expressed using a Taylor 214 expansion around the bulk melting point (µw,0, µm,0, T f,0 ):

D R A F T P H , γ WH µ H w , µ H m P L , γ WL µ L w , µ L m ∆Tf Ω H = Ω L A, V A, V Dp Mineral Hydrate Liquid 1.8 2.0 2.2
215 P H = P0 + 1 v H 0 ∆T f s H 0 + ∆µw + ∆µm x H m,0 1 -x H m,0 P L = P0 + 1 v L 0 ∆T f s L 0 + ∆µw + ∆µm x L m,0 1 -x L m,0 [2]

216

where P Φ = P Φ (µw, µm, T f ) is the pressure of phase Φ at (µw, 217 µm, T f ) and P0 = P H 0 (µw,0, µm,0, T f,0 ) = P L 0 (µw,0, µm,0, T f,0 )

218 is the pressure at (µw,0, µm,0, T f,0 ) [corresponding to Dp ∼ ∞ 219 in the Laplace equation]. ∆T f = T f -T f,0 is the melting 220 temperature shift of confined methane hydrate with respect 221 to its bulk counterpart. ∆µi = ∆µ Φ i = µ Φ i -µ Φ i,0 is the differ- 222 ence of chemical potential µ Φ i of species i [i = methane (m), water (w)] for phase Φ at T f (at this point µ H i = µ L i =
N Φ m,0 /N Φ w,0 = x Φ m,0 /(1 -x Φ m,0
) where x Φ m,0 is the methane mole 237 fraction in phase Φ at (µw,0, µm,0, T f,0 ) and that x L m,0 ∼ 0 for the liquid phase, Eqs. ( 1) and ( 2) lead to:

∆T f T f,0 = v H 0 ∆h f,0 (γ WL -γ WH ) 2 Dp + v L 0 v H 0 -1 P L -P0 + ∆µm v H 0 x H m,0 1 -x H m,0 [3]
where ∆h f,0 = s L f,0 -s H f,0 T f,0 is the bulk molar enthalpy of melting at (T f,0 , P0). At this stage, it is often assumed that (1) the hydrate and liquid molar volumes are similar (v H 0 ∼ v L 0 ) and that (2) the methane chemical potential contribution is negligible (∆µm ∼ 0). With these approximations, one arrives at the Gibbs-Thomson equation which relates ∆T f to Dp:

∆T f T f,0 = 2v H 0 (γ WL -γ WH ) ∆h f,0 1 Dp [4] Introducing Young's equation, γ WL -γ WH = γ LH cos θ with θ = π, Eq. (4) leads to the classical formulation: ∆T f /T f,0 = -2γ LH v H 0 /∆h f,0 Dp.
In contrast to this simplified equation, Eq. ( 3) is a revised version that accounts for the effect of methane gas in the free energy balance and the significant density difference between the liquid and hydrate.

In what follows, we determine the following parameters using independent molecular simulations to check the quantitative validity of Eq. (3): molar volumes v H and v L , molar enthalpy of melting ∆h f = h L w + h V m -h H m,w from the methane hydrate, liquid water, and methane vapor molar enthalpies. Fig. 3(b) shows v H , v L , h H , and h L at different bulk equilibrium conditions: (233 K, 1 atm), (262 K, 10 atm), and (286 K, 100 atm). We obtain ∆h f = 8.35 kJ•mol -1 , v L = 1.85 × 10 -5 m 3 •mol -1 , and v H = 2.28 × 10 -5 m 3 •mol -1 at T = 286 K and P = 100 atm. Such enthalpy of melting ∆h f leads to an entropy of melting ∆s f = ∆h f /T f,0 = 29.3 J•K -1 •mol -1 which is comparable to that reported in [START_REF] Jacobson | Can Amorphous Nuclei Grow Crystalline Clathrates? The Size and Crystallinity of Critical Clathrate Nuclei[END_REF]. In addition to these

324

This result shows that accounting for the methane contribu-

325

tion is important as it represents a non-negligible contribution 326 to the expected Gibbs-Thomson shift in the hydrate melting point. As a last remark, even if the extended equation correctly predicts the melting point depression as a function of pore space D * p available to the fluid molecule center of mass, we note that the classical version of this equation leads to inferred pore sizes that do not strongly differ from the pore size Dp corresponding to the distance between opposite surfaces (from center of solid atoms on opposite walls). This result is important as it shows that the use of the classical Gibbs-Thomson equation to relate pore size and melting point shift remains reasonable -especially for large pores when the difference between Dp and D * p becomes negligible.

Free energy and phase stability. To probe the formation and dissociation kinetics of confined methane hydrate, we combined GCMC simulations with free energy calculations (details are provided in the Methods section). Using the umbrellasampling technique, one can determine the free energy profile Ω(x) as a function of an order parameter x, which describes the transition from the liquid to the crystal phases [Fig. 4(a)]. The local order parameter x = Q6 was used to identify methane hydrate (Q6 ∼ 0.55) and liquid water (Q6 ∼ 0.35). A perfect methane hydrate of dimensions Lx = Ly = Lz ∼ 2.36 nm was prepared for bulk and confined methane hydrate in a slit pore of width Dp ∼ 2.86 nm. Then, a harmonic potential U B (Q6) is added to the energy calculated in the GCMC simulations to force the system to sample states corresponding to the given Q6 value [This strategy is illustrated in Fig. 4(a) using the color squares]. Such a biased potential allows one to determine the probability distribution P B (Q6) from the biased molecular simulations. Finally, the free energy profile is obtained by subtracting the biased potential energy from the biased probability distribution: Ω(Q6) = -kBT ln[P B (Q6)] -U B (Q6).

Fig. 4(b) shows the free energy Ω/kBT for the bulk (I) and confined (II, Dp ∼ 2.86 nm) hydrate/liquid transitions as a function of the local order parameter Q6 at different T . For the bulk transition [I in Fig. 4(b)], methane hydrate is found to be stable for T < T f,0 . For the confined transition [II in Fig. 4(b)], all free energy calculations are performed above the expected melting temperature T f of confined methane hydrate as the umbrella sampling technique at lower T failed to converge (as discussed below, this is due to the slow formation/dissociation kinetics for methane hydrate). As expected, liquid water is the favorable phase at these T ; indeed, free energy difference between methane hydrate and liquid water at these T is negative. To estimate the melting temperature of bulk and confined methane hydrate, the free energy difference ∆Ω as a function of temperature T was extracted from the free energy profiles Ω(Q6) shown in Fig. 4. A polynomial fit was performed to describe each free energy profile and estimate the location and free energy minimum corresponding to the liquid and hydrate phase. As shown in Fig. S9, we find that ∆Ω depends linearly on T with the melting temperature corresponding to the temperature for which ∆Ω ∼ 0. From such data, the bulk and confined melting temperatures, which are consistent with those found using the DCM, are T f,0 = 302 ± 16 K and T f = 257 ± 10 K (the error bars are estimated from the uncertainty in extrapolating ∆Ω ∼ 0).

To gain physical insights into the molecular mechanisms leading to methane hydrate formation/dissociation in confinement, we show in Fig. S10 the density profile of water ρw and methane ρm in the direction normal to the pore surface for [START_REF] Arjun | Unbiased atomistic insight in the competing nucleation mechanisms of methane hydrates[END_REF]. For confined nucleation/dissociation, our microscopic data above -obtained independently using free energy calculations and simple Monte Carlo simulations -suggest that nucleation occurs through direct crystallization as the solid surface leads to significant surface ordering of water which, in turn, promotes hydrate cage formation. However, as described above, the fact that hydrate forms not in direct contact with the solid surface but in the region adjacent to the adsorbed water layer at the pore surface can be seen as parallel to the methane hydrate formation seen for bulk structures. Indeed, the adsorbed water layer coexisting with methane hydrate is playing a role somewhat equivalent to that of the amorphous phase observed under some specific thermodynamic conditions for bulk methane hydrate. In the next section, we develop a simple nucleation model in which crystallization is driven by the free energy balance between surface and volume contributions. Despite its simplicity, this model accounts for the disordered water layer at the pore surface since the interfacial tension between methane hydrate and the porous solid is calculated at coexistence conditions (which specifically include all atomistic/structural details of the confined phases).

Formation and dissociation kinetics. Despite its robustness, umbrella sampling fails to provide a rigorous estimate for the nucleation barrier ∆Ω * because the inferred value depends on the chosen order parameter and finite system size. To apprehend the nucleation process involved in the formation and dissociation of confined methane hydrate, we extend our investigation using a mesoscale description relying on thermodynamic ingredients identified in the molecular approach. In more detail, this coarse-grained strategy relies on the classical nucleation theory which describes nucleation at the mesoscopic level -therefore neglecting microscopic aspects including the role of molecular surface defects on crystallization/melting. Moreover, by relying on simple scalar parameters which are considered as homogeneous and constant, this approach fails to capture complex processes involving an irregular critical nucleus. Yet, as shown in the rest of this section, such a simple model is a robust framework to compare different methane hydrate nucleation mechanisms while relying on thermodynamic reaches the opposite pore surface. Depending on the situation/parameters considered, the bridge grows with a free energy that is either lower or higher than the maximum free energy reached in the hemispherical nucleus configuration. (left) For large pores, the hemispherical nucleus reaches its stability limit before the hydrate nucleus reaches the opposite surface. In this case, the free energy barrier is given by the value reached at the instability point (this mechanism is strictly equivalent to surface nucleation). ingredients derived from our molecular simulation approach.

D R A F T a c R V H A LH A LH R V H q A WH R V H q A WH A LH Bulk Surface

468

As illustrated in Fig. 5(a), considering a porous medium coex-469 isting with a bulk phase made up of gas methane and liquid of the nucleus geometry, its grand free energy can be expressed

487 as Ω = -P H V H -P L V L + γ WH A WH + γ WL A WL + γ LH A LH .

488

The excess grand free energy ∆Ω = Ω -Ω L needed to 489 generate such a nucleus with respect to the liquid phase

490

(Ω L = -P L V + γ LW A) writes:

491 ∆Ω = (P L -P H )V H + (γ WH -γ WL )A WH + γ LH A LH = ∆P V H -γ LH cos θA WH + γ LH A LH

495

For each formation mechanism, the nucleation barrier ∆Ω * 496 corresponding to the critical nucleus is given by the maximum 497 free energy along the nucleation path ∆Ω(V H ) (52).

Fig. 5(c) compares the free energy barriers for bulk, surface, and confined nucleation at a given shift T -T f,0 with respect to the bulk formation temperature T f,0 (all analytical calculations can be found in SI text). As expected, ∆Ω * diverges at T f,0 for bulk and surface nucleation and at T f (Dp) for confined nucleation. The free energy barriers for the bulk and surface mechanisms are ∆Ω * = 4/3πγ LH R * 2 and ∆Ω * = γ LH πR * 2 (2 -3 cos θ + cos 3 θ)/3. As for confined nucleation, the exact critical nucleus and associated free energy barrier depend on T -T f,0 as illustrated in Fig. 5(b). For small T -T f,0 [red line, right of the red square in Fig. 5(c)], the critical nucleus corresponds to a hydrate bridge which forms as the hemispherical cap initiated at the pore surface becomes unstable. For intermediate T -T f,0 [red line, between the red square and circle in Fig. 5(c)], ∆Ω * corresponds to a critical nucleus consisting of the hemispherical cap becoming unstable as it reaches the opposite pore surface. For large T -T f,0 [red line, left of red circle in Fig. 5(c)], the critical nucleus consists of the hemispherical cap becoming unstable before reaching the opposite pore surface (as a result, in this temperature range, confined nucleation merges with surface nucleation). As shown in the inset Fig. 5(c), if we consider ∆Ω * at a given ∆T with respect to T f,0 or T f (Dp), ∆Ω * is much lower for confined nucleation than for bulk and surface nucleation. This result suggests that the hydrate/liquid transition is much faster when confined in porous rocks, which is consistent with previous experimental data [START_REF] Linga | Enhanced rate of gas hydrate formation in a fixed bed column filled with sand compared to a stirred vessel[END_REF][START_REF] Phan | Aqueous Methane in Slit-Shaped Silica Nanopores: High Solubility and Traces of Hydrates[END_REF]. We also compare in Fig. S12(a) the critical nucleus size R * for bulk, surface, and confined nucleation. As expected, R * = -2γ LH /∆P for bulk and surface nucleation are identical as can be inferred from the corresponding free energy approach described in SI text. On the other hand, for all T -T f,0 , R * for confined nucleation is lower or equal to that for surface nucleation.

D R A F T

In line with the discussion above, when the critical nucleus As noted in Ref. [START_REF] Knott | Homogeneous Nucleation of Methane Hydrates: Unrealistic under Realistic Conditions[END_REF], this result suggests that the formation 591

Discussion

A molecular simulation approach is used to show that the thermodynamics of methane hydrate confined down to the nanoscale conforms the classical macroscopic picture as described by the Gibbs-Thomson equation. When confined in pores of a few molecular sizes, even if water and methane possess different wetting properties towards the host surface, they form a gas hydrate phase with structural properties very close to their bulk counterpart. A negative shift in the melting point is observed as a result of the pore/liquid interfacial tension being smaller than the pore/hydrate interfacial tension. Beyond thermodynamic aspects, by showing that confinement leads to facilitated methane hydrate formation and dissociation, our free energy calculations and mesoscopic nucleation model provide a theoretical support for the faster kinetics observed in experimental studies.

Together with the large body of experimental data available, the present results provide a unifying picture of methane hydrate in confined environments or in the vicinity of surfaces. While the description of confined methane hydrate adopted here is simplified compared to real systems, our findings remain meaningful to physical situations such as methane hydrate trapped in rocks on Earth or in other planets, comets, etc. They are also relevant to important energy and environmental aspects such as the expected impact of methane hydrate dissociation on global warming, the use of methane hydrate as energy storage devices, the formation of methane hydrate in pipelines, etc. Despite the simple approach used in this work, in addition to the melting point depression being consistent with available experimental data, the predicted methane hydrate formation and dissociation times provide a robust molecular scale picture of experiments in this field. In particular, considering the mechanisms and times associated with homogeneous and heterogeneous nucleation, methane hydrate is expected to form in the vicinity of pore surface and then to extend to the bulk external phase. In this context, beyond the temperature-dependent crossover between confined and surface nucleation, the predicted scaling between the induction time and pore size τ ∼ Dp is an important result which sheds light on complex experimental behavior. Despite these results, there is a number of aspects which were not considered explicitly in our approach. As shown by Borchardt and coworkers [START_REF] Casco | Experimental evidence of confined methane hydrate in hydrophilic and 842 hydrophobic model carbons[END_REF], even when very similar porous samples are considered, the detailed surface chemistry is expected to play a key role on the thermodynamics and kinetics of methane hydrate formation in confinement. In particular, while our approach only considered the formation (dissociation) of methane hydrate from (towards) liquid water, these authors observed using advanced structural characterization tools more complicated mechanisms involving the formation and distribution of hexagonal and/or disordered ice within the sample porosity. In this regard, pore sizes commensurate with the hydrate structure were selected so that we cannot rule out more complex mechanisms in case of important mismatch between pore space and crystal parameter. Moreover, while diffusion limitations of water and methane through the porosity is expected to affect the kinetics of methane hydrate formation and dissociation (56), they are not taken into account in our statistical mechanics approach which only considers thermodynamic aspects since mass transport and diffusion are not explicitly treated. However, in the classical nucleation D R A F T theory, such diffusion limitations only affect the prefactor and the typical nucleation time remains mostly driven -as a first order approximation -by the free energy metastability barrier.

More importantly, despite the drawbacks identified above, the fact that our molecular simulation approach provides a picture consistent with experiments in almost every point makes our approach suitable to tackle issues related to methane hydrate in confined geometries.

Materials and Methods

Models. Methane is described as a Lennard-Jones (LJ) site with OPLS-UA force field parameters [START_REF] Jorgensen | Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids[END_REF]. Water is described using the TIP4P/Ice water model ( 58

LJ ij (r) = 4ε ij [(σ ij /r) 12 -(σ ij /r) 6 ]
where r is the distance separating i and j while ε ij and σ ij are the corresponding parameters. The LJ interactions are truncated beyond rc = 1.19 nm. The like-atom parameters are given in Table S1 while the unlike-atom parameters are obtained using the Lorentz-Berthelot mixing rules. The intermolecular potential between i and j also includes the Coulomb contribution: u C ij (r) = q i q j /(4πε 0 r)

where q i and q j denote the atomic charges (Table S1). Ewald summation is used to correct for the finite system size with an accuracy 10 -5 . Considering the small box size used here, possible finite size effects cannot be ruled out. All calculations for bulk and confined hydrate were performed for the same lateral extension to minimize the impact of such potential effects. As will be shown below, we stress that the melting point extrapolated to an infinite pore matches accurately the bulk phase transition point as assessed using the direct coexistence method (DCM) and free energy calculations. Considering that finite size effects are expected to affect phase transitions for 2D and 3D systems to a different extent, this result suggests that our comparison between data for confined and bulk methane hydrate is relevant. Possible formation/dissociation hysteresis is another important technical limitation that should be discussed. First, we recall that both the DCM technique and free energy calculations are supposed to be free of such metastability issues. Indeed, by explicitly probing direct coexistence between two phases and/or by assessing accurately the free energy of each phase, these two methods allow overcoming the classical problems in molecular modeling applied to phase transitions (especially, crystallization). This is confirmed by the fact that these two methods lead to consistent results for both bulk and confined methane hydrate. the metastability free energy which can be calculated according to:

Molecular

N k P U k (Q 6 ) exp - U B k (Q 6 ) -Ω k (Q 6 ) k B T [8]

40 that

 40 the shift in melting point qualitatively follows the Gibbs-41 Thomson equation with ∆T f scaling linearly with the recipro-42 cal of the pore size Dp, i.e. ∆T f ∼ 1/Dp [e.g. Refs. (24, 25) 43 and Ref. (20) for a recent review]. This classical macroscopic

44

 44

Fig. 1 .

 1 Fig. 1. Phase equilibrium of methane hydrate-liquid water in nanoporous rocks. (a) Illustration of confined liquid water (cyan) and methane hydrate (blue) in a porous rock described with a stone aspect. (b) Typical molecular configurations for methane hydrate (left) and liquid water (right). The red and white spheres are the oxygen and hydrogen atoms of water. The gray spheres are the methane molecules which are located inside the cages formed by water molecules in methane hydrate. (c) Bulk pressure-temperature (P -T ) phase diagram for ice (I), liquid (L), hydrate (H), and vapor (V). The dashed line denotes the ice melting line while the solid line A-Q-B corresponds to the hydrate melting line: for a given pressure P , H exists at low T while V exists at high T . Confinement is expected to shift the melting point P , T f (P ) to lower T and larger P .

  136

Fig. 2 (DFig. 2 .

 22 Fig. 2(c) shows the methane mole fraction xm for the slit 137

-

  Thomson formalism. Independent molecular simulations were used to determine parameters to assess the quantitative validity of the Gibbs-Thomson equation. Capillary crystallization is described at the macroscopic level using this equation -which is analogous to the Kelvin equation but for liquid/solid transitions. Such thermodynamic model relies on Laplace equation (which links the pressure difference in the confined crystal and liquid to their interfacial tension) combined with the chemical potential equality between the 179 different phases. Two important assumptions are usually made 180 when deriving the Gibbs-Thomson equation: (1) the crystal 181 and liquid have the same molar volume (v C ∼ v L ) and (2) 182 Young's equation holds for liquid/solid systems with a contact 183 angle θ (γ WL -γ WC = γ LC cos θ). While the robustness of the 184 first assumption can be assessed, the second assumption is 185 key as the concept of solid/liquid contact angle in confinement 186 remains unclear. A third assumption, which only pertains 187 to multi-component phases such as hydrates, consists of ne-188 glecting the impact of the lowest mole fraction component 189 on stability. For methane hydrate, this assumption consists 190 of deriving phase stability without considering the chemical 191 potential contribution from methane vapor. In what follows, 192 we formally derive the Gibbs-Thomson equation without in-193 voking the contact angle and extends its applicability to binary 194 solids by including the methane contribution into the stability 195 condition. Then, using independent molecular simulations, we 196 estimate the different ingredients to discuss the validity of the 197 Gibbs-Thomson equation when applied to confined hydrate 198 (all derivation steps can be found in SI Text).

Fig. 3 .

 3 Fig. 3. Gibbs-Thomson effect. (a) L-H-V equilibrium in pores: liquid water (cyan) and hydrate (blue) are confined in a slit pore (gray) of a width Dp, surface area A, and volume V . Ω H and Ω L are the hydrate and water grand potentials at pressures P H and P L . (b) Molar volume v and enthalpy h for liquid water (circles), full-occupancy hydrate (squares), and empty-occupancy hydrate (triangles). These parameters are determined along bulk L-H-V equilibrium at 233 K/1 atm (red), 262 K/10 atm (blue), and 286 K/100 atm (green). The right panel shows the pore/liquid γ WL (black) and pore/hydrate γ WH (gray) interfacial tensions as a function of 1/D * p [D * p = Dp -2σ * is the pore size accessible to the molecule center of mass]. The dashed lines are fitted against γ(Dp) ∼ γ 0 + a exp[-Dp/b]. (c) Melting shift ∆T f /T f,0 for confined hydrate T f with respect to bulk T f,0 at P = 100 atm as a function of 1/D * p . The solid circles are obtained from the direct coexistence method, while the dashed lines are determined using the Gibbs-Thomson equations: the gray line corresponds to the classical equation [Eq. (4)] while the blue line corresponds to the equation in which the methane chemical potential contribution is accounted for [Eq. (3)]. For each equation, the color shaded area denotes the confidence interval (error bars).

Fig. 3 (Fig. 3 (

 33 Fig. 3(b) shows γ WH and γ WL as a function of pore size. Each 272

Fig. 4 .

 4 Fig. 4. Free energy calculations for methane hydrate and liquid water. (a) Umbrella sampling technique: the black solid line shows the free energy Ω as a function of the local order parameter Q6. With this technique, a biased potential energy U B (Q6) (see text) is used to force the system to explore configurations around a given Q6 (as indicated by the color squares corresponding to different Q6). Methane hydrate (Q6 ∼ 0.53) and liquid water (Q6 ∼ 0.38) correspond to the two minima observed in Ω(Q6). The biased simulations provide estimates for the biased probability distribution P B (Q6). The unbiased free energy Ω(Q6) can be computed using the biased distribution minus the biasing potential energy: Ω(Q6) = -k B T ln[P B (Q6)] -U B (Q6). The red and white spheres are the oxygen and hydrogen atoms of water while the gray spheres are the methane molecules. (b) Free energy profile Ω(Q6) for bulk (I) and confined (II) methane hydrates at different T (as labelled in the graph). All free energies are normalized to k B T .

Fig. 5 .

 5 Fig. 5. Nucleation mechanisms for confined hydrate. (a) Schematic nucleus along bulk, surface, and confined nucleation mechanisms. Methane hydrate and liquid are shown as dark and light blue phases while the solid surface is shown in grey. V H and R are the volume and radius of the nucleus while A LH (orange line) and A WH (pink line) are the hydrate/liquid and pore/hydrate interfaces. θ is the contact angle where the liquid, hydrate, and solid intersect. (b) Different schematic grand free energy along the confined nucleation mechanism. In all cases, the nucleation barrier ∆Ω * is shown by the green arrow. The red line shows the free energy change ∆Ω as a function of the volume of the hydrate nucleus V H . Starting from a hemispherical nucleus, three growth paths are possible. (center, right) The nucleus grows and transforms into a bridge nucleus as it

  (c) Nucleation free energy barriers ∆Ω * calculated using the mesoscopic model for bulk (black), surface (blue) and confined (red) nucleation as a function of the temperature shift T -T f,0 with respect to T f,0 . The blue and red dashed lines indicate the formation temperature for bulk (T f ) and confined (T f,0 ) hydrate, respectively. The inset shows the same data ∆Ω * as a function of the temperature shift with respect to the formation temperature: ∆T = T -T f,0 for bulk/surface nucleation and ∆T = T -T f (Dp) for confined nucleation. Dp = 2.86 nm is used for confined nucleation and all calculations are performed with γ LH = 32 mJ/m 2 as estimated in Ref. (54).

  470 water, three possible mechanisms must be considered: (1) bulk 471 nucleation in the external liquid phase, (2) surface nucleation 472 at the external surface of the porous solid, and (3) confined 473 nucleation within the porous solid. Bulk nucleation corre-474 sponds to homogeneous nucleation through the formation of a 475 spherical nucleus which grows until it becomes unstable. Both 476 surface and confined nucleation correspond to heterogeneous 477 mechanisms that initiate at the solid surface. Surface nucle-478 ation involves the formation of a surface hemispherical cap 479 that grows until it becomes unstable at the critical nucleus size. its height reaches the opposite pore surface [Fig. 5(b)]. Using 483 the classical nucleation theory (48, 49), the exact nucleus shape 484 for each mechanism is obtained by minimizing the surface en-485 ergy at fixed methane hydrate volume V H (50, 51). Regardless 486

  492where V = V H + V L and A = A WH + A WL are the total pore 493 volume and surface. The second equality is obtained by using 494 ∆P = P L -P H and Young equation γ WL -γ WH = γ LH cos θ.

539

  of bulk methane hydrate observed in real conditions never 540 occurs through homogeneous nucleation. This interpretation 541 is supported by our data which point to a smaller ∆Ω * for 542 surface nucleation under the same conditions. 543 To provide a quantitative estimate of the forma-544 tion/dissociation kinetics for confined hydrate, let us con-545 sider the experimental situation of a porous medium filled 546 with liquid water and gas methane in equilibrium with a bulk 547 reservoir. We assume that the system is brought to tempera-548 ture T within a time insufficient to allow hydrate formation 549 within the bulk phase (justified as typical experiments are 550 performed on much shorter timescales than the characteristic 551 nucleation time for bulk hydrate). With such conditions, both 552 the confined and bulk phases are metastable when reaching 553 T < T f (Dp). Using the nucleation theory, the nucleation 554 rate 1/τ = J × A for surface and confined nucleation can 555 be expressed as the surface area A developed by the porous 556 material multiplied by the number of nuclei formed per unit 557 of surface and unit of time, J = ρsDZ exp [-∆Ω * /kBT ] (ρs 558 is the nucleus surface density, Ds the diffusivity, and Z the 559 Zeldovich factor which accounts for the probability that a 560 nucleus at the top of the free energy barrier does not melt 561 and actually ends up forming hydrate). Assuming the same 562 prefactor ρsDZ for surface and confined nucleation, the latter 563 expression allows comparing the nucleation times for these 564 two mechanisms from the porous surface area Ap and external 565 surface area Ae. For weak metastability [small T f (Dp) -T ], 566 the confined nucleation time τc is much faster than at the 567 external surface τs because ∆Ω * c < ∆Ω * s and Ap Ae. Sim-568 ilarly, for strong metastability [large T f (Dp) -T ], even if 569 ∆Ω * c = ∆Ω * s , τc τs because Ap Ae. These results are 570 consistent with the work by Page and Sear (55) who investi-571 gated in-pore nucleation by means of Monte Carlo simulations 572 for a 2D spin lattice. These authors reported evidence for 573 a two-step nucleation process where a new thermodynamic 574 phase first nucleates within the pore and then occurs in the 575 bulk external phase. By considering the different mechanism 576 rates, it was also shown that the effective nucleation rate is 577 optimal for a well-defined pore size. Regardless of the ex-578 act confined nucleation mechanisms (hemispherical cap versus 579 bridge), considering that Ap is proportional to the pore surface 580 to volume ratio and, hence, to the reciprocal pore size, the 581 nucleation time is expected to scale as τc ∼ Dp. In particular, 582 the simple model proposed in our work predicts that the ratio 583 between the nucleation times for confined and non-confined 584 hydrate is proportional to the ratio of the external to pore 585 specific surface area, τc/τs ∝ Ae/Ap. Despite its simplicity, 586 this semi-quantitative picture is consistent with experimental 587 observations showing that confinement reduces the formation 588 time from weeks to hours (20, 21). These calculations depend 589 on the choice made for θ but our conclusions remain valid for 590 all angles [Fig. 12(b) in SI].

  ) which contains 4 sites: an LJ site on the oxygen, two point charges on the hydrogen atoms, and a site M corresponding to the O negative charge at a distance d OM = 0.1577 Å from the oxygen toward the hydrogen atoms along the H-O-H angle bisector. In this model, the O-H bond length is 0.9572 Å and the H-O-H angle 104.52 • . As shown in Refs. (41, 59), TIP4P/Ice water combined with OPLS-UA methane accurately reproduces the experimental phase diagram of bulk methane hydrate. We use the stochastic procedure by Buch et al. (60) to generate a configuration with sI structure (27, 41, 60). The slit pores are built as follows. Each pore surface is of lateral dimensions Lx = Ly ∼ 2.38 nm and made up of 11 × 11 squares whose vertices and centers are occupied by a solid atom. Solid atoms are maintained frozen in all simulations. The potential energy includes no intramolecular interactions as we consider a rigid water model, frozen solid, and united-atom methane model. The intermolecular potential between two atoms i and j includes a short-range repulsion and attractive dispersion described using the LJ potential: u

  Dynamics and Monte Carlo simulations. Molecular Dynamics (MD) in the NPT ensemble is used to determine the molar volumes (v L , v H ) and enthalpies (h L w , h V m , h H m,w ) at bulk coexistence conditions (T f,0 , P 0 ). The velocity-Verlet algorithm is used to integrate the equation of motion with an integration timestep of 1 fs. The temperature and pressure are controlled using Nosé-Hoover thermostat and barostat with a typical relaxation time of 2 ps for both NVT and NPT ensembles. All MD simulations were performed using Lammps (61). Monte Carlo simulations in the grand canonical ensemble (GCMC) are used in our DCM to determine the confined transition temperature T f . In this ensemble, the system has a constant volume V , methane and water chemical 722 potentials µm and µw, and temperature T . Monte Carlo moves 723 in the grand canonical ensemble include molecule rotations, trans-724 lations, insertions and deletions for both water and methane. A 725 move from an old (o) to a new (n) microscopic states is accepted or 726 rejected using a Metropolis scheme with an acceptance probability 727 Pacc = min{1, p n µmµw V T /p o µmµw V T } where p µmµw V T corresponds 728 to the density of states in the grand canonical ensemble: 729p µmµw V T (s N ) ∝ V N i = h/ √ 2πm i k B Tthe thermal wavelength for molecule i (m, 731 w). Nw and Nm are the numbers of water and methane molecules, 732 s N is the coordinate set of the N = Nw + Nm molecules in a 733 microscopic configuration whose intermolecular energy is U (s N ). 734 Starting from initial molecular configurations, isobaric-735 isothermal ensemble molecular dynamics is performed for at least 6 736 ns to ensure equilibration is reached. Thermodynamic properties are 737 then averaged using configurations taken in the last 2 ns to estimate 738 the molar volume and enthalpy for methane hydrate, liquid water 739 and methane vapor. As for the GCMC simulations, energy, number 740 of methane and water molecules are averaged over configurations 741 taken every 50000 Monte Carlo moves (where one move consists of a 742 molecule translation, rotation, insertion or deletion). Each GCMC 743 simulation consists of at least 10 10 moves (corresponding to at least 744 10 6 moves per water/methane molecule). Translation, rotation, in-745 sertion and deletion are attempted with the following probabilities: 746 40%, 40%, 10% and 10%. With such numbers, the typical Monte 747 Carlo acceptance probability in the Metropolis scheme is about 35% 748 for translation/rotation and 0.2% for insertion/deletion. While the 749 latter number is rather low, the very large number of Monte Carlo 750 moves considered provide efficient sampling and equilibration (as 751 evidenced by the fact that the DCM calculations lead to results 752 consistent with those obtained using free energy calculations). 753 Free energy calculations. Umbrella sampling is used to determine 754 the free energy Ω as a function of the local order parameter Q 6 . 755 This quantity, which allows identifying liquid water and methane 756 hydrate (62), is determined for a given oxygen O i as: 758where Q 6m,i are complex vectors defined asQ 6m,i = 759 1/N b,i Y 6m (r ij ).The summation over j = 1 to N b,i runs for all 760 neighbor oxygen within a distance r c = 0.35 nm) while Y 6m (r ij ) 761 are the spherical harmonics that depend on the vector r ij = r j -r i . 762The free energy was determined using umbrella sampling with 763 biased grand canonical Monte Carlo simulations. Starting from 764 methane hydrate, we force its transformation into liquid water 765 using a Q 6 dependent biasing potential U B (Q 6 ). Such biased 766 simulations yield the biased probability distribution P B (Q 6 ) which 767 can be corrected from the biased potential energy to determine 768 the free energy as described in the main text. To sample the 769 entire domain Q (0.3 -0.6), we consider Ns = 61 windows with 770 a spacing of 0.05 so that Ns biased Monte Carlo simulations are 771 carried out with the corresponding Q 0 6,k . In more detail, for the 772 k-th window, the following biasing harmonic potential is used: 773 U B k (Q 6 ) = 1/2K[Q 6 -Q 0 6,k ] 2 where K = 5 × 10 7 K is the force 774 constant and Q 0 6,k the central value. We use the weighted average 775 of the unbiased probability distribution of each window P U k (Q 6 ) to 776 determine the full-unbiased probability distribution P U (Q 6 ): 777 P U (Q 6 ) = Ns k=1

  778where P U k (Q 6 ) and N i are the unbiased probability distribution and 779 the number of samples, respectively, in the k-th window.Ω k (Q 6 ) is D R A F T

  the porosity involves a critical nucleus 100 that corresponds either to a hemispherical cap or a bridge 101 nucleus. In both cases, within the classical nucleation theory, 102 we derive a single expression for the nucleation rate 1/τ which 103 is found to scale with the material specific surface area or, 104 equivalently, with the reciprocal pore size, 1/τ ∼ 1/Dp.

		105
	Results	106
	Melting point in confinement. Fig. 2(a) illustrates the 107
	direct coexistence method (DCM) extended to the Grand 108
	Canonical ensemble to investigate the effect of confinement 109
	on L-H-V equilibrium (SI Text). A hydrate molecular con-110
	figuration coexisting with liquid water in a slit pore of a size 111
	Dp is prepared (Dp ∼ 1.67 nm, 2.86 nm, 5.23 nm, and 7.61 112
	nm are considered). Using this starting configuration, for 113
	a given pressure P , Monte Carlo simulations in the Grand 114
	Canonical ensemble (GCMC) are performed for different T to 115
	determine the melting temperature T f as follows. Methane 116
	hydrate melts into liquid water for T > T f [I in Fig. 2(b)] while 117
	liquid water crystallizes into methane hydrate for T < T f [II 118
	in Fig. 2(b)]. The Grand Canonical ensemble ensures that 119
	L-H-V coexistence is simulated de facto (39, 40); because the 120
	system is in equilibrium with an infinite reservoir at chemical 121
	potentials corresponding to bulk liquid water µw(P, T ) and 122
	methane vapor µm(P, T ) at the same P and T , DCM in this 123
	ensemble is equivalent to simulating a three phase coexisting 124
	system [µm(P, T ) and µw(P, T ) are taken from (41)]. It was 125
	verified that for T < T f the methane hydrate occupancy is 126
	equivalent to that for bulk methane hydrate by measuring the 127
	methane mole fraction in the pore center (through integration 128
	of the water/methane density profiles in the region occupied 129
	by the methane hydrate except the water/methane layer in 130
	contact with the surface). Moreover, the confined hydrate 131
	structure was found to correspond to the expected sI structure. 132
	As shown in Fig. S1, the order parameter profile Q6(z) and 133
	hydrogen bonds per water molecule N HB (z) measured along 134
	the position normal to the surface match those obtained for 135
	the sI structure.

mesoscopic description are faster than for bulk hydrate. By 93 considering different nucleation paths (homogeneous nucle-94 ation, surface nucleation, and in-pore nucleation), we identify 95 the critical nucleus leading to the formation of methane hy-96 drate at the surface of the host porous material. This mesoscale 97 thermodynamic approach unravels a crossover in the nucle-98 ation process; depending on temperature, the formation of 99 methane hydrate within

| www.pnas.org/cgi/doi/10.1073/pnas.XXXXXXXXXX Jin and Coasne

Jin and Coasne PNAS | September 28, 2021 | vol. XXX | no. XX | 3

| www.pnas.org/cgi/doi/10.1073/pnas.XXXXXXXXXX Jin and Coasne

| www.pnas.org/cgi/doi/10.1073/pnas.XXXXXXXXXX Jin and Coasne

| www.pnas.org/cgi/doi/10.1073/pnas.XXXXXXXXXX Jin and Coasne

| www.pnas.org/cgi/doi/10.1073/pnas.XXXXXXXXXX Jin and Coasne

Confined nucleation corresponds to the formation and growth 481 of a surface hemispherical cap which transforms into a bridge as consists of the hemispherical cap becoming unstable before 533 reaching the opposite surface, R * for confined nucleation is 534 equal to its value predicted for surface nucleation. For bulk 535 hydrate, in agreement with previous conclusion by Molinero 536 and coworkers [START_REF] Knott | Homogeneous Nucleation of Methane Hydrates: Unrealistic under Realistic Conditions[END_REF], the very large free energy barrier at mod-

782 Starting from Eq. ( 8) with Ω k (Q 6 ) = 0, we iterate self-consistently 783 between Eqs. ( 9) and ( 8) until convergence is reached.

784

While umbrella sampling probes unambiguously the free energy 785 of different phases, they can be prone to biases related to the choice 786 made for the order parameter when estimating free energy barriers.

787

On the one hand, the free energy minima found for each involved 788 phase is independent of the chosen order parameter -therefore