Unbiased atomistic insight in the competing nucleation mechanisms of methane hydrates.

Methane hydrates have important industrial and climate implications, yet their formation via homogeneous nucleation under natural, moderate conditions is poorly understood. Obtaining such understanding could lead to improved control of crystallization, as well as insight into polymorph selection in general, but is hampered by limited experimental resolution. Direct molecular dynamics simulations using atomistic force fields could provide such insight, but are not feasible for moderate undercooling, due to the rare event nature of nucleation. Instead, we harvest ensembles of the rare unbiased nucleation trajectories by employing transition path sampling. We find that with decreasing undercooling the mechanism shifts from amorphous to crystalline polymorph formation. At intermediate temperature the 2 mechanisms compete. Reaction coordinate analysis reveals the amount of a specific methane cage type is crucial for crystallization, while irrelevant for amorphous solids. Polymorph selection is thus governed by kinetic accessibility of the correct cage type and, moreover, occurs at precritical nucleus sizes, apparently against Ostwald's step rule. We argue that these results are still in line with classical nucleation theory. Our findings illuminate how selection between competing methane hydrate polymorphs occurs and might generalize to other hydrates and molecular crystal formation.

Methane hydrate is an ice-like solid substance that encages methane inside a hydrogen-bonded water network. While methane gas is hydrophobic under ambient conditions, a mixture of water and methane gas can spontaneously nucleate into hydrates at low temperatures and high pressures (1). Methane hydrates are abundantly present at the ocean floors and have an estimated double the amount of energy stored compared to conventional resources of fossil fuels (2). With the growing demand of natural gas (3), methane hydrates can act as a potential energy resource (4). On the other hand, the formation of hydrates is of concern to the petroleum industry as they can clog oil pipelines (5), causing flow assurance problems. With this background, interest within the scientific community to comprehend the formation mechanism of hydrates has increased tremendously in the past decade.

The primary theoretical framework to understand hydrate formation is based on the widely used classical nucleation theory (CNT), which postulates a growing spherical solid nucleus in a supercooled liquid phase. As an unfavorable interface is created, this nucleus needs to overcome an activation barrier and reach a critical size before it spontaneously further grows into a bulk crystalline phase. The thermodynamically stable phase of methane hydrate is the structure type I (sI) crystal composed of 2 standard cages, (12 pentagons and 2 hexagons, Fig. 1) and cages in a ratio 3:1 (4). This form is predominant on the Earth’s surface. However, methane hydrates could also form via a nonclassical mechanism (6), in which an amorphous metastable intermediate (composed of nonstandard cages and with many fewer cages) precedes formation of a crystalline state. Alternatively, the less common crystalline sII phase can form, consisting of and cages with ratio 1:2 (4).

Fig. 1.

(A) Schematic view of TPS 1-way shooting (17). An initial trajectory (blue line) crosses the barrier (dashed curve). Trial trajectories either move backward (green line) to the liquid state or forward (red line) to the solid state. A growing nucleus along a path illustrates the nucleation process and is made up of (red) and (blue) cages. (B) Path length distribution for 280 K separates into amorphous (red) and crystalline (blue) paths, suggesting 2 competing nucleation routes. (C) Average path length as a function of temperature. The amorphous paths are ns longer than the crystalline paths at any temperature, indicating a different formation mechanism.

A conceptual understanding of the nucleation mechanism of methane hydrates is still an open question and, in fact, much debated. Current experimental techniques have limited spatiotemporal resolution and cannot easily provide atomistic details about the mechanism of hydrate nucleation. Molecular dynamics simulations can in principle overcome this limitation and provide such insight. However, under moderate conditions, the free energy barrier associated with nucleation is large, resulting in long induction times before the rare nucleation event occurs. This rare event behavior renders straightforward molecular dynamics simulations impractical (e.g., a barrier of would require a wall clock time exceeding the age of the Universe). As the barrier is much reduced when the driving force is higher, e.g., at stronger undercooling or supersaturation, many simulation studies (1, 6–10) were performed at lower temperatures or supersaturated concentration, to observe nucleation within a feasible simulation time. However, recent experimental work has shown that nucleation at high driving force might differ significantly from the CNT prediction (11). Other approaches have used biasing techniques (12) and/or coarse-grained water models (13) to model the methane hydrate nucleation behavior. While all these approaches aim to tackle the rare event problem, they bias the dynamics, which could severely influence the mechanism. For instance, the most common methane hydrate phase in these simulations tends to be amorphous (or sII) and not the experimentally found thermodynamically stable crystalline (sI) phase. This implies that the previously used methods and/or models lack sufficient accuracy, and therefore true, unbiased dynamics using an accurate atomistic force field are required. Moreover, some recent studies (7, 14) indicated that methane hydrate nucleation can follow multiple channels but how this process takes place is yet unexplained.

In this work, we focus on the mechanism of methane hydrate nucleation for moderate undercooling where competition between the formations of crystalline and amorphous solid phases is expected. To avoid all of the abovementioned problems, we employ an accurate atomistic force field (Tip4P/ICE (15)) in combination with the powerful transition path sampling (TPS) (16–18) methodology. TPS harvests an ensemble of unbiased molecular dynamics trajectories connecting the initial reactant and the final product state, by performing a random Monte Carlo walk in trajectory space. By focusing on the rare transition only, the method avoids the induction time involved in the nucleation process, rendering the method exponentially more efficient than straightforward molecular dynamics (MD). It has been stated before (19) that path sampling is necessary to generate a statistically significant set of reactive trajectories to unravel hydrate nucleation. Nevertheless, the cost of even a single trajectory exhibiting a nucleation event is substantial (), and while TPS was previously used to study nucleation processes (20–23), these computational costs have prevented use of TPS for atomistic force fields until recently.

The mechanistic aspects of methane hydrate nucleation are studied in a methane/water mixture with a sI stoichiometry, at experimentally relevant temperatures (24) of 270, 275, 280, and 285 K and a pressure of 500 bar. The melting point for the model used is 303 2 K (25). In the metastable liquid state the system is phase separated: A saturated water–methane mixture is in equilibrium with a methane gas reservoir. Starting from an initial path exhibiting the nucleation transition (), we performed extensive path sampling simulations, amounting to an aggregate simulation time of more than 1 ms of molecular dynamics, and obtained path ensembles for the barrier crossing from the liquid to the solid phase for different temperatures. Note that we label the product state as solid, because a successful pathway can end up in either an amorphous or a crystalline phase.

Fig. 1 visualizes the scheme for path sampling across a free energy barrier. Starting from an initial path, new pathways are generated which reach either the liquid or the solid state. The individual path length fluctuates as each path can take any route in the free energy landscape and is halted when it reaches the solid or liquid state. Depending on the temperature of sampling, methane hydrate nucleation can occur via 1 of 2 main channels. At low temperatures (), the nucleation process is likely to yield an amorphous solid, whereas at higher temperatures (), the critical nucleus becomes larger, and the preferred nucleating mechanism shifts to the formation of a sI crystal phase. Remarkably, at intermediate temperatures (280 K) both mechanisms are sampled in the TPS run. The path length distribution for 280 K in Fig. 1 shows that pathways reaching the amorphous phase are (on average) longer than crystalline paths. Fig. 1 shows that the average duration for a path to reach a solid state from the liquid state increases with decreasing temperature as molecules (especially water) slow down at lower temperature, hampering the ordering of the molecules (). Moreover, at any particular temperature (e.g., 280 K), transition paths that end in the crystalline state are, on average, shorter ( ns) compared to amorphous paths ( ns) (). This faster commitment to the solid phase is also visible as a difference in growth rates ().

Mechanistic Interpretation of the Path Ensembles

The obtained ensemble of pathways contains information about the atomistic mechanism and underlying reaction coordinate (RC) of the processes. To obtain further insight, we need to describe the trajectories in terms of meaningful collective variables (CVs). One of the main advantages of TPS is that it allows direct and unbiased extraction of the reaction coordinate (). The first CV is the nucleus size, expressed as the number of methanes in the growing solid nucleus, also known as the mutually coordinated guest (MCG) (26). The second CV is the cage ratio (CR) defined as the ratio of cages with respect to cages. Both CVs are chosen based on extensive reaction coordinate analysis () which shows that the size of the nucleus (MCG) and the number of cages are key ingredients in the reaction coordinate. The MCG order parameter can distinguish liquid from solid, while the cage ratio distinguishes amorphous (CR 1) from crystal (CR 1) structures.

To interpret the path ensembles, we present in Fig. 2 all of the sampled (extended; ) pathways as a 2D path density plot (27), in the plane of the pertinent variables MCG and CR. (Note that we chose CR above the cages itself, as the latter would make the plots very skewed.) Paths start in the liquid phase at MCG 2, and undergo nucleation until the solid nucleus has become postcritical, and will eventually form a system-spanning structure. At low temperature the path density shows 1 dominant channel (Fig. 2). Most paths (80%) end up in the amorphous phase with a cage ratio (CR 1). The most common cage type observed throughout the trajectories is the for the entire path ensemble, typical for amorphous solid phases. As the system mostly gets trapped in a metastable amorphous structure, the kinetics dominate over thermodynamics at this temperature, and a direct relaxation to the thermodynamically stable state is avoided. This behavior should also be observed for any molecular dynamics simulation of methane hydrate at temperatures below 270 K (28). The amorphous nucleation mechanism might be related to the blob hypothesis (6), which supposes a (metastable) intermediate cluster phase forms before the crystal.

Fig. 2.

Path densities (PD) from TPS as function of CR and MCG cluster size for the 4 temperatures: (A) 270 K, (B) 275 K, (C) 280 K, and (D) 285 K. The white dots overlaying the PD represent the LCP in each path ensemble. The circle diagrams show the overall cage fractions in the LCP, with red color representing cages; blue cages; and gray color nonstandard , , , , and cages. The black vertical dotted line (at CR = 1) demarcates the crystalline and amorphous regions. Both 270 K and 275 K show dominance of cages in the barrier region. At 280 K, the percentage of () cages increases and at 285 K, this cage type becomes most abundant. This visible change in the topology of the barrier region reflects the change in mechanism with increasing temperature.

A few trajectories end in a sII hydrate-type structure, characterized by the presence of more than cages in the last frame of the trajectory. Using Raman spectroscopy, Schicks and Ripmeester (29) showed the sII structure can exist as a kinetic intermediate state before reaching the thermodynamically stable sI type at moderate conditions. We did not observe that rearrangement, due to the long timescales required for the significant reorganization of hydrogen bonds in the solid hydrate.

As the temperature increases, the molecular kinetics become faster (), leading to shorter nucleation trajectories and allowing the system to explore more crystalline structures, and the path density shifts to a higher cage ratio. At the highest temperature of 285 K all pathways are dominated by cages throughout the transition, with CR , and none of the sampled pathways ended in a sII crystal phase. The path density peaks around CR , very different w.r.t. lower temperatures. While the nucleus grows, the monomers ( and cages) arrange themselves in the sI symmetry, a formation mechanism clearly not in agreement with the blob hypothesis. The fast molecular kinetics at this temperature allow easy arrangement in the most stable hydrate form. Indeed, trajectories are now shorter ( ns) (Fig. 1).

At the intermediate temperature of 280 K, the path density is broad (Fig. 2), suggesting that 2 reaction channels coexist at this temperature, 1 toward the amorphous region and 1 toward the crystalline region. Pathways with initially more cages mostly end in an amorphous state while pathways with a higher cage fraction solidify toward a sI methane hydrate (). The path sampling switches back and forth between amorphous and crystalline channels, suggesting the free energy barrier for the 2 routes/channels is roughly comparable at this temperature.

The least changed path (LCP) ( and ref. 17) obtained from the path ensemble is a proxy for the barrier region and hence the critical nucleus. The LCP is plotted as white dots on top of the path densities in Fig. 2. For 270 K the critical nucleus as defined by the LCP has a size between 30 and 40 MCG units, too small to encode much structural information about the solid phase, as there are very few fully formed cages. This can also be inferred from the large degree of discretization in the cage ratio in the LCP, typical for the ratio of small numbers. The size of the critical nucleus grows with temperature to about 120 to 140 MCG units for 285 K.

A structural analysis of the critical nucleus at the different temperatures is given by the circle diagrams in Fig. 2 . At 270 and 275 K, the nucleus is predominantly composed of cages and thus has a strong amorphous character. The circle diagram in Fig. 2 shows that at 285 K the critical nucleus is primarily composed of cages, with CR , indicating that it is already crystalline. Clearly, homogenous nucleation of methane hydrate can directly lead to a thermodynamically stable state, contrary to what was found in previous work (8, 12, 30). At 280 K pathways ending in crystalline and amorphous phases have a similar critical nucleus size, something that has been observed before in the nucleation of ice (31), which indicates that the nucleation barrier for both the amorphous and crystal formations is also located at a similar nucleus size (although at different composition). The LCP switches between and as the dominant cage type. The cage composition in Fig. 2 shows substantially more cages compared to 275 K and, in fact, is bimodal.

This bimodality becomes even clearer when we separately plot the amorphous and crystalline pathways (based on their endpoints) for the subcritical region of the nucleation transition in Fig. 3. Both amorphous and crystalline paths have the same kind of cage ratio distribution for MCG 50 (Fig. 3). As the nucleus grows, the amorphous and crystalline paths diverge and the distribution is bimodal in the critical region (Fig. 3). In the postcritical region, the distribution is uniquely separated (Fig. 3). The path densities in Fig. 3 show that the growing nuclei in the 2 channels branch out in the precritical region. The amorphous paths shift to low CR values, due to an increasing number of cages, before moving back to CR , whereas the crystalline pathways exhibit a steady increase in the amount of cages, resulting in a cage ratio above unity in the critical nucleus regime. Clearly, the cage distributions already diverge in the precritical regime. Such selection might be induced by templating, as forming another becomes easier if there is already a cage with 2 hexagons present.

Fig. 3.

(A and B) Path densities for amorphous (A) and crystalline (B) paths obtained at 280 K. The 2 routes already diverge in the precritical region. (A and B, Left) Typical nuclei along the path are depicted, with cages in red, cages in blue, and methanes in the MCG cluster with incomplete cages as green spheres. (C–E) Cage ratio distributions for specific nucleus size ranges for the amorphous and crystalline trajectories in the 280-K path ensemble. (C) In the precritical domain, MCG 50, both the crystalline and amorphous distributions overlap, indicating that 2 channels are merged in this region. (D) In the critical region, 80 MCG 100, the amorphous and crystalline pathways show unimodal distribution peaking at different values of cage ratio with large overlap. (E) The 2 distributions are separated in the postcritical region, 250 MCG 300, with minor overlap, showing 2 separate channels.

Connection to Classical Nucleation Theory

The critical nucleus size at each temperature can also be estimated by averaging the MCG value of the shooting points in each TPS simulation (Fig. 4 and ). Using CNT we can then estimate the behavior of the nucleation free energy barrier for the formation of the crystalline and amorphous phases. In we show how the driving force can be estimated from the overpressure of methane in the gas phase and can be related to the critical nucleus size in Fig. 4. Together with the previously established surface tension of methane hydrates (32, 33), this in turn determines the CNT free energy barriers, shown in Fig. 4. Note that at 280 K the barriers for amorphous and crystalline nucleation are indistinguishable. Note also that the CNT barriers are generally lower than previous estimates (12), possibly due to the different water model used or the higher driving force induced by the spherical methane gas reservoir.

Fig. 4.

(A) The temperature dependence of the critical nucleus size. Red points are the critical sizes based on the TPS results. The fit to the amorphous phase (blue curve) and the fit to the crystalline phase (orange curve) cross at 280 K. (B) The CNT free energy barriers for amorphous and crystalline nucleation for 270, 275, 280, and 285 K. Note that the barrier height for both amorphous and crystal phases is equal around 280 K, as expected. (C) The findings can be summarized in an idealized CNT free energy surface as a function of size and cage ratio for 275 K (blue), 280 K (red), and 285 K (green). The arrows indicate schematically pathways traversing from liquid to solid (dashed arrows, to amorphous phase; solid arrows, to crystalline phase). While at low temperature (e.g., 275 K), the free energy barrier to nucleate the amorphous solid is lowest, the trend is reversed at higher temperature (e.g., 285 K), where sampled pathways mostly end up in the crystalline phase. At 280 K both mechanisms are accessible.

Fig. 4 combines the above findings in a schematic CNT-based free energy landscape. For low temperature conditions the free energy landscape is tilted toward the amorphous side. When the temperature increases, the free energy landscape tilts to make the barrier at higher crystallinity the preferred one. While we stress that this is only a qualitative analysis and mechanistic details will change with the guest molecule, this switching mechanism may be common in other gas hydrates and water nucleation. Further research using enhanced sampling methods can give a more accurate quantitative estimate.

Discussion

In summary, using extensive transition path sampling simulations (exceeding 1 ms simulation time) we have shown that for decreasing undercooling, homogeneous methane hydrate nucleation switches from a mechanism in which an amorphous solid is formed to a mechanism in which a crystalline solid is formed. In line with CNT, the free energy landscape has 2 underlying channels with 2 nucleation barriers that depend on the temperature. At low temperatures, where the molecular kinetics are slow, the system nucleates into an amorphous intermediate through a rough energy landscape, while at higher temperatures, kinetics is faster and pathways are more guided by thermodynamic driving forces and reach the stable sI structure. This agrees with simulations performed in the NVE ensemble (14), which found enhancement of crystallinity in the nucleus, most likely due to the release of latent heat.

At intermediate temperatures the 2 mechanisms compete, with 1 channel leading to the thermodynamically stable sI hydrate, while the other leads to the amorphous phase or to metastable sII structures. This 2-channel picture is further confirmed by reaction coordinate analysis, which shows that at high temperature, the size and the crystal composition are important, while at low temperature only the size is important. For the high-temperature regime the reaction coordinate close to the dividing surface is characterized by both the nucleus size and the cage type. In contrast, the low-temperature regime shows no dependence on cage type, but only on the nucleus size. Moreover, the optimal dividing surface for the crystallization shows an anticorrelation between the number of large cages and the cluster size. This finding corroborates with the earlier prediction for simple Lennard-Jones systems (20, 23) that crystal quality is a determinant in critical nuclei: Critical nuclei can be small but well formed or large and of lesser quality. The similar trend found here suggests that this is a more general principle in crystal nucleation.

One might think that the slow kinetics at low temperature is the reason that prevents the cage ratio from becoming sufficiently high, but even at low temperature the cage ratios can fluctuate quite a bit in the precritical regime. It is already in this regime that the selection for the final structure is made. Recent experiments on protein crystallization (34) also indicate that polymorph selection can occur precritical. This behavior seems not to follow Ostwald’s step rule (35), which states that a metastable state is formed before the thermodynamically stable state is formed. However, one might argue that all paths, including the crystalline paths, begin in the metastable amorphous state before they shift to the crystalline phase, which would be in line with the step rule.

Our results are consistent with several previously reported hypotheses for methane hydrate nucleation: the labile cluster hypothesis, local structuring, cage adsorption, the blob hypothesis, and the funnel picture. All these hypotheses are aiming at describing the nucleation from different perspectives. Guo et al. (36) proposed the cage adsorption hypothesis in which methane molecules are strongly attracted to hydrate-like cages. As small cages attract methane molecules by a force similar to hydrophobic interaction, these cages should be the first to appear. Indeed, in our path density plots all pathways start in the lower left region. The labile cluster hypothesis (37) states that hydrated cages form and then associate into larger nuclei. The local structuring mechanism (38) pointed out that the guest molecules would have to come together to locally structure the surrounding water. Indeed, we observe only fully formed cages in the core of the growing nucleus. The blob hypothesis (6) is a combination of both these concepts, proposing an initial formation of amorphous aggregate that lacks long-range order of molecules, eventually growing into an amorphous clathrate, which then would transform into a crystalline SI clathrate. We observe this hydrate nucleation mechanism at 270 K and 275 K. However, for 280 K and certainly at 285 K, the amorphous blob is replaced by a more crystalline nucleus even at the initial (precritical) stage (as shown in path density plots). The formation of a blob and subsequently an amorphous phase can be related to kinetic trapping in the rugged nucleation free energy landscape. Kusalik and coworkers (7) have interpreted this behavior in terms of a potential energy funnel model as a conceptual description of molecular crystallization. It is possible to interpret our results in this light, with 2 underlying energy funnels guiding the system in either the amorphous or the crystalline phases.

Our findings contribute to answering the long-pending question about how methane hydrates nucleate at moderate conditions and give a clearer view on polymorph selection in hydrate nucleation. This insight may assist future studies that aim to improve materials synthesis and dissolving strategies of hydrates. Another result of this work is that the common notion of 2-step nucleation (such as the blob hypothesis) is not always applicable to hydrate systems. While the detailed mechanisms depend on the guest molecules, these systems can exhibit a 2-channel nucleation free energy landscape, in which the low temperature 2-step nucleation mechanism is replaced by a direct 1-step crystallization mechanism at higher temperature.

Materials and Methods

System Setup and Force Field.

Water was represented by the TIP4P/Ice model (15) and methane was modeled using united atom Lennard-Jones interactions ( = 1.22927 kJ/mol and = 3.700 Å). This combination was shown previously to mimic experimentally determined properties very well (25). Most of the MD simulations were performed using OpenMM 7.1.1 (39). As an initial configuration, a stable methane hydrate sI structure was generated from the unit cell definition (40) and equilibrated at a pressure of 500 bar (for details see ).

Order Parameters.

For analysis we track various CVs or order parameters (OPs). The nucleus size is measured by the MCG parameter (26), which counts methane molecules involved in the largest solid nucleus in the system, using a home-written code. In addition, we determine several other size-related OPs, such as the number of waters in the nucleus (see for details). The structure of the growing nucleus is characterized by its cage types. We analyze the cage type for each methane in the MCG-based cluster with an algorithm similar to that in ref. 41, using a home-written code. We identified 7 main types of cage structure, , , , , , , and , where the superscript indicates the number of the polygons made by the hydrogen-bonded water in the cage facet. The base number gives the type of polygon (4, square; 5, pentagon; 6, hexagon). The ratio of the number of and cages, denoted the cage ratio (CR), is used as an indicator of sI crystallinity and has been employed previously (10, 12, 14, 42). This cage ratio is for a perfect sI, and lower than unity, , for an amorphous or sII structure.

Transition Path Sampling.

We perform TPS (16), employing the spring shooting algorithm (17) to generate trial trajectories. We define the liquid stable state by the absence of any cluster larger than a dimer, MCG 2. The liquid is in equilibrium with a spherical methane reservoir, formed by methanes that do not dissolve. We define the solid stable state by MCG . For efficiency we halt trajectories at MCG 200, where all paths are postcritical and fully committed toward the solid phase. To determine the final fate of these trajectories for MCG we extended each successful pathway by regular MD until MCG 420, to solidify the entire system. A trajectory is classified as amorphous (or sII), if the transition path in the 250 MCG 300 range has more frames with CR 1 than CR 1, and crystalline otherwise.

An initial trajectory was generated by melting the sI crystal at 400 K and time reversing the trajectory. This unphysical trajectory was brought to the required temperature by shooting trajectories from selected intermediate frames. Two trajectories from the same frame that reached the solid and liquid, respectively, were glued together and used as an initial trajectory for TPS. This path was relaxed subsequently (). All TPS runs were executed using the OpenPathSampling (18) package. All trial paths were stored for later analysis.

Reaction Coordinate Analysis by Likelihood Maximization.

One of the major advantages of TPS is that it allows an unbiased evaluation of the reaction coordinate. We employed the likelihood maximization (LM) approach of Peters and Trout (43) to identify pertinent ingredients in the reaction coordinate. Interpreting the TPS shooting point data as instances of a committor calculation, the LM method predicts the best model to represent those data. Starting from a predefined set of candidate order parameters, the LM method finds linear combinations of these OPs with the highest likelihood to reproduce the shooting point data. For more information see .