Effect of methane concentration on the formation pathways of methane hydrate near hexagonal ice surfaces

Abstract

Molecular dynamics simulations are performed to investigate the heterogeneous nucleation of methane hydrate near ice surfaces over a range of initial methane concentrations (3.0–15.0 mol%) at 250 K and 50 MPa. Across all concentrations studied, the presence of ice enhances methane hydrate nucleation, albeit via distinct mechanisms depending on methane availability. At low initial methane concentrations (<5.8 mol%), ice growth precedes hydrate formation; the advancing ice front concentrates methane in the remaining liquid, triggering hydrate nucleation once the local concentration approaches some threshold (4–6 mol%). The formation of hydrate then consumes methane, enabling further ice growth. Therefore, we observed a coupled growth mechanism regulated by local methane availability. At intermediate concentrations (5.8–6.3 mol%), ice growth is suppressed, yet nucleation remains promoted due to modifications in the nucleation potential near the ice interface, which lower the free-energy barrier and reduce the critical nucleus size. At high concentrations (>6.3 mol%), spontaneous hydrate formation emerges via collective stochastic nucleation, occurring randomly throughout the aqueous phase. In all cases, hydrate nuclei preferentially form at a distance from the ice interface or in the bulk liquid; direct nucleation on the ice surface is rare and observed only in the presence of cubic ice domains arising from stacking faults in the hexagonal lattice. These results reveal the complex interplay between methane concentration, ice growth, and hydrate nucleation, providing mechanistic insights into the dynamic behavior of clathrate formation near ice interfaces.

---

Introduction

Clathrate hydrates, formed under conditions of low temperature and high pressure, are crystalline solids composed of water and small guest molecules. These hydrates exhibit various structures that can accommodate different sizes of guest molecules such as H₂, CH₄, CO₂, C₃H₆, etc. The guest molecules are enclosed within polyhedral cages of water molecules held together by hydrogen bonds.¹ The abundant reserves of clean energy, particularly in the form of methane, stored in naturally occurring gas hydrates, have garnered significant interest from the scientific community.² These hydrates have the potential for diverse technological applications, including cold energy storage,³ CO₂ capture,^4,5 CO₂ sequestration,⁶ desalination,^7,8 seafloor stability, and transportation of natural gases in remote areas.^9,10 The utilization of hydrate-based technologies enhances their value and suitability, especially for countries lacking conventional domestic energy resources.^11,12

Clathrate hydrate formation can proceed via homogeneous or heterogeneous nucleation. Homogeneous nucleation occurs in a uniform system, such as an aqueous methane solution, where nucleation can happen spontaneously under hydrate forming conditions. In contrast, heterogeneous nucleation is facilitated by interfaces or particles that lower the free-energy barrier, providing sites for hydrate embryo formation.^13–15 Experimental studies^16–21 suggest that the homogeneous nucleation rates are often much lower as compared to the heterogeneous ones.^22,23 The reduction in critical nuclei size and energy barriers was considered to be the reason for the rapid nucleation of hydrate, heterogeneously.¹⁹ Therefore, it is believed that gas hydrate formation in nature occurs via a heterogeneous nucleation mechanism. One of the classic examples of gas hydrate formation is in the oil/gas pipelines (flow-assurance problem), causing severe risk of blockages at both onshore and offshore regions.^24,25 Although the formation of hydrates in oil/gas pipelines occurred at the gas–water interface heterogeneously, the solid surfaces (pipeline) also play a significant role. The agitation, mixing (creating more water and gas interfaces), the roughness of the pipe, and the presence of sand or scale can affect the rate of hydrate nucleation. Understanding these interfacial mechanisms and their influence on both nucleation and growth is essential for developing effective hydrate-control strategies and advancing hydrate-based technologies.

A number of studies have been performed to understand the heterogeneous nucleation of gas hydrates.^15,26–37 Different surfaces such as silica, clay, carbon steel, activated carbon, and glass beads were used with a supersaturated aqueous solution of methane to explore the nucleation process in heterogeneous systems.^{15,26,27,29–31,34,38} Babu et al.³⁴ investigated the effect of particle size, pore space, and water saturation (50 and 100%) levels in their experimental study, using silica sand (hydrophilic surface) and activated carbon (hydrophobic surface). In their conclusions, formation and dissociation of a transient layer of hydrate occurred at the surface of activated carbon grains, while in the case of silica sand, hydrate was formed between the interstitial pores of the silica sand particles. A similar kind of experimental study has been performed by Kerkar et al.³⁹ using 5.0 wt% BaCl₂ solutions with glass beads. The author used synchrotron X-ray computed microtomography (CMT) technology to report that methane hydrate was nucleated away from the surface or into the pores of glass beads but not at the surface. Zylyftari et al.⁴⁰ investigated cyclopentane hydrate nucleation in the presence of ice, salt, and oil, suggesting that the ice-oil-aqueous phase contact line was the heterogeneous nucleation site for cyclopentane hydrate formation. Even the presence of an additive (ethanol) in ice powder accelerates the methane hydrate formation process.⁴¹ In the aforementioned experimental studies, the key focus was on where the nucleation occurred in the presence of different surfaces and supersaturated solution. However, experimental observations can be the result of coupled effects of the surface and confinements, making it difficult to isolate the surface effects separately. The deficiency of experiments motivates us to explore the system by molecular dynamics simulation to improve our understanding about the nucleation mechanism, induction time, and nucleation rate calculations by comparing the homogeneous and heterogeneous systems, while using variable guest concentrations.

Several molecular dynamics simulation studies have been conducted to understand the nucleation of gas hydrate in the presence of a surface at the molecular level. S. Liang and P. G. Kusalik¹⁵ studied the methane hydrate nucleation behavior in the presence of hydroxylated silica confinement that was 80 Å apart. The supersaturated aqueous methane solution was prepared with a temperature impulse applied to the initial methane hydrate crystal. The melting of methane hydrate left the amorphous hydrate-like structure in the liquid phase. They observed that the silica surface was able to provide nucleation sites despite the structural mismatch and the appearance of intermediate cages facilitated methane hydrate nucleation at the silica surface. In conclusion, the intermediate half-cage structures (amorphous hydrate-like solid structures) formed initially at the silica surface were able to connect better with the methane hydrate structure. The CO₂ hydrate nucleation in the presence of the silica surface using supersaturated CO₂ concentration (14.81 mol%) was studied by Bai et al.²⁶ They observed an ice-like structure first formed at the silica surface. The nucleation of the first CO₂ hydrate layer developed from the intermediate structure, which that acted as a nucleation seed and formed the second layer. They concluded and suggested a three-stage nucleation mechanism for CO₂ hydrate formation and proposed that the silica surface does not directly help the nucleation of CO₂ hydrate, but the intermediate layer does. Later, Bai et al.²⁹ further explored the effect of hydrophilicity and crystallinity of the silica surface on CO₂ hydrate nucleation. They concluded that the hydrophilicity of a surface can change the nucleation mechanism from a three-step to a two-step by adjusting the local structure of the water layer. Yan and coworkers³¹ performed NVT simulations, using the clay surface (Na-montmorillonite) to study methane hydrate nucleation and growth. Methane molecules migrate and accumulate near the clay pores, which allows the formation of “interlayer hydrate” and the “surface hydrate”. The effects of hydrophobic and hydrophilic surfaces on methane hydrate nucleation were investigated by He et al.,³² using MD simulations. They reported contrasting effects of both the surfaces (graphite and silica) on hydrate nucleation using a supersaturated solution (14.81 mol% CH₄). The graphite (hydrophobic) surface adsorbs the methane molecules and forms a nano-bubble with a flat or a curved (low aqueous CH₄ concentration) surface, and hydrate nucleation was not observed until 2.5 μs. An interfacial water layer has been formed between the nano-bubble and graphite surface, which prevents direct contact between them. In contrast, the silica (hydrophilic) surface helps to form a cylindrical nano-bubble (high aqueous CH₄ concentration) in the solution that helps to nucleate methane hydrate into the bulk but not at the silica surface. A brief summary of the literature work has been given in Table 1.

Table 1 A summary of previous studies on heterogeneous gas hydrate formation and comparison among different systems

Experimental or MD study	Substrate (hydrophilic or hydrophobic)	Guest molecule	T and P	Gas concentration	Surface or bulk nucleation
Note: TIP4P/Ice, six-site, TIP4P, SPC/E, and TIP4P/2005 force field models for water; OPLS-AA and OPLS-UA force field models for methane; EPM2 force field model for CO2; CHARMM, hydroxylated silica model, silanol/silane models, Lopes models for silica surface; CLAYFF model for clay surface; CHARMM27 model for graphite surface.
Experimetal^34 (2013)	Silica sand (hydrophilic)	CH4	277.15 K, 8 MPa	—	Interstitial pores
	Activated carbon (hydrophobic)	CH4	277.15 K, 8 MPa	—	Surface with a transient layer
Experimental^39 (2014)	Glass beads	CH4	272.16 K, 6.76 MPa	14.0%	Bulk
Experimental^40 (2014)	Ice	cyclopentane	273.15, 1 atm	0.05% (v/v)	Ice–oil-aqueous phase contact line
Experimental^41 (2019)	Ice powder	CH4	−15 to −1 °C 16.55 MPa	—	Follows shrinking core model
Experimental^42 (2023)	Active ice	CH4	272 K, 6 MPa	—	Follows shrinking core model
Steady-state MD^15 (2011)	Hydroxylated silica	CH4	250 & 255 K, 50 MPa	8.84%	Surface with intermediate half cage-structure
MD^26 (2011)	Silica	CO2	275 K, 25 MPa	14.81%	Surface with an ice-like intermediate layer
MD^29 (2015)	Silica with varying hydrophilicity	CO2	265 K, 15 MPa	4.2%	Bulk
MD^31 (2016)	Clay surface (Na-montmorillonite)	CH4	260 K, 10 MPa	14.81%	“Interlayer hydrate” and the “surface hydrate”
MD^32 (2017)	Silica	CH4	250 K, 50 MPa	14.81%	Bulk
	Graphite	CH4		14.81%	No nucleation observed
Steady-state MD^14 (2013)	Hexagonal ice	CH4	265 K, 50 MPa, 100 MPa	5.0%, 10.0%	Surface
MD^36 (2017)	Hexagonal ice	CH4	NVE; starting 257 K, 100 MPa	∼12.0%	Surface and bulk
MD^43 (2023)	Hexagonal ice	CO2	NPT 260 K, 5 MPA	19.0%	Whole ice is melted, bulk
MD This work	Hexagonal ice	CH 4	250 K, 50 MPa	3.0–15.0%	Surface and bulk

---

Ice is recognized as a catalytic surface that accelerates gas hydrate nucleation. Experimental work has shown that the presence of an ice interface markedly shortens the induction time for both cyclopentane hydrates⁴⁰ and methane hydrates,^41,42 confirming its promotional role. Ice, which provides a hydrophilic surface, has also been a subject of several studies. Pirzadeh and Kusalik¹⁴ studied methane hydrate formation using steady-state MD simulation, where one ice interface was melting, and the other was growing. In their study, a saturated aqueous solution of methane with 5.0 and 10.0 mol% was used at 50 and 100 MPa. The induced-promote-nucleation (IPN) mechanism was observed, i.e., accumulation of CH₄ molecules near the ice surface and defective cage formation leading to nucleation. Zhang and Guo³⁶ performed NVE simulations to illustrate the nucleation of methane hydrate using hexagonal ice (basal, prism I, and prism II surfaces) with methane as a bubble. Their results indicate that ice melts in the beginning and later hydrate nucleation can occur either homogeneously (into the bulk) or heterogeneously (at the ice surface). The structure of the ice–clathrate interface was studied by Nguyen et al.⁴⁴ using MD simulation. The importance of an interfacial transition layer (ITL), half-cage connecting structures, and the importance of interfacial free energy were explained in detail. The different kinds of partial cages’ formation (5–8 membered, 5–7 and 5–8 membered and 5–6 membered rings)^14,36,44 were reported at the ice interface. Lu et al.⁴³ recently carried out similar simulations of ice, solution, and a CO₂ nanobubble but under NPT conditions where ice melts rapidly and so CO₂ hydrate formed in the bulk solution near the nanobubble. Poon and Peters⁴⁵ developed a stochastic model to explain that a growing interface (ice) can enhance the solution concentration to many folds and significantly accelerated the nucleation process.

Despite extensive experimental and molecular dynamics (MD) simulation studies on heterogeneous gas hydrate nucleation, the detailed relationship between methane concentration, ice morphology, and hydrate formation remains insufficiently understood. Prior investigations have often focused separately on structural motifs or interfacial phenomena, without fully addressing how dynamic changes in methane concentration modulate ice–hydrate interactions at the molecular level. To bridge this gap, we systematically examine the nucleation of methane hydrate in both the presence (heterogeneous) and absence (homogeneous) of ice across a range of initial methane concentrations (3.0 to 15.0 mol%) at 250 K and 50 MPa. Our results demonstrate that ice facilitates methane hydrate formation through distinct mechanisms depending on the methane concentration. Most nucleation events occur near the ice interface rather than directly on it; true surface nucleation is only observed when the growing ice front transitions into cubic ice via stacking faults. These findings provide molecular-level insights into how methane concentration governs heterogeneous nucleation pathways in the presence of a dynamically evolving ice front.

Computational details

The initial models were prepared using Material Studio⁴⁶ as shown in Fig. 1. The perfect ice Ih structure was created by replicating the orthorhombic unit cell in a 6 × 3 × 3 arrangement, resulting in a simulation box containing 432 H₂O molecules with the dimension of 27 Å × 23.4 Å × 22 Å. The construction of ice followed both the Hayward⁴⁶ and Bernal–Fowler ice rule⁴⁷ with zero dipole moment. The ice crystal is then cut open in the z-direction and a slab of liquid water with dissolved methane is added. The overall dimensions of the box (ice + aqueous methane solution) is about 27 Å × 23.4 Å × 124 Å. The basal face of hexagonal ice (Ih) was exposed to the aqueous phase. The basal face of the hexagonal ice is the most stable and slowest growing face among the other faces (primary and secondary prisms) and has lower interfacial free energy (34.5 mJ m⁻²) as compared to the primary and secondary prism faces (35.1 and 35.2 mJ m⁻²).⁴⁸ Similarly, the homogeneous systems were prepared, where only an aqueous methane solution was present with variable methane concentrations. In both homogeneous and heterogeneous systems, a total of 18 models were prepared using variable methane concentrations, as listed in Table 2.

Fig. 1 The initial simulation configurations of (A) the reference homogeneous aqueous methane solution and (B) the heterogeneous system of the solution in contact with the basal surface of hexagonal ice Ih (6 × 3 × 3 lattice, 432 H₂O molecules). Water and methane molecules are presented in their respective colors. All systems are run at 250 K and 50 MPa and a variable initial aqueous methane mole fraction from 3.0 to 15.0%.

Table 2 Description of simulation models (homogeneous and heterogeneous system, total 179 simulation runs)

CH4 concentration (mol%)	Homogeneous system (# of H2O and CH4)	Heterogeneous system [432-ice] (# of H2O and CH4)	Initial box dimensions (Å × Å × Å)	Independent #runs (homo, hetero)
3.0	2554 and 80	2027 and 63	27 × 23.4 × 124	15, 7
4.0	2531 and 105	2041 and 85		15, 5
4.5	2518 and 120	2031 and 96		12, 5
5.0	2508 and 132	2022 and 106		20, 20
5.8	2484 and 153	2005 and 124		20, 20
6.25	2478 and 165	1998 and 133		5, 5
8.33	2428 and 221	1958 and 178		5, 5
10.0	2389 and 265	1926 and 214		5, 5
15.0	2168 and 381	1835 and 323		5, 5

---

Molecular dynamics simulation settings and force fields

All molecular dynamics simulations were performed using the Gromacs (version 4.5.5) and (version 2022.1) software package.⁴⁹ The force fields used in this study are TIP4P/Ice⁵⁰ for water and OPLS-AA⁵¹ for methane, and both have been used in several past studies.^51–53 A spherical cutoff of 10 Å was used for short-range Coulomb interactions and Lennard-Jones potential.⁵⁴ The long-range Coulomb and van der Waals contributions to energy and pressure were calculated by applying Particle–Mesh–Ewald (PME)⁵⁵ and dispersion corrections. The Lorentz–Berthelot (LB) combining rules were used to calculate the parameters of the cross-interaction terms between water and methane. These force fields have been shown to provide good accuracy for various thermodynamic properties.^50,56–58 The initial structures were first energy minimized using a steepest descent algorithm, followed by a short 20 ps NVT run at 200 K to relieve initial stresses and then an NPT heating run at 50 MPa for the next 100 ps from 200 K to the final temperature of 250 K at the rate of 0.5 K ps⁻¹ using a Nose–Hoover thermostat⁵⁹ and Parrinello–Rahman barostat.⁶⁰ MD runs were performed for up to 2000 ns using a leap-frog integrator with a 1 fs time step, saving frames every 5 ps for analysis.

The simulation conditions (temperature, pressure, and aqueous methane concentrations) as shown in Table 2 were carefully selected to permit methane hydrate formation while preventing melting of the ice slab. Preliminary simulations were conducted to establish the equilibrium conditions of relevant systems using the employed force field. At 50 MPa, the melting points of ice and sI methane hydrate were determined to be approximately 270 K and 298 K, respectively, consistent with the experimental values of 268 K⁵⁰ and 290 K.⁵⁷ Thus, the chosen simulation temperature of 250 K lies well below the melting point of both phases (see Table S1 in the ESI†). The equilibrium methane concentration in the liquid–hydrate two-phase system is found to be approximately 2 mol%. Notably, the presence of dissolved methane depresses the melting point of ice by ∼15–20 K. At a methane concentration of ∼5.8 mol%, the ice melting point is reduced to approximately 250 K. Consequently, simulations performed at 250 K and 50 MPa across a range of methane concentrations (3.0–15.0 mol%) allow for the study of the competition between ice growth and hydrate nucleation. The upper bound of 15 mol% was chosen because spontaneous methane bubble formation is observed at this concentration within the first nanosecond of simulation.

Results and discussion

The mechanism of methane hydrate formation varies distinctly with initial methane concentration. Accordingly, we categorize and discuss the formation pathways in three concentration regimes: low (<5.8 mol%), intermediate (5.8 to 6.3 mol%), and high (>6.3 mol%).

Case I: low methane concentrations

For the lower range of initial methane concentrations (3.0 to 5.0 mol%), we present the results using 3.0 mol% methane concentration as shown in Fig. 2 and 5.0 mol% methane as shown in Fig. S1 of the ESI.† These figures provide an overview of the simulation outcomes for heterogeneous systems at low methane concentrations. The overall process of methane hydrate formation can be divided into three stages: Stage I, the pre-nucleation phase; Stage II, during which nucleation occurs; and Stage III, characterized by the growth of the hydrate phase.

Fig. 2 Snapshots taken from the 3.0 mol% methane hydrate nucleation simulation showing the formation of partial cages at the ice interface (A), the appearance of a non-dissolving nucleus at 302 ns (B), the growth of methane hydrate till 600 ns (C), and the final structures of both ice and hydrate at the end of 2000 ns (D). The time evolution of the number of ice (blue) and hydrate-like water molecules (orange) along with the mole fraction of methane (purple) is shown in (E).

Prior to nucleation. As the simulation begins, ice grows rapidly (blue curve in Fig. 2E), resulting in concentrating methane in the aqueous phase (purple curve in Fig. 2E). A layered distribution of methane near the ice interface (∼12 Å) is observed (see Fig. S2, ESI†), consistent with previous studies.^14,36 The locally enriched methane concentration at the ice interface facilitates the formation of hydrate-like partial cages (Fig. 2A), as also reported by Pirzadeh and Kusalik.¹⁴ In all independent runs under varying methane concentrations, transient hydrate-like structures were repeatedly observed at the ice interface. These structures typically had short lifetimes (1–2 ns) and detached from the ice interface due to hydrogen bond breakage. Once separated, they migrated into the bulk liquid as clusters of the hydrate-like water structure, as illustrated in Fig. 2B. These detached structures were not stable. Some partially dissolved in the bulk, while others persisted and acted as precursors or “seeds” for methane hydrate nucleation near the ice interface.

Nucleation stage. As the ice grows, the methane concentration in the surrounding aqueous phase increases, thereby enhancing the thermodynamic driving force for methane hydrate nucleation. Conversely, the elevated methane concentration inhibits further ice growth. The details about the inhibitory effect on ice growth by methane enrichment are given in the ESI† (Fig. S7). As a result, ice growth gradually slows down and nearly ceases when the local methane concentration reaches approximately 4.0 mol% (around 250 ns shown in Fig. 2E). Shortly, methane hydrate nucleation is observed at approximately 302 ns (Fig. 2E). The nucleation time is identified using the MCG1 order parameter,⁶¹ which tracks the number of methane molecules in the largest hydrate-like cluster. In the metastable liquid, this parameter remains low but increases sharply once nucleation begins (see Fig. S3 and S4 in the ESI†). The nucleation time is recorded as the point at which the largest stable cluster first exceeds a threshold of 16 methane molecules. This time agrees well with results based on the F4 order parameter (orange curve in Fig. 2E) and visual inspection of simulation trajectories. Compared to homogeneous systems with the same initial methane concentration, the nucleation rates in the heterogeneous system are significantly higher (see Table 3). This enhancement is attributed to the elevated methane supersaturation resulting from ice growth. Notably, nuclei tend to form near the ice interface (Fig. 2C and Fig. S1, ESI†) but not directly on the ice surface.

Table 3 Summary of the nucleation characteristics of methane hydrate nucleation under various initial aqueous methane concentrations in a homogeneous (without ice) and a heterogeneous (with ice) setting

CH4 conc. (mol%)	BN/SN	Min. distance (Å)	BSA-z (Å)	NR (cm^−3 s^−1) (HeS)	NR (cm^−3 s^−1) (HoS)	Nucleation time (ns) for homo	Nucleation time (ns) for hetero	NR enhancement (HeS/HoS)
Note: bulk nucleation (BN), surface nucleation (SN), minimum distance from ice interface (Min. distance), bulk space available (BSA) for nucleation (z-length), nucleation rate (NR), heterogeneous system (HeS), homogeneous system (HoS); ‘—/—‘ indicates the case of collective stochastic nucleation, where it was difficult to locate the first nuclei formation location. NA stands for not available due to no nucleation within 2000 ns (4 mol%) or 3000 ns (3 mol%).
3.00	7/0	23 ± 9	67 ± 17	4.1 × 10^25	NA	>3000	595 ± 225	NA
4.00	5/0	26 ± 8	83 ± 4	4.2 × 10^25	NA	>2000	441 ± 26	NA
4.50	5/0	31 ± 10	90 ± 8	4.4 × 10^25	3.6 × 10^24	855 ± 82	328 ± 147	12.2
5.00	20/0	25 ± 11	96 ± 4	5.2 × 10^25	1.7 × 10^25	470 ± 51	315 ± 24	3.1
5.80	20/0	24 ± 8	99 ± 1	1.6 × 10^26	4.7 × 10^25	228 ± 32	101 ± 15	3.4
6.25	4/1	35 ± 13	100 ± 0	4.0 × 10^26	2.6 × 10^26	46 ± 14	38 ± 8	1.5
8.33	—/—	—	100 ± 0	5.9 × 10^26	4.5 × 10^26	27 ± 7	26 ± 5	1.3
10.00	—/—	—	100 ± 0	12.2 × 10^26	9.8 × 10^26	13 ± 2	12 ± 4	1.2
15.00	—/—	—	100 ± 0	17.2 × 10^26	12.4 × 10^26	10 ± 2	7 ± 1	1.4

---

Growth stage. Once the nucleus reaches the critical size (MCG1 = 16), it continues to grow irreversibly. Here, a continuous supply of methane molecules is essential to sustain this growth (Fig. 2C and D). As the hydrate nucleus grows, methane is consumed more rapidly than water because sI hydrate contains approximately 14.8 mol% methane, significantly higher than the surrounding aqueous solution. This methane depletion in the liquid phase can relieve the inhibitory effect on ice growth, allowing the ice phase to resume growth. Consequently, methane concentration in the bulk rises again, promoting further hydrate growth. This feedback mechanism produces an oscillatory pattern in methane concentration and a coupled growth behavior between the hydrate and ice phases.

The pace of coupled ice–hydrate growth is set by how quickly methane can be shuffled ahead of the advancing solid fronts. Between 300 ns and 800 ns, the ice front incorporates about 400 water molecules, equivalent to ∼20 Å of new crystalline thickness (22 Å layer contains 432 water molecules). This translates to a linear growth rate of ∼0.041 Å ns⁻¹. As each ice layer forms, methane is rejected into the liquid, causing the bulk CH₄ mole fraction climbs from 4.0 mol% to 6.4 mol%, steadily raising the local supersaturation. Over the same 500 ns window, the hydrate nucleus enlarges from 16 to 36 cages, a gain of 20 cages (≈0.04 cages ns⁻¹). Its radius increases from 9.4 Å to 12.3 Å, giving a radial growth rate of ∼0.006 Å ns⁻¹, roughly one-seventh the rate of the ice front. In contrast, methane diffusion is much faster: with a diffusivity of D = 0.0452 × 10⁻⁵ cm² s⁻¹ (4.52 Ų ns⁻¹), the root-mean-square displacement is , nearly two orders of magnitude faster than the hydrate interface and an order of magnitude faster than the ice front. The result is a coupled growth scenario, in which the advancing ice front acts as a piston that concentrates methane ahead of it, thereby slowing its own growth while promoting hydrate formation. The growing hydrate consumes methane, lowering the local supersaturation and thus the thermodynamic driving force for its continued growth, which in turn allows ice growth to resume.

Notably, two distinct oscillation frequencies in methane concentration are observed in Fig. 2E. A faster oscillation, with a period of approximately 40–50 ns, is associated with the coupled dynamics of ice and hydrate growth and persists throughout the growth phase. A slower oscillation, with a larger period of ∼70–100 ns (from 1100 to 2000 ns), emerges during the later growth stage. Further details about the oscillation of methane concentration are given in the ESI.† This slower fluctuation reflects the progressive depletion of liquid water and increased methane enrichment at both ice and hydrate surfaces, leading to more pronounced variations in measured methane concentration. (The method for determining methane concentration is detailed in the ESI†)

Case II: intermediate methane concentrations

In the intermediate methane concentration range (5.8–6.25 mol%), the growth of ice is effectively suppressed, resulting in a largely stationary ice front. Fig. 3 illustrates the nucleation process of the initial 5.8 mol% scenario. As can be seen, the ice front does not change before the nucleus develops to the size of MCG1 = 16 at 80 ns (comparing Fig. 3A to B). Nevertheless, methane molecules accumulate at the interface, forming a layered distribution pattern as in the low initial concentration cases. (see also Fig. S2, ESI†). The locally elevated methane concentration at the interface facilitates the initial formation of hydrate-like structures directly on the ice surface. However, due to lattice mismatch between the hydrate and the underlying ice crystal, these nascent half-cages are destabilized and tend to detach from the surface. A portion of these detached structures subsequently dissolves into the bulk liquid, while others persist and contribute to the nucleation of stable hydrate clusters near the ice interface. Such structure order formation and transmission are similar in the low initial concentration cases. Kumar et al.⁶² recently showed that in the presence of an ice interface, both the critical nucleus size and the free-energy barrier for nucleation are reduced, thereby facilitating heterogeneous nucleation compared to the homogeneous system. Once the nucleus is formed, the coupled growth mechanism of ice and methane hydrate is observed as in the low initial methane concentration cases (see Fig. S1, ESI†).

Fig. 3 The promoted nucleation of methane hydrate from an initial 5.8 mol% CH₄ aqueous solution. The snapshot taken before the nucleation (A), with nucleus reaching the critical size MCG1 = 16 (B) and coupled growth of ice and hydrate at the end of 300 ns (C).

Case III: higher methane concentrations

At elevated methane concentrations (8.33–15.0 mol%), hydrate-like clusters initiate simultaneously at multiple, spatially uncorrelated sites within the supersaturated aqueous phase, making it infeasible to identify a single dominant nucleation locus. This multi-site emergence of hydrate embryos, often occurring within just a few nanoseconds after the start of the simulation, suggests a collective stochastic nucleation (Fig. 4A) in all runs with ≥8.33 mol% methane (see Fig. S5 (ESI†) for the 10.0 mol% case and Fig. S6 (ESI†) for the 8.3 mol% case). The elevated methane concentration increases supersaturation and thus the thermodynamic driving force Δμ, which reduces the nucleation barrier according to Classical nucleation theory,⁶³

Fig. 4 The case of collective stochastic nucleation using 15.0 mol% methane concentration: (A) multiple nuclei formations are observed at the same time (t = 1 ns) and later (B) convert into the crystalline structure (end of 500 ns). Three layers of ice have been melted, indicating a shrinking core mechanism when methane is excessive.⁶⁴

Therefore, the nucleation process may transition from a rare-event regime to collective stochastic nucleation or even spinodal-like behavior at high methane concentrations. These nascent structures appear randomly throughout the aqueous phase: some located within a few angstroms of the ice interface, and others deep within the bulk liquid. Due to the mismatch of hydrate and ice lattice, there appears to be a quasi-liquid layer (QLL)⁴⁴ between the ice and hydrate nuclei formed near the interface. Therefore, no surface nucleation is observed in these cases, even though the nuclei form very close to the interface. The previous study, where nucleation is observed in the interfacial regime¹⁴ and later separated from the ice surface to form crystal, higher concentration scenario, is consistent with our observations.

At the highest concentration studied (15.0 mol%), which exceeds the sI-hydrate stoichiometry of 14.81 mol%, we also observe methane bubble formation (Fig. 4A) and partial melting of the ice substrate (Fig. 4B) in growth stage. Here, the surplus methane drives further hydrate growth by drawing water from the melting ice, which is consistent with a shrinking-core model under gas supersaturation conditions.^36,64

Nucleation rate calculation. The nucleation rate (J) is calculated using the following equation⁶⁵where N_R is the number of reactive (nucleating) trajectories, N_NR is the number of nonreactive trajectories, τ_i is the induction time for the ith reactive trajectory, τ_j is the total simulation time for the jth nonreactive trajectory, and V_liq is the volume of the aqueous phase. The induction time was calculated using a threshold value of MCG1 = 16, explained in a previous section. In general, the above nucleation rate formula was defined for the homogeneous system. While in the case of the heterogeneous system, the surface area is very important and should be considered. Since the surface area of ice is the same in all the heterogeneous systems in this study, we used the same formula (eqn (2)) to calculate the nucleation rate without normalizing the surface area. This makes the comparison with the homogeneous cases easier.

In general, the nucleation rate falls in the range of 4.1 × 10²⁵ to 17.2 × 10²⁶ cm⁻³ s⁻¹ with the heterogeneous systems having a greater nucleation rate (up to 12 times greater) compared to the corresponding homogeneous system with the same initial methane concentrations. The nucleation rate increases with increasing initial methane concentration for both homogeneous and heterogeneous systems. Our results are similar to the work of Zhang et al.⁶⁶ and DeFever,⁶⁷ where the nucleation rate (5 × 10²⁵ cm⁻³ s⁻¹) was calculated at 250 K and 50 MPa using the MFPT method and 230 K and 50 MPa (1.3 × 10²⁶ cm⁻³ s⁻¹) using forward flux sampling (FFS), respectively.

From Table 3, the nucleation rate of methane hydrate increases significantly in the presence of ice, and this enhancement can be attributed to two primary mechanisms. First, in heterogeneous systems with low initial methane concentrations (3.0 to 5.0 mol%), ice growth during the induction period concentrates methane in the surrounding liquid, thereby increasing the thermodynamic driving force for hydrate nucleation. However, at intermediate initial methane concentrations (5.8 and 6.25 mol%), the ice front remains nearly stationary, yet an enhancement in hydrate nucleation is still observed. This interfacial catalytic effect, beyond mere concentration enhancement, has been shown by Kumar et al.⁶² to involve a modification in the nucleation free-energy landscape near the ice surface, i.e., reducing both the critical nucleus size and free-energy barrier. Notably, this interfacial reduction in the nucleation barrier is expected to contribute across the full range of methane concentrations, complementing the concentration-driven mechanism at lower methane loadings.

Location of first stable nuclei formation. The location of the first non-dissolving nucleus (i.e., the nucleus that eventually develops to form a crystal) is an interesting parameter in heterogeneous systems, which shows the preferred location of the first stable nucleus formation. Fig. 5 shows the probability distribution of such locations at different concentrations (3.0 mol% to 6.25 mol%). It was observed that in most cases, methane hydrates preferred to nucleate near but not at the ice surface. However, there was one case where hydrate nucleation was observed at the ice interface (discussed in the next section). For the case of 8.33 mol% and above, it becomes difficult to identify a single non-dissolving nucleus due to the formation of multiple nuclei at a time (collective stochastic nucleation).

Fig. 5 Histogram (# events) of the separation distance of the first stable nucleus from the ice surface under different initial methane concentrations.

The probability of occurrence of first nuclei formation peaks at 24 ± 4 Å distance and starts to decrease towards or away from the surface. Therefore, it is insightful that the preferred location of first nuclei is at a short distance from the ice interface. We have simulated 7 runs for 3.0 mol% and 5 runs for each concentration of methane (4.0 and 4.5 mol% CH₄), and 20-runs each for 5.0 mol% and 5.8 mol% concentrations of methane. No surface nucleation was observed in 3.0 to 4.5 mol%, while in the case of 6.25 mol%, surface nucleation was observed. Therefore, the nucleation of methane hydrates at the ice interface is a rare event. Fig. 6 illustrates the surface-nucleation mechanism with oxygen atoms colored by their F4 order parameter, calculated from averaged trajectory coordinates. Here, F4  =  −0.95 indicates cubic ice (red), −0.50 indicates hexagonal ice (silver), and +0.98 indicates sI hydrate (blue). The stacking faults within the hexagonal lattice intermittently produce a cubic–ice layer, and surface nucleation (panel (A)) only takes place once a contiguous region of cubic ice is formed, and a sufficient methane concentration exists at the interface. Panels (B) and (C) then capture the ensuing growth of methane hydrate at the cubic–ice interface, with stable connectivity clearly established by 75.35 ns and 100 ns, respectively.

Fig. 6 (A) Surface nucleation case (6.25 mol% CH₄) at 67.25 ns, where cubic ice (red color) is formed to initiate nucleation and (B) and (C) further growth at 75.35 ns onwards and hydrates stick to the ice surface. Hexagonal ice, cubic ice and methane hydrate are shown in this figure by silver, red and blue colors, respectively. Methane gas molecules are shown in green color.

Conclusions

This study systematically investigates how initial methane concentration (3.0–15.0 mol %) influences the nucleation and growth of methane hydrate in the presence of an ice interface at 250 K and 50 MPa. Across the entire concentration range, heterogeneous systems exhibit significantly accelerated nucleation and growth compared to their homogeneous counterparts, underscoring the catalytic role of ice. Three primary mechanisms are identified to account for this enhancement:

1. Methane enrichment by ice growth

At low methane concentrations (<5.8 mol %), advancing ice fronts reject methane into the surrounding aqueous phase, progressively enriching it. As the local methane concentration approaches a critical threshold (∼5.8–6.0 mol %), hydrate nucleation is initiated near the ice interface, coinciding with the slowing or cessation of further ice growth.

2. Thermodynamic barrier reduction

The presence of the ice interface modifies the nucleation free-energy landscape by lowering the nucleation barrier and reducing the critical nucleus size. This facilitation effect persists at all methane concentrations and becomes prominent when the ice front is stationary (intermediate concentrations) and under collective stochastic nucleation conditions at higher methane concentrations (≥8.3 mol %).

3. Coupled ice–hydrate growth

During the growth phase, the ice front and the forming hydrate phase engage in a feedback loop: ice growth increases local methane concentration, which in turn promotes hydrate formation; hydrate growth then depletes methane, allowing the ice front to advance again. This dynamic coupling results in oscillatory supersaturation cycles and synchronized advancement of both phases.

Furthermore, hydrate nuclei consistently emerge within an “active layer” approximately 20–40 Å from the ice interface, rather than directly on the ice surface. This spatial preference arises from the structural mismatch between the hydrate and hexagonal ice lattices, as well as the reduced nucleation barrier in the interfacial region. Direct surface nucleation is rare and is only observed under specific structural conditions, such as the presence of cubic ice domains formed via stacking faults. Overall, this study provides a clearer molecular-level understanding of how methane concentration influences hydrate formation near ice surfaces. These insights have important implications for controlling hydrate formation in both natural environments and engineered systems.

Conflicts of interest

There are no conflicts to declare.

Data availability

The input files for reproducing the data in the work are available at Github at https://github.com/dalipntu/Effect-of-Methane-Concentration-on-the-Formation-Pathways-of-Methane-Hydrate-near-Hexagonal-Ice/tree/main.

Acknowledgements

This research was partially supported by the National Science and Technology Council of Taiwan (NSTC 111-2113-M-001-038-MY2 and 112-2221-E-002-029-MY3) and the National Taiwan University (114L892001). The computational resources from the National Center for High-Performance Computing of Taiwan and the Computing and Information Networking Center of National Taiwan University are acknowledged.

References

1. E. D. Sloan and C. A. Koh, Clathrate Hydrates of Natural Gases, CRC Press, 2008.
2. R. Boswell, Is Gas Hydrate Energy Within Reach?, Science, 2009, 325(5943), 957–958,  DOI:10.1126/science.1175074.
3. Z. Yin, J. Zheng, H. Kim, Y. Seo and P. Linga, Hydrates for cold energy storage and transport: A review, Adv. Appl. Energy, 2021, 2, 100022,  DOI:10.1016/j.adapen.2021.100022.
4. H. Dashti, L. Zhehao Yew and X. Lou, Recent advances in gas hydrate-based CO2 capture, J. Nat. Gas Sci. Eng., 2015, 23, 195–207,  DOI:10.1016/j.jngse.2015.01.033.
5. Z. W. Ma, P. Zhang, H. S. Bao and S. Deng, Review of fundamental properties of CO₂ hydrates and CO₂ capture and separation using hydration method, Renewable Sustainable Energy Rev., 2016, 53, 1273–1302,  DOI:10.1016/j.rser.2015.09.076.
6. J. Zheng, Z. R. Chong, M. F. Qureshi and P. Linga, Carbon Dioxide Sequestration via Gas Hydrates: A Potential Pathway toward Decarbonization, Energy Fuels, 2020, 34(9), 10529–10546,  DOI:10.1021/acs.energyfuels.0c02309.
7. P. Babu, A. Nambiar, T. He, I. A. Karimi, J. D. Lee, P. Englezos and P. Linga, A Review of Clathrate Hydrate Based Desalination To Strengthen Energy–Water Nexus, ACS Sustainable Chem. Eng., 2018, 6(7), 8093–8107,  DOI:10.1021/acssuschemeng.8b01616.
8. M. N. Khan, C. J. Peters and C. A. Koh, Desalination using gas hydrates: The role of crystal nucleation, growth and separation, Desalination, 2019, 468, 114049,  DOI:10.1016/j.desal.2019.06.015.
9. A. K. Sum, C. A. Koh and E. D. Sloan, Clathrate Hydrates: From Laboratory Science to Engineering Practice, Ind. Eng. Chem. Res., 2009, 48(16), 7457–7465.
10. C. A. Koh, E. D. Sloan, A. K. Sum and D. T. Wu, Fundamentals and Applications of Gas Hydrates, Annu. Rev. Chem. Biomol. Eng., 2011, 2(2), 237–257.
11. N. Nakicenovic, N. Grübler and A. McDonald, Global Energy Perspectives, Cambridge University Press, Cambridge, 1998, WEC/IIASA.
12. IEA. Key World Energy Statistics 2011, International Energy Agency, Paris, 2012.
13. S. A. Bagherzadeh, P. Englezos, S. Alavi and J. A. Ripmeester, Influence of Hydrated Silica Surfaces on Interfacial Water in the Presence of Clathrate Hydrate Forming Gases, J. Phys. Chem. C, 2012, 116(47), 24907–24915.
14. P. Pirzadeh and P. G. Kusalik, Molecular Insights into Clathrate Hydrate Nucleation at an Ice-Solution Interface, J. Am. Chem. Soc., 2013, 135(19), 7278–7287,  DOI:10.1021/ja400521e.
15. S. Liang and P. G. Kusalik, The nucleation of gas hydrates near silica surfaces, Can. J. Chem., 2015, 93(8), 791–798,  DOI:10.1139/cjc-2014-0443.
16. P. Thoutam, S. R. Gomari, F. Ahmad and M. Islam, Comparative Analysis of Hydrate Nucleation for Methane and Carbon Dioxide, Molecules, 2019, 24, 6.
17. H. K. Abay and T. M. Svartaas, Multicomponent Gas Hydrate Nucleation: The Effect of the Cooling Rate and Composition, Energy Fuel, 2011, 25(1), 42–51,  DOI:10.1021/ef1011879.
18. K. Lekvam and P. Ruoff, A Reaction Kinetic Mechanism for Methane Hydrate Formation in Liquid Water, J. Am. Chem. Soc., 1993, 115(19), 8565–8569,  DOI:10.1021/ja00072a007.
19. L. Jensen, K. Thomsen and N. von Solms, Propane hydrate nucleation: Experimental investigation and correlation, Chem. Eng. Sci., 2008, 63(12), 3069–3080,  DOI:10.1016/j.ces.2008.03.006.
20. S. Takeya, A. Hori, T. Hondoh and T. Uchida, Freezing-memory effect of water on nucleation of CO₂ hydrate crystals, J. Phys. Chem. B, 2000, 104(17), 4164–4168,  DOI:10.1021/jp993759.
21. N. Maeda and X. D. Shen, Scaling laws for nucleation rates of gas hydrate, Fuel, 2019, 253, 1597–1604,  DOI:10.1016/j.fuel.2019.05.096.
22. S. Garrault-Gauffinet and A. Nonat, Experimental investigation of calcium silicate hydrate (C-S-H) nucleation, J. Cryst. Grow., 1999, 200(3–4), 565–574.
23. B. Kvamme, S. A. Aromada and N. Saeidi, Heterogeneous and homogeneous hydrate nucleation in CO2/water systems, J. Cryst. Grow., 2019, 522, 160–174.
24. H. Fang and M. Duan, Offshore Operation Facilities, Gulf Professional Publishing, 2014, pp. 537–686.
25. H. Fang and M. Duan, Chapter 4 - Special Problems of Deep-Sea Oil and Gas Engineering, In Offshore Operation Facilities, ed. H. Fang and M. Duan, Gulf Professional Publishing, 2014, pp. 537–686.
26. D. S. Bai, G. J. Chen, X. R. Zhang and W. C. Wang, Microsecond Molecular Dynamics Simulations of the Kinetic Pathways of Gas Hydrate Formation from Solid Surfaces, Langmuir, 2011, 27(10), 5961–5967,  DOI:10.1021/la105088b.
27. S. Liang, D. Rozmanov and P. G. Kusalik, Crystal growth simulations of methane hydrates in the presence of silica surfaces, Phys. Chem. Chem. Phys., 2011, 13(44), 19856–19864,  10.1039/c1cp21810g.
28. A. Phan, D. R. Cole and A. Striolo, Aqueous Methane in Slit-Shaped Silica Nanopores: High Solubility and Traces of Hydrates, J. Phys. Chem. C, 2014, 118(9), 4860–4868,  DOI:10.1021/jp500081t.
29. D. S. Bai, G. J. Chen, X. R. Zhang, A. K. Sum and W. C. Wang, How Properties of Solid Surfaces Modulate the Nucleation of Gas Hydrate, Sci. Rep., 2015, 5, 12747,  DOI:10.1038/srep12747.
30. M. E. Casco, J. Silvestre-Albero, A. J. Ramirez-Cuesta, F. Rey, J. L. Jorda, A. Bansode, A. Urakawa, I. Peral, M. Martinez-Escandell and K. Kaneko, et al., Methane hydrate formation in confined nanospace can surpass nature, Nat. Commun., 2015, 6.
31. K. F. Yan, X. S. Li, Z. Y. Chen, Z. M. Xia, C. G. Xu and Z. Q. Zhang, Molecular Dynamics Simulation of the Crystal Nucleation and Growth Behavior of Methane Hydrate in the Presence of the Surface and Nanopores of Porous Sediment, Langmuir, 2016, 32(31), 7975–7984,  DOI:10.1021/acs.langmuir.6b01601.
32. Z. J. He, P. Linga and J. W. Jiang, CH4 Hydrate Formation between Silica and Graphite Surfaces: Insights from Microsecond Molecular Dynamics Simulations, Langmuir, 2017, 33(43), 11956–11967,  DOI:10.1021/acs.langmuir.7b02711.
33. S. Lee, Y. Lee, J. Lee, H. Lee and Y. Seo, Experimental Verification of Methane-Carbon Dioxide Replacement in Natural Gas Hydrates Using a Differential Scanning Calorimeter, Environ. Sci. Technol., 2013, 47(22), 13184–13190,  DOI:10.1021/es403542z.
34. P. Babu, D. Yee, P. Linga, A. Palmer, B. C. Khoo, T. S. Tan and P. Rangsunvigit, Morphology of Methane Hydrate Formation in Porous Media, Energy Fuels, 2013, 27(6), 3364–3372.
35. S. Liang and P. G. Kusalik, Nucleation of Gas Hydrates within Constant Energy Systems, J. Phys. Chem. B, 2013, 117(5), 1403–1410.
36. Z. C. Zhang and G. J. Guo, The effects of ice on methane hydrate nucleation: a microcanonical molecular dynamics study, Phys. Chem. Chem. Phys., 2017, 19(29), 19496–19505.
37. Y. Hu, K. H. Lee, B. R. Lee and A. K. Sum, Gas hydrate formation from high concentration KCl brines at ultra-high pressures, J. Ind. Eng. Chem., 2017, 50, 142–146,  DOI:10.1016/j.jiec.2017.02.007.
38. S. A. Bagherzadeh, P. Englezos, S. Alavi and J. A. Ripmeester, Molecular Modeling of the Dissociation of Methane Hydrate in Contact with a Silica Surface, J. Phys. Chem. B, 2012, 116(10), 3188–3197.
39. P. B. Kerkar, K. Horvat, K. W. Jones and D. Mahajan, Imaging methane hydrates growth dynamics in porous media using synchrotron X-ray computed microtomography, AGU100, 2014, 15, 4759–4768.
40. G. Zylyftari, A. Ahuja and J. E. Morris, Nucleation of cyclopentane hydrate by ice studied by morphology and rheology, Chem. Eng. Sci., 2014, 116, 497–507.
41. Y. A. Chen, L. K. Chu, C. K. Chu, R. Ohmura and L. J. Chen, Synthesis of Methane Hydrate from Ice Powder Accelerated by Doping Ethanol into Methane Gas, Sci. Rep., 2019, 9.
42. P. Xiao, J. J. Li, W. Chen, W. X. Pang, X. W. Peng, Y. Xie, X. H. Wang, C. Deng, C. Y. Sun and B. Liu, et al., Enhanced formation of methane hydrate from active ice with high gas uptake, Nat. Commun., 2023, 14(1), 8068,  DOI:10.1038/s41467-023-43487-6.
43. Y. Lu, Q. P. Li, X. X. Zhang, X. Lv, L. X. Zhang, J. F. Zhao and L. Yang, Molecular Simulation on Hydrate Nucleation in the Presence of Initial Ih Ice and Nanobubble, Energy Fuels, 2023, 37(4), 3307–3313.
44. A. H. Nguyen, M. A. Koc, T. D. Shepherd and V. Molinero, Structure of the Ice-Clathrate Interface, J. Phys. Chem. C, 2015, 119(8), 4104–4117,  DOI:10.1021/jp511749q.
45. G. G. Poon and B. Peters, A Stochastic Model for Nucleation in the Boundary Layer during Solvent Freeze-Concentration, Cryst. Growth Des., 2013, 13(11), 4642–4647,  DOI:10.1021/cg401172z.
46. A. Lottermoser and E. Giese, Measurement of the surface tension of solutions of potassium and lithium salts in fatty acids with the cycle method under exclusion of carbonic acid in air, paying particular attention to further disruptive influences I, Kolloid-Z., 1935, 73(2), 155–174.
47. J. D. Bernal and R. H. Fowler, A theory of water and ionic solution, with particular reference to hydrogen and hydroxyl ions, J. Chem. Phys., 1933, 1(8), 515–548,  DOI:10.1063/1.1749327.
48. M. Ambler, B. Vorselaars, M. P. Allen and D. Quigley, Solid-liquid interfacial free energy of ice Ih, ice Ic, and ice 0 within a mono-atomic model of water via the capillary wave method, J. Chem. Phys., 2017, 146, 7.
49. S. Pronk, S. Páll, R. Schulz, P. Larsson, P. Bjelkmar, R. Apostolov, M. R. Shirts, J. C. Smith, P. M. Kasson and D. van der Spoel, et al., GROMACS 4.5: a high-throughput and highly parallel open source molecular simulation toolkit, Bioinformatics, 2013, 29(7), 845–854.
50. J. L. Abascal, E. Sanz, R. Garcia Fernandez and C. Vega, A potential model for the study of ices and amorphous water: TIP4P/Ice, J. Chem. Phys., 2005, 122(23), 234511.
51. W. L. Jorgensen, D. S. Maxwell and J. TiradoRives, Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids, J. Am. Chem. Soc., 1996, 118(45), 11225–11236.
52. A. M. Fernandez-Fernandez, M. M. Conde, G. Perez-Sanchez, M. Perez-Rodriguez and M. M. Pineiro, Molecular simulation of methane hydrate growth confined into a silica pore, J. Mol. Liq., 2022, 362.
53. J. Y. Wu, L. J. Chen, Y. P. Chen and S. T. Lin, Molecular Dynamics Study on the Growth Mechanism of Methane plus Tetrahydrofuran Mixed Hydrates, J. Phys. Chem. C, 2015, 119(34), 19883–19890,  DOI:10.1021/acs.jpcc.5b05393.
54. J. E. Jones, On the determination of molecular fields - II From the equation of state of a gas, Proc. R. Soc. London, Ser. A, 1924, 106(738), 463–477.
55. A. Y. Toukmaji and J. A. Board, Ewald summation techniques in perspective: A survey, Comput. Phys. Commun., 1996, 95(2–3), 73–92.
56. H. Nada and J. P. J. M. van der Eerden, An intermolecular potential model for the simulation of ice and water near the melting point: A six-site model of H2O, J. Chem. Phys., 2003, 118(16), 7401–7413.
57. P. Gayet, C. Dicharry, G. Marion, A. Graciaa, J. Lachaise and A. Nesterov, Experimental determination of methane hydrate dissociation curve up to 55 MPa by using a small amount of surfactant as hydrate promoter, Chem. Eng. Sci., 2005, 60(21), 5751–5758.
58. L. K. Wang, G. J. Chen, G. H. Han, X. Q. Guo and T. M. Guo, Experimental study on the solubility of natural gas components in water with or without hydrate inhibitor, Fluid Phase Equilib., 2003, 207(1–2), 143–154.
59. W. G. Hoover, Canonical dynamics: Equilibrium phase-space distributions, Phys. Rev. A: At., Mol., Opt. Phys., 1985, 31(3), 1695–1697.
60. M. Parrinello and A. Rahman, Polymorphic Transitions in Single-Crystals - a New Molecular-Dynamics Method, J. Appl. Phys., 1981, 52(12), 7182–7190.
61. B. C. Barnes, G. T. Beckham, D. T. Wu and A. K. Sum, Two-component order parameter for quantifying clathrate hydrate nucleation and growth, J. Chem. Phys., 2014, 140(16), 164506,  DOI:10.1063/1.4871898.
62. D. Kumar, L. Y. Huang, D. T. Wu and S. T. Lin, Accelerated Methane Hydrate Nucleation near Ice: Free Energy Landscape and Kinetic Enhancement, J. Phys. Chem. C, 2025 DOI:10.1021/acs.jpcc.5c03969.
63. D. A. Nicholson and G. C. Rutledge, Analysis of nucleation using mean first-passage time data from molecular dynamics simulation, J. Chem. Phys., 2016, 144(13), 134105,  DOI:10.1063/1.4945256.
64. A. Falenty, A. N. Salamatin and W. F. Kuhs, Kinetics of CO2-Hydrate Formation from Ice Powders: Data Summary and Modeling Extended to Low Temperatures, J. Phys. Chem. C, 2013, 117(16), 8443–8457,  DOI:10.1021/jp310972b.
65. A. Tamhane and D. D. Dunlop, Statistics and Data Analysis: From Elementary to Intermediate, Prentice-Hall, 2000.
66. Z. C. Zhang, C. J. Liu, M. R. Walsh and G. J. Guo, Effects of ensembles on methane hydrate nucleation kinetics, Phys. Chem. Chem. Phys., 2016, 18(23), 15602–15608,  10.1039/c6cp02171a.
67. R. S. DeFever and S. Sarupria, Nucleation mechanism of clathrate hydrates of water-soluble guest molecules, J. Chem. Phys., 2017, 147(20), 204503,  DOI:10.1063/1.4996132.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d5cp02466h

---