Molecular dynamics study of unexpected hydrate formation from supersaturated methane solution at low pressure

doi:10.21203/rs.3.rs-2224894/v1

Download PDF

Article

Molecular dynamics study of unexpected hydrate formation from supersaturated methane solution at low pressure

https://doi.org/10.21203/rs.3.rs-2224894/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Using molecular dynamics the evolution of a metastable solution for "methane + water" was studied for various methane concentrations at T = 270 K and P = 1 bar. We have found the intriguing behavior of the system: the formation of hydrate-like structures is observed in the concentration range. This formation was studied on different time scale levels. We have shown that this process is a collective process of H-bond network reorganization in solution at the initial stage (0–1 ns) until the first fluctuating cavities appear. This is associated with a preceding increase of tetrahedrality the water molecules. The molecules with a high degree of tetrahedrality are distributed uniformly throughout the volume. The formation of small and large fluctuation cavities present in the structure I (sI) appears before steady hydrate growth begins. Later, the constantly present hydrate cavities become the foundation for a few independent hydrate nucleation centers. This new mechanism of hydrogen bond network reorganization depends on the entropy of the cavity arrangement of the guest molecules in the hydrate lattice and leads to the hydrate growth.

Physical sciences/Physics/Chemical physics

Physical sciences/Physics/Condensed matter physics

Gas hydrates are crystalline compounds with hydrogen-bonded water molecules acting as the host network and guest molecules placed in the host cavities. Typically for gas hydrate formation high pressure and low temperature are required. The most common hydrate structures in nature are the cubic structure I (sI), the cubic structure II (sII), and the hexagonal structure H (sH). The formation of the different hydrate structures depends on the size of the guest molecule [1].

Gas hydrate formation is an important process for a number of different applications: gas separation, energy storage, energy transport, cold energy storage, sequestration and desalination. Hydrate formation from the liquid phase involves the crucial step of hydrate nucleation, which is very important for developing efficient hydrate practical applications [2].

The typical time evolution of gas consumption during the hydrate formation process can be divided into three main phases: dissolution phase, metastable supersaturated phase, and growth phase [2]. The dissolution phase starts with increasing gas uptake due to the dissolution into the liquid phase.

Experimental molecular-level description of the hydrate nucleation process in water is challenging not only due to the extremely small length/time scale involved but also due to the stochastic nature of the nucleation process, which is only possible with state-of-the-art measurements. It should also be noted that the anomalous behavior of water is also important. Density fluctuations of water with dimensions of ~ 0.6 nm have been observed [3] and later confirmed by theoretical investigations [4]. The properties of water are determined by its dynamic hydrogen bond network [5], which influences the stochastic nature of the process of hydrate nucleation.

Over the past few decades numerous theoretical molecular dynamics (MD) simulations have been performed in which the authors unraveled crucial aspects of crystal nucleation in liquids (see reviews [2, 6, 7]). These MD simulations have been used to obtain a detailed understanding of the nucleation process in different systems, where the time evolution of a liquid into a crystal is more or less accurately reproduced.

One hypothesis for hydrate formation [1] involves hydrate precursor structures, in which the water molecules form individual clathrate cages (hydrate precursors) around the dissolved guest molecules prior to the hydrate formation process. Then, if favorable conditions are met, these cages will agglomerate and form a solid clathrate hydrate.

Using MD simulations of the hydration of hydrophobic molecules, volumetric numerical data describing details of the process have been obtained. According to the resulting data, the nonpolar solutes change the water hydrogen bond network’s structure, resulting in the creation of structured hydration shells around the dissolved hydrophobic molecules [8–19]. These clusters (often called “postcursors”) are preserved for an extended period of time in the liquid phase [20]. The degree of supercooling influences these shells and more amorphous structures are formed at lower temperatures [21]. Conversely, in other studies [22, 23] the authors did not indicate of enhanced structure around the hydrophobic solutes and did not find clathrate cages in solutions with long lifetimes. Experimental studies were carried out to reveal the water structure around non-polar solutes [24–28]. These results are consistent with previous computer simulation data, where there are loosely organized hydration cages around the hydrophobic solutes, but without evidence of the regular clathrate hydrate cages in solution that are present in water. For example, an NMR investigation of tetrahydrofuran (THF) hydrate formation observed THF hydration shells with no evidence of hydrate precursor formation [29].

There have been many research studies of the mechanism of gas hydrate nucleation [30–36]. Currently there are two main molecular mechanisms for hydrate nucleation. The first is the labile cluster hypothesis (LCH) [30] and the second is the local structure hypothesis (LSH) [31]. According to the LCH, prior to hydrate formation, water molecules form labile, relatively long-lived unstable clusters resembling the solid hydrate cavities. Then, under favorable conditions, these cages agglomerate and recombine into a solid clathrate hydrate phase [30]. Conversely, LSH suggests that guest molecules are arranged in a configuration similar to that of the crystal due to concentration fluctuations, and then the water molecules rearrange around the guest molecules to form a hydrate structure proposed by Radhakrishnan and Trout [31]. The concept of LSH, in which the hydrate can appear first with the ordering of guest molecules followed by water molecules, differs from LCH. Moreover, the authors of LSH argue against LCH by showing that it is thermodynamically favorable for the labile clusters to disintegrate rather than agglomerate to form a larger cluster.

In the hydrate phase the concentration of hydrophobic gas molecules such as methane is more than two orders of magnitude higher than the solubility of methane in liquid water [36]. In order to increase methane solubility it is possible to lower the temperature and/or increase the pressure of the mixture. In one study [37] the authors examined metastable methane solubility using MD simulations. It was shown that to increase methane concentration in the water phase by two orders in magnitude at fixed temperature. It is also necessary to increase the pressure by more than two orders. It follows from the results of another study [38], at 258.5 K and P = 1.2 bar, the equilibrium concentration of methane in water is x_M = 0.0005 molar fraction, and at P = 945 bar, x_M = 0.05. In that study the critical solubility, above which metastable solutions spontaneously form a hydrate phase was found to be about 0.05 molar fraction, comparable to the concentration of methane in the hydrate (x_M = 0.07) with half filling of all the hydrate cavities [38].

The very first stage of clathrate hydrate formation is characterized by gas consumption into the water phase with the formation of metastable supersaturated water. This stage is often called the induction time. Increasing gas concentration is observed until the critical nuclei of a new hydrate phase are formed. A critical nucleus can be defined as the minimum amount of a new phase that is capable of existing independently [38]. In water + gas systems, relaxation processes of the structural redistribution of dissolved guest molecules occur during the induction time, leading to the initial stages of growth phase formation. The metastable supersaturation phase remains a key step toward a complete understanding of hydrate nucleation [39].

Direct production of a supersaturated aqueous hydrate gas solution is still a challenge. However, indirect methods of obtaining supersaturated solutions do exist. Thus, for example, solid amorphous solutions of N2, O2, CO and Ar have been obtained [40, 41], from which subsequently gas hydrates have been produced at unexpectedly low pressures.

Study of the initial steps of hydrate formation [42] showed that a dodecahedral water cage (DWC) immersed in the bulk water phase can efficiently adsorb dissolved methane molecules, which causes the DWC lifetime to increase exponentially with the number of methane molecules adsorbed. At the same time adsorbed methane molecules reduce the number of hydrogen bonds that tend to break the DWC, and change the hydrogen bond network surrounding the DWC from a liquid water structure toward a cage-like hydrate structure. Such behavior argues in favor of considering hydrate nucleation as a cooperative process, where the guest molecules can stabilize the hydrate cages as well as the hydrate structure itself [15, 42]. Another approach to describe the initial evolution of water + gas systems toward hydrate formation was reported by Guo et al. [43], who proposed the cage adsorption hypothesis (CAH). This theory predicts a critical concentration of ~ 0.04 molar fraction for methane aqueous solution, above which methane hydrate begins to nucleate [38].

Several MD techniques and methods have been developed to obtain deeper insight into the mechanisms behind the initial stages of the metastable supersaturation phase and growth of the hydrate phase under a high driving force (high excess pressure and low temperature supercooling). The key parameters found are order parameters [44, 45]. Nucleation pathways associated with both host and guest ordering and the mutual influence of guest and host ordering are being explored. The resulting parameters are used to describe the dynamics and mechanisms of hydrate formation and other key properties.

A number of different simulation studies [46–50] have revealed that nano-bubbles in the gas phase trapped inside the water phase could act as centers of hydrate nucleation at the liquid-gas interface that functions as a heterogeneous nucleation center. Additionally, studies from the Uchida group [51, 52] have shown that the induction time of C₂H₆ and Xe hydrates is significantly lower if the water solution is made using a hydrate with nano-/micro- bubbles that act as nucleation centers.

In MD simulations an applied thermostat dissipates the released energy from the hydrate formation process and restores select thermodynamic properties of the system. Under these conditions, the formation of hydrate cages and the assembly of amorphous hydrate-like clusters that relax over time to more stable uniform or mixtures of hydrate phases have been observed [53].

To produce hydrate structure formation with the presence of a gas-water interface, MD simulation studies are carried out at high pressures (from 500 to 2000 bar) and/or low temperatures [14–19, 23, 32–34, 37, 38, 54–57]. In the referenced works, the nucleation time varied from several hundreds of nanoseconds up to several microseconds depending on the exact thermodynamic (P-T) conditions, selection of guest type, and water molecule model.

A different approach was applied by Sarupria and Debeneditti [58] to study a supersaturated homogeneous aqueous methane solution. The authors simulated a methane solution with 0.07 mol fraction of methane at T = 240 K and P = 200 bar and reported a several hundred nanosecond nucleation time for 0.07 molar concentration of methane, and found structures similar to those reported in the case of a gas-water interface at similar temperature [14–19, 32–34, 37, 38]. The nucleation time is drastically reduced at higher methane concentration (0.08) [59]. Moreover, increasing the degree of subcooling leads to the formation of more ordered hydrate structures. Thus, homogeneous supersaturated methane aqueous solutions make it possible to model hydrate growth under lower overpressure or undercooling [21, 59, 60] in contrast to models with separated phases.

Later simulations showed four phases of nucleation [61]. There is an increase in solvated methane concentration in the aqueous domain via diffusion through the methane − water interface, the formation of unstructured clusters of methane molecules solvated in water, decreasing water content of the nucleus to a value compatible with the type II methane clathrate hydrate composition and final reordering of solvated methane and water as in the ‘blob’ hypothesis [34].

The structures of methane, propane, and methane + propane hydrates formed from water mixtures with high structural order and low computational time have been obtained using MD simulations [59, 62]. In these works, a moderate subcooling level has been applied and overall thermodynamic conditions considered were within the P and T range of practical interest. Other important steps to revealing the hydrate formation mechanism were made by other studies [35, 54, 55, 58, 63] where the initial stages of the metastable supersaturation phase and growth of the hydrate phase at high excess pressure and low temperature supercooling with an aqueous phase prepared by melting the hydrate phase has been studied. The study of hydrate nucleation in a “methane + water + THF” system [64, 65] showed the dependence of induction time, growth rate, and cage structure on the P-T conditions and THF concentration.

We later studied a dynamic hydrogen bond network in the liquid water phase [5, 66]; based on this, we developed efficient methods to analyze hydrogen bond networks and a method for determining the hydrate nucleus. In previous work [67] we showed the possibility of methane hydrate growth from a bulk water + methane solution or in the presence of sea salt and hydrate seed at excess pressure without notable induction time. However, the nucleation mechanism was not studied.

In the present work, we have focused on an examination of methane hydrate growth outside its thermodynamic stability region (without excess pressure and/or supercooling applied) using MD simulation methods from a supersaturated methane solution. The initial solutions are based on a uniform distribution of gas molecules inside the bulk water phase, which results in an ideally stirred solution. The aim is to study the dynamics of an H-bond network and the influence of the fluctuation of hydrate nuclei on the hydrate growth process in supersaturated methane solution with an emphasis on studying the behavior of water molecules with a high degree of local order.

Five homogeneous methane + water mixture models with different compositions (see Table 1) were constructed. The selected methane concentrations (C) correspond to 20%, 40%, 60%, 80% and 100% rates of single filling of small (5¹²) and large (5¹²6²) cavities of sI hydrate composed of 10001 water molecules. In contrast with previous works [35, 54, 55, 58, 63], where the aqueous phase was (homogeneously) supersaturated with gas and prepared by melting the hydrate phase, in the current study the initial supersaturated solution was similar to a different study [59]. The water molecules were placed at the nodes of a virtual cubic lattice with a step of 3 Å and a random orientation in space. The number of lattice sites slightly exceeded the number of water molecules, so some sites remained empty, but the model remained cubic. The gas molecules were also randomly arranged at the nodes of another independent virtual lattice with the same step size. For verification purposes, two series of simulations (A and B) were carried out with different distributions of water and gas molecules. By default, we present the results corresponding to series B, and qualitative differences between series A and B are described. Before initiating the simulation, molecular non-overlap was confirmed.

Table 1

Designation of concentration, C, for methane (N_M) + water (N_W) solution compositions and correspondence to molar concentration (x_M).
C	N_M	N_W	x_M
0.2	348	10001	3.36 mol%
0.4	696	10001	6.5 mol%
0.6	1043	10001	9.45 mol%
0.8	1391	10001	12.2 mol%
1.0	1739	10001	14.8 mol%

All simulations were conducted at T = 270 K and P = 1 bar in the NPT-ensemble with 1 fs time-step using LAMMPS software [68]. The time constants of the Nose-Hoover thermostat and barostat [69, 70] implemented in LAMMPS were 300 fs and 3000 fs, respectively. Details of selected parameters are presented in Supplementary Materials. Long-range interaction was calculated using the PPPM method [71]. The cut-off distance of the Coulomb and van der Waals forces was set to 13 Å. The SHAKE algorithm [72] was used in order to maintain the geometry of the water molecules. Periodic boundary conditions were applied for every model.

Snapshots for analysis purposes were written every 100 ps for all long-term simulations (more than 10 ns). In order to study the pre-nucleation processes and early steps of hydrate nucleation itself for the large 0.6 model (10001 H₂O + 1043 CH₄) during the first 10 ns, snapshots were taken every 100 fs. In order to reveal structure changes due to nucleation, the pure water model containing 10001 H₂O molecules and the bulk methane hydrate model containing 6 x 6 x 6 sI unit cells (9936 H₂O + 1728 CH₄) were simulated at the same conditions over 10 ns and 1 ns, respectively.

The TIP4P/Ice potential [73] (σ = 3.1668 Å, ε = 0.21095 kcal/mol, q_O = − 1.1794 |e|, q_D = + 0.5897 |e|, O-M distance = 0.1577 Å, O-H distance = 0.9572 Å, and H-O-H angle = 104.52°) was used to describe the water molecules due to its high melting point accuracy in comparison to other simple 3-point and 4-point potentials. Water molecules were considered as rigid. Methane molecules were considered as single-point particles with potential parameters described elsewhere [74].

The structures of the homogeneous mixture (series A and B) and the simulation run parameters with force-field data are presented in Supplementary Materials.

The analysis of the hydrogen bond (H-bond) network was carried out as described in our previous works [5, 66]. The geometric criterion for H-bonding was chosen as follows: R_OO ≤ 3.2 Å and ∠HOH ≤ 30°. The local (short range order) tetrahedrality order parameter F₃ [44] and torsion angle order parameter F₄ [45] were used in this study. Reference points of the F₃ parameter for the ideal water crystal and bulk liquid are F_{3_crystal} = 0 and F_{3_liquid} ≈ 0.1. Reference points of the F₄ parameter for the hydrate, bulk liquid and ice are F_{4_hydrate} ≈ 0.7, F_{4_liquid} ≈ 0 and F_{4_ice} ≈ − 0.3.

Analysis of polygons (the closed chains of 4, 5, 6 H₂O molecules connected by H-bonds) for the presence of common edges made it possible to find various cavities (closed polyhedrons). The hydrate with structure I (sI) consists of 5¹² and 5¹²6² cavities only, however the 5¹²6³, 5¹²6⁴, 4¹5¹⁰6², 4¹5¹⁰6³ and 4¹5¹⁰6⁴ cavities were studied as well.

Fluctuation character of the H-bond network structure at the hydrate pre-nucleation interval

The significant effect of guest molecules on the properties of the H-bond network leading to the formation of small 5¹² fluctuation cavities at the initial stage of the hydrate nucleation process in a methane solution at low pressure has been shown.

The characteristic time dependence of the number of large (5¹²6²) and small (5¹²) sI hydrate cavities in solution with C = 0.6 and P = 1 bar is presented in Fig. 1 for a single trajectory at which the formation and dissipation of sI cavities were observed over the time period 0-10 ns. The simulation could be divided into several intervals according to the features of small and large cavity formation. In the first interval (up to 1 ns) no cavities were formed. The second interval, 1-4.5 ns, is characterized by the formation and dissipation of single hydrate cavities (fluctuation cavities) in water solution (see Fig. 1b) with lifetimes of about one picosecond and which could be the predecessors of hydrate nuclei. In the next interval (4.5-6 ns) at least one cavity exists in the volume but because of fluctuations it is possible to observe up to 3 different separate small cavities (see Fig. 1c). Then after 6 ns, transformation of the hydrogen bond network is observed, which leads to the formation of 2-3-4 small cavities as shown in Fig. 1d, and only after a simulation time of 7 ns, the fluctuation of large sI cavities becomes detectable and then persistent after 7.5 ns, although their number fluctuates (see Fig. 1e). After a 10 ns simulation time, steady growth of the hydrate is observed. The number of cavities grows, and fluctuations do not lead to the disappearance of a significant fraction of the cavities. During modeling of system transformations there are no particular spatial ordering of the cavities and the entire hydrogen bond network of water molecules is involved. In this simulation we focused only on small (5¹²) and large (5¹²6²) cavities that form the sI hydrate, but as noted earlier other cavity types are present in this system and show similar behavior.

Figures 2a,b shows the time dependence of the order parameter F₃ and F₄ in the first 10 ns of simulation and comparisons to the bulk hydrate, bulk water and cavity forming water molecules in solution. The difference between F₃ and F₄ for bulk water and F_{3_liquid} and F_{4_liquid} is due to the fact that the TIP4P model series does not have an H-O-H angle equal to the right tetrahedral angle of 109.47 degrees, which is used in some other potentials. The value of the order parameters for bulk water does not change over time because the water is in an equilibrium state. The sharp decrease of the F₃ parameter (an increase in tetrahedrality compared to pure water) occurs during the first 1 ns, which shows steady H-bond network reorganization (the difference is ~0.01 and gradually increases). Qualitatively similar results are presented for the F₄ local torsion angle order parameter: whole solution order is close to pure water; however, the order parameters of cavity forming molecules in solution for both F₃ and F₄ are very close to those of a crystal in the same temperature and pressure region. Thus, the scattered cavities formed in the H-bond network are topologically close to hydrates. In general, the presence of dissolved methane molecules has an ordering effect on the structure of water. This water structure ordering in the presence of high concentration of guest molecules is similar to the behavior of water molecules described in ref. [75] where clathrate-like structures were shown to form under ambient conditions in the presence of confinement and a hydrophobic environment. Authors noted that as hydrophobic compounds, many different molecules can be chosen, such as alkanes like methane, ethane, propane, etc.

Figure 2c shows the time dependence of the portion of highly ordered water molecules for which the parameter F₃ is lower than 0.025 (highly ordered water molecules, HOWMs). This value was chosen because the water molecules forming the cavities have an average value of F₃ < 0.025. As expected, the order parameters of these molecules are much more "crystalline" than for other molecules. Increase in the number of HOWMs compared to bulk water confirms the ordering of the structure prior to the formation of regular hydrate cavities. The fluctuations in the number of HOWMs indicate the dynamic nature of the local environment of all molecules united by the H-bond network.

Figure 3 shows the visualization of HOWMs during the first nanoseconds. At t = 0 ns (Fig 3a) the number of HOWMs is insignificant; however, at t = 0.1 ns (Fig 3b) this number is about ~20-25% of all molecules. It is important to note that HOWMs are delocalized in this time period (Fig 3a-f) and do not form a single local cluster. This allows the cavity to fluctuate forms by any group of H-bond network HOWMs.

Radial distribution function of HOWMs.

The radial distribution functions (RDF) were calculated to characterize the reorganization of the H-bond network structure. Figure 4a-d shows the oxygen–oxygen RDF for HOWMs in methane solution, bulk water and bulk hydrate systems. The solution model peaks transform from liquid type with 3 peaks to hydrate type with 4 peaks over the time period. The position of the third peak shifts closer to the hydrate type. The high value of the first peak for the solution could indicate that HOWMs form groups being uniformly distributed and form hydrate-like structure/fragments. At a larger time scale, the further convergence of the RDF lines for the HOWMs and the general structure of water can be expected. The gas-oxygen RDF is discussed in Supplementary Materials and is presented in Fig. S3.1.

Hydrate growth outside the stability region

Using the MD method, we examined the evolution of two sets of independently generated water + methane systems with different methane concentrations in the system at P = 1 bar and T = 270 K. The total number of water molecules was 10001. Chosen concentrations are much higher than the equilibrium methane concentration in water, and thus the supersaturated homogeneous solutions obtained here are metastable. The pressure and temperate fluctuations are presented in Fig. S4.1a,b and Fig. S4.2a,b in Supplementary Materials. The simulation for C = 0.2 over 100 ns (Fig. 5a-d) shows neither hydrate formation nor phase separation but a metastable supersaturated solution with homogeneously distributed gas molecules. At methane concentrations C = 0.4 (Fig. 5e-h) and C = 0.6 (Fig. 5i-l), the simulation results in the appearance of a volume region with a different methane concentration. The region with higher methane concentration corresponds to a hydrate phase. The other region is a water solution with greatly decreased methane concentration in comparison with the initial supersaturated phase. Increasing concentration (C = 0.8, Fig. 5m-p) leads to the appearance of the third region of high methane concentration – gas nanobubble. Formation of a gas phase at C = 0.8 is a competing process with the hydrate formation process: unlike the other systems where all runs provide similar results, for C = 0.8 in the first run behavior like C = 1.0 is observed, whereas in the second run we observed a hydrate growth process but with a nano-bubble formation. At C = 1.0 (Fig. 5q-t), no hydrate formation was observed. Complete 100 ns evolutions of metastable supersaturated solutions with C = 0.2, 0.4, 0.6, 0.8, 1.0 are presented in Supplementary Videos S1, S2, S3, S4, S5 respectively. According to a previous study [37] nanobubbles formed in the water phase at C = 0.8 and C = 1.0 could be considered as hydrate nucleation centers at the microsecond timescale [76].

Experimental obtaining of supersaturated liquid solutions at atmospheric pressure is a task that does not have a direct solution, but the formation of a solid solution seems to be possible through the deposition of water vapour and methane on a cooled substrate [40,41].

The chemical potential of water molecules in the hydrate phase depends not only on the interaction between water molecules in the hydrate lattice, but also on the number of the guest molecules in the hydrate cavities and exists as an entropy contribution as shown in our previous work [77]. The value of this contribution at a certain guest molecule concentration exceeds the difference between the chemical potential of water molecules in the liquid phase and the empty hydrate lattice framework. This leads to hydrate stabilization. In this case, due to a sufficiently large number of dissolved methane molecules, the filling of the hydrate is large enough to stabilize it with respect to the liquid phase even at atmospheric pressure. We observed this phenomenon at concentrations of 0.4 (6.5 mol%) and 0.6 (9.45 mol%) during MD simulations. When the concentration increases to more than 0.8 (12.2 mol%), phase separation occurs and this is in agreement with theoretical work [59]. This leads to a significant decrease in methane concentration in the water phase and disruption of this hydrate formation mechanism. With a decrease to 0.2 (3.36 mol%) or lower methane concentration, the amount of dissolved methane molecules in the water phase is insufficient for an adequately large entropy contribution to the water chemical potential, which was confirmed by the presence of a minimum methane concentration in the solution to initiate hydrate formation [38].

Molecular dynamics simulations using the TIP4P/Ice water model [54] showed the growth of methane hydrate at T = 245 K and P = 0.2 bar, which is outside the methane hydrate stability area [78]. The authors hypothesized that unexpected hydrate growth is due to the spherical shape of the initially separated gas phase, which can create additional pressure and cause stable hydrate growth [54].

Previous results [21,59,60] assume that using a supersaturated methane solution as the initial structure allows one to observe hydrate formation at moderate conditions such as lower pressure and/or lower sub-cooling. Our results show that the effect of supersaturation of an aqueous solution of methane has greater importance than was previously suggested and allows hydrate formation at 1 bar.

The dependence of the number of cavities on time for systems with different methane concentrations is presented in Fig. 6. Methane forms an sI hydrate that consists only of 5¹² (Fig. 6a) and 5¹²6² (Fig. 6b) cavities. Nevertheless, during the hydrate formation process other cavity types could be formed because of fluctuations, surface effects and other effects. These cavity types were also drawn in the plots (Fig. 6c-g). For systems with concentrations of 0.4, 0.6, and 0.8, the number of observed cavities increases with time. It should be noted that the formation of large (5¹²6²) cavities occurs after the onset of formation of small (5¹²) ones. An increase in the number of different types of cavities is strongly correlated with methane concentration in the water phase. The system with C = 0.6 shows the fastest speed and, therefore, this concentration is the most favorable for hydrate formation under these conditions. We believe that at different pressure/temperature conditions, the most favorable methane concentration should be nearly the same. In other words, hydrate formation is faster if the gas/water molar relation is within a certain range. Similar dependencies for series A and several smaller models are given in Supplementary Materials (Figures S6.1 and S6.2).

The sI hydrate unit cell consists of 6 large and 2 small cavities, but in our calculations, this relation is not met due to a large number of other cavity types (5¹²6³, 5¹²6⁴, 4¹5¹⁰6², 4¹5¹⁰6³, and 4¹5¹⁰6⁴) that reconfigure in later hydrate formation phases. Thus, the 4¹5¹⁰6² cavity topologically is very close to the large 5¹²6² cavity present in the sI hydrate and can easily reconfigure into a proper large 5¹²6² cavity. Using such an assumption the resulting relation is not far from the target 3/1 value. Disintegration (change of structure/size) of the cavities is observed only as fluctuations in time. In the time interval of 0–100 ns, growth of the hydrate phase (0.4) reaching saturation (0.6) was observed. The case of C = 0.8 is exceptional, where slow hydrate decomposition is caused by gas phase separation. However, the resulting phases (0.4 and 0.6) are not thermodynamically stable but under the given p-T conditions, they are metastable. Additional driving forces (such as supercooling, high pressure, promoters, or a helper gas) could be used to stabilize the hydrate structure.

The dependence of the normalized H-bond number on time and concentration is shown in Fig. 7a and confirms previous results. The growth of N_H-bond/N_mol indicates an increasing number of molecules involved in the formation of 4 H-bonds that are connected with crystal-like ordering, where each molecule forms exactly 2 donor and 2 acceptor H-bonds (N_H-bond/N_mol = 2). At C = 0.2 the local structure of the solution does not change in terms of mean values. The decrease in N_H-bond/N_mol at C = 1.0 clearly indicates that a high proportion of molecules form the phase boundary. The N_H-bond/N_mol steadily increases for solutions with initial C = 0.4 and 0.6, and in the case of C = 0.8, during the first 25 ns (hydrate growth). In this case, the growth rate of the hydrate from solution with C = 0.8 is higher than at 0.6, and much higher than at 0.4.

The distribution of average N_H-bond/N_mol in dependence on methane concentration makes (Fig. 7b) it possible to estimate the methane concentration boundaries for successful hydrate formation. At C > 0.2 and C < 0.8 a stable increase of N_H-bond/N_mol number is observed. Thus, this range could be used for hydrate formation at low pressure.

The F₃ and F₄ order parameters were calculated for a 0-100 ns time interval (Figs. 7c and 7d, respectively) for the considered solutions. Solutions with C = 0.4 and C = 0.6 show an increase in the gradual ordering of the H-bond network during the entire simulation time. The difference from the ideal value is due to the fact that a significant fraction of water remained in the liquid phase. The solution with C = 0.8 revealed ordering until t = 23 ns; however the number of hydrate cavities increased. The discrepancy between these results is explained by an increase in methane nano-bubbles and the influence of the bubble surface. The disordering of H-bond network solution with C = 1.0 caused by methane bubble formation is expressed as an increase in the F₃ parameter and a decrease in the F₄ parameter. The solution with C = 0.2 shows no significant change over time. The value F₄ ≈ –0.03 for concentrations of 0.2 and 1.0 corresponds to a liquid (bulk water), and the small difference from F_{4_liquid} is not significant and can be caused by the features of the geometry of the chosen water model.

The obtained results indicate the cooperative nature of the hydrate nucleation process and it is determined by the rearrangement of the hydrogen bond network, which is associated with an individual change in the tetrahedral order parameter of each water molecule in solution.

A change in the distribution of gas molecules also indicates a partial transition of the solution to the hydrate phase. Due to the initial uniform mixing of gas molecules, the peaks in the regions of 3-5 Å (first peak) and 6-8 Å (second peak) are present (Fig. 8a-e) at t = 1 ns. The comparison with bulk hydrate and the increase in the first and second peaks of g(r) with time as well as the appearance of the third peak clearly confirm hydrate growth in solutions with C = 0.4 (Fig. 8b), 0.6 (Fig. 8c) and 0.8 (Fig. 8d). The value of g(r) for the second peak exceeds the height of the first peak (0.4 and 0.6) or reaches a close value (0.8); however, in the case of 0.4 and 0.6, the distribution becomes constant with distance, indicating the presence of a significant amount of methane mixed with water even after t = 100 ns. However, in a system with a methane concentration of 0.8, due to the formation of a bubble (phase separation), the distribution value does not reach a constant value with increasing gas-gas distance. Moreover, at t = 20 ns, the growth of the second peak is insignificant. This is in agreement with theoretical works for carbon dioxide [79] and methane [61] hydrate formation.

The second peak of g(r) for the C = 1.0 solution (Fig. 8e) smooths over 20 ns, which is faster than the disordering of order parameters taking place during 30-40 ns indicates the separation of phases is nearly complete, but residual methane segregates slowly. A decrease of g(r) with r clearly confirms the phase separation.

For the solution with C = 0.2 (Fig. 8a), where no hydrate growth or phase separation is observed over 100 ns, no qualitative differences appear over time. The constant value of g(r) at high r values indicates that the methane concentration is not high enough to break the H-bond network.

The collective nature of the methane hydrate nucleation process in the H-bond network of TIP4P/Ice water outside the region of thermodynamic stability was considered. At T = 270 K and P = 1 bar the process of hydrate formation depended on the methane concentration in water. At methane concentrations of 0.4 (6.5 mol%) and 0.6 (9.45 mol%), growth of the hydrate phase was observed, and at a concentration of 0.8 (12.2 mol%) the competing process of the formation of methane nanophases was found. It can be assumed that the nature of this process is connected with a sufficiently high entropy contribution of the guest molecules to the chemical potential of the water molecules in the hydrate phase. This contribution allows the transformation of the H-bond network of metastable liquid solution into the H-bond network of the stable hydrate lattice.

The collective character of H-bond network reorganization in methane solution was revealed as an increase in water short-range order tetrahedrality at the initial stage of simulation (0–1 ns) before the appearance of the first fluctuation of sI hydrate cavities. Molecules with a high degree of tetrahedrality (low F₃ order parameter values) are distributed uniformly. A significant change in F₃ preceded the formation of cavities, and thus the entire network of hydrogen bonds was involved in the hydrate formation process. Fluctuating cavities appeared in the interval 1.5–6 ns. Later, the fluctuating and stable cavities were constantly present in the system, and an increase in their number indicated the initial stage of the formation of nanohydrate phases.

The formation of hydrate cavities and hydrate nanostructures under the described conditions is possible using real-time Raman spectrometry as in ref [80].

Acknowledgements

The reported study was funded by the Russian Foundation for Basic Research (project number 21-52-52001) jointly with the Ministry of Science and Technology of Taiwan (grant no. 110-2923-E-002-008-MY3). We are grateful to the team at the Center of the Computational Materials Science, Institute for Materials Research, Tohoku University (project no. 2012SC0504) for permanent support of their supercomputer facilities. The authors are also thankful to Prof. T.M. Briere for carefully reading this manuscript.

Data availability

The datasets generated and analyzed during the current study are available from the corresponding author by request.

Author contributions

V.R.B. supervised the project, formulated the tasks, designed the research and wrote the paper. K.V.G. and R.K.Z. performed the simulations and developed the software used for analysis. Y.Y.B., R.V.B., L.-J.C., Y.K. performed the data analysis. All the authors commented on and revised the manuscript.

Competing interests

The authors declare no competing interests.

Sloan, E. D., Jr. & Koh, C. A. Clathrate Hydrates of Natural Gases (CRC Press, 2007).
Khurana, M., Yin, Z. & Linga, P. A Review of Clathrate Hydrate Nucleation. ACS Sustain. Chem. Eng. 5, 11176−11203 (2017).
Bosio, L., Teixeira, J. & Stanley, H. E. Enhanced density fuctuations in supercooled H2O, D2O, and ethanol-water solutions: evidence from small-angle x-ray scattering. Phys. Rev. Lett. 46, 597–600 (1981).
Geiger, A., Mausbach, P., Schnitker, J., Blumberg, R. L. & Stanley, H. E. Structure and dynamics of the hydrogen bond network in water by computer simulations. J. Phys. Colloq. 45, C7-13–C7-30 (1984).
Belosludov V. et al. The nano-structural inhomogeneity of dynamic hydrogen bond network of tip4p/2005 water. Sci. Rep. 10, 7323; 10.1038/s41598-020-64210-1 (2020).
Sosso, G. C. et al. Crystal Nucleation in Liquids: Open Questions and Future Challenges in Molecular Dynamics Simulations. Chem. Rev. 116, 7078−7116 (2016).
English, N. J. & MacElroy, J. M. D. Perspectives on molecular simulation of clathrate hydrates: progress, prospects and challenges. Chem. Eng. Sci. 121, 133–156 (2015).
Owicki, J. C. & Scheraga, H. A. Monte Carlo calculation in isothermal-isobaric ensemble.2. Dilute aqueous solution of methane. J. Am. Chem. Soc. 99, 7413-7418 (1977).
Swaminathan, S., Harrison, S. W. & Beveridge, D. L. Monte Carlo studies on the structure of dilute aqueous solution of methane. J. Am. Chem. Soc. 100, 5705-5712 (1978).
Geiger, A., Rahman, A. & Stillinger, F. H. Molecular dynamics study of the hydration of Lennard-Jones solutes. J. Chem. Phys. 70, 263-276 (1979).
Guillot, B. & Guissani, Y. A computer simulation study of the temperature dependence of the hydrophobic hydration. J. Chem. Phys. 99, 8075-8094 (1993).
Chau, P. L. & Mancera, R. L. Computer simulation of structural effect of pressure on the hydrophobic hydration of methane. Mol. Phys. 96, 109-122 (1999).
Long, J. & Sloan, E. D. Quantized water cluster around apolar molecules. Mol. Simul. 11, 145-161 (1993).
Skipper, N. T. Computer simulation of methane-water solutions. Evidence for a temperature-dependent hydrophobic attraction. Chem. Phys. Lett. 207, 424-429 (1993).
Mancera, R. L., Buckingham, A. D. & Skipper, N. T. The aggregation of methane in aqueous solution. J. Chem. Soc., Faraday. Trans. 93, 2263-2267 (1997).
Silverstein, K. A. T., Haymet, A. D. J. & Dill, K. A. Molecular model of hydrophobic solvation. J. Chem. Phys. 111, 8000-8009 (1999).
Silverstein, K. A. T., Haymet, A. D. J. & Dill, K. A. A simple model of water and the hydrophobic effect. J. Am. Chem. Soc. 120, 3166-3175 (1998).
Mancera, R. L. Computer simulation of effect of salt on the hydrophobic effect. J. Chem. Soc., Faraday Trans. 94, 3549-3559 (1998).
Lazaridis, T. Solvent reorganization energy and entropy in hydrophobic hydration. J. Phys. Chem. B. 104, 4964-4979 (2000).
Lederhos, J. P., Long, J. P., Sum, A., Christiansen, R. L. & Sloan, E. D. Effective kinetic inhibitors for natural gas hydrates. Chem. Eng. Sci. 51, 1221-1229 (1996).
Arjun, Berendsen, T. A., & Bolhuis, P. G. Unbiased atomistic insight in the competing nucleation mechanisms of methane hydrates. PNAS 116, 19305-19310 (2019).
Meng, E. C. & Kollman, P. A. Molecular dynamics studies of the properties of water around simple organic solutes. J. Phys. Chem. 100, 11460-11470 (1996).
Hernandez-Cobos, J., Mackie, A. D. & Vega, L. F. The hydrophobic hydration of methane as a function of temperature from histogram reweighting Monte Carlo simulations. J. Chem. Phys. 114, 7527-7535 (2001).
Broadbent, R. D. & Neilson, G. W. The interatomic structure of argon in water. J. Chem. Phys. 100, 7543-7547 (1994).
De Jong, P. H. K., Wilson, J. E., Neilson, G. W. & Buckingham, A. D. Hydrophobic hydration of methane. Mol. Phys. 91, 99-104 (1997).
Koh, C. A., Wisbey, R. P., Wu, X., Westacott, R. E. & Scoper, A. K. Water ordering around methane during hydrate formation. J. Chem. Phys. 113, 6390-6397 (2000).
Filipponi, A., Bowron, D. T., Lobban, C. & Finney, J. L. Structural determination of hydrophobic hydration shell of Kr. Phys. Rev. Lett. 79, 1293-1296 (1997).
Bowron, D. T., Filipponi, A., Roberts, M. A. & Finney, J. L. Hydrophobic hydration and the formation of a clathrate hydrate. Phys. Rev. Lett. 81, 4164-4167 (1998).
Gao, S., Chapman, W. G. & House, W. NMR and Viscosity Investigation of Clathrate Hydrate Formation and Dissociation. Ind. Eng. Chem. Res. 44, 7373-7379 (2005).
Sloan, E. D. & Fleyfel, F. A. A molecular mechanism for gas hydrate nucleation from ice. AIChE J. 37, 1281-1292 (1991).
Radhakrishnan, R. & Trout, B. L. A new approach for studying nucleation phenomena using molecular simulations: Application to CO₂ hydrate clathrates. J. Chem. Phys. 117, 1786-1796 (2002).
Moon, C., Taylor, P. C. & Rodger, P. M. Molecular Dynamics Study of Gas Hydrate Formation. J. Am. Chem. Soc. 125, 4706-4707 (2003).
Hawtin, R. W., Quigley, D. & Rodger, P. M. Gas hydrate nucleation and cage formation at a water/methane interface. Phys. Chem. Chem. Phys. 10, 4853-4864 (2008).
Jacobson, L. C., Hujo, W. & Molinero, V. Amorphous Precursors in the Nucleation of Clathrate Hydrates. J. Am. Chem. Soc. 132, 11806-11811 (2010).
Vatamanu, J. & Kusalik, P. G. Observation of two-step nucleation in methane hydrates. Phys. Chem. Chem. Phys. 12, 15065-15072 (2010).
Guo, G.-J., Zhang, Y.-G., Liu, C.-J. & Li, K.-H. Using the face-saturated incomplete cage analysis to quantify the cage compositions and cage linking structures of amorphous phase hydrates. Phys. Chem. Chem. Phys. 13, 12048-12057 (2011).
Walsh, M. R., Koh, C. A., Sloan, E. D., Sum, A. K. & Wu, D. T. Microsecond simulations of spontaneous methane hydrate nucleation and growth. Science 326, 1095-1098 (2009).
Guo, G.-J. & Rodger, P. M. Solubility of Aqueous Methane under Metastable Conditions: Implications for Gas Hydrate Nucleation. J. Phys. Chem. B 117, 6498-6504 (2013).
Rodger, P. Methane hydrate: melting and memory. Ann. N. Y. Acad. Sci. 912, 474-482 (2000).
Hallbrucker, A., & Mayer, E. Unexpectedly stable nitrogen, oxygen, carbon monoxide and argon clathrate hydrates from vapour-deposited amorphous solid water: an X-ray and two-step differential scanning calorimetry study. J. Chem. Soc. Faraday Trans. 86, 3785-3792 (1990).
Hallbrucker, A., & Mayer, E. Unexpectedly stable clathrate hydrates formed from microporous vapor-deposited amorphous solid water at low “external” guest pressures and their astrophysical implications. Icarus 90, 176-180 (1991).
Guo, G.-J., Zhang, Y.-G. & Liu, H. Effect of Methane Adsorption on the Lifetime of a Dodecahedral Water Cluster Immersed in Liquid Water: A Molecular Dynamics Study on the Hydrate Nucleation Mechanisms. J. Phys. Chem. C 111, 2595-2606 (2007).
Guo, G.-J., Li, M., Zhang, Y.-G. & Wu, C.-H. Why can water cages adsorb aqueous methane? A potential of mean force calculation on hydrate nucleation mechanisms. Phys. Chem. Chem. Phys. 11, 10427-10437 (2009).
Baez, L. A. & Clancy, P. Computer Simulation of the Crystal Growth and Dissolution of Natural Gas Hydrates. Ann. N. Y. Acad. Sci. 715, 177-186 (1994).
Rodger, P. M., Forester, T. R. & Smith, W. Simulations of the methane hydrate/methane gas interface near hydrate forming conditions conditions. Fluid Ph. Equilibria, 116, 326-332 (1996).
Yagasaki, T., Matsumoto, M., Andoh, Y., Okazaki, S. & Tanaka, H. Effect of bubble formation on the dissociation of methane hydrate in water: a molecular dynamics study. J. Phys. Chem. B 118, 1900-1906 (2014).
Bagherzadeh, S. A., Alavi, S., Ripmeester, J. A. & Englezos, P. Formation of methane nano-bubbles during hydrate decomposition and their effect on hydrate growth. J. Chem. Phys. 142, 214701; 10.1063/1.4920971 (2015).
Knott, B. C., LaRue, J. L., Wodtke, A. M., Doherty, M. F. & Peters, B. Communication: bubbles, crystals, and laser-induced nucleation. J. Chem. Phys. 134, 171102; 10.1063/1.3582897 (2011).
Matsumoto, M., Wada, Y., Oonaka, A. & Onoe, K. Polymorph control of glycine by antisolvent crystallization using nitrogen minute-bubbles. J. Cryst. Growth 373, 73-77 (2013).
Zhang, Y. et al. Effect of Nanobubble Evolution on Hydrate Process: A Review. J. Therm. Sci. 28, 948-961 (2019).
Uchida, T., Yamazaki, K. & Gohara, K. Generation of micro-and nano-bubbles in water by dissociation of gas hydrates. Korean J. Chem. Eng. 33, 1749-1755 (2016).
Uchida, T., Yamazaki, K. & Gohara, K. Gas nanobubbles as nucleation acceleration in the gas-hydrate memory effect. J. Phys. Chem. C 120, 26620-26629 (2016).
Ripmeester, J. A. & Alavi, S. Some current challenges in clathrate hydrate science: Nucleation, decomposition and the memory effect. Curr. Opin. Solid State Mater. Sci. 20, 344-351 (2016).
Walsh, M. R. et al. Methane Hydrate Nucleation Rates from Molecular Dynamics Simulations: Effects of Aqueous Methane Concentration, Interfacial Curvature, and System Size. J. Phys. Chem. C 115, 21241-21248 (2011).
Walsh, M. R. et al. The Cages, Dynamics, and Structuring of Incipient Methane Clathrate Hydrates. Phys. Chem. Chem. Phys. 13, 19951-19959 (2011).
Liang, S. & Kusalik, P. G. Crystal growth simulations of H₂S hydrate. J. Phys. Chem. B 114, 9563-9571 (2010).
Liang, S. & Kusalik, P. G. Exploring nucleation of H₂S hydrates. Chem. Sci. 2, 1286-1292 (2011).
Sarupria, S. & Debenedetti, P. G. Homogeneous Nucleation of Methane Hydrate in Microsecond Molecular Dynamics Simulations. J. Phys. Chem. Lett. 3, 2942-2947 (2012).
Jimenez-Angeles, F. & Firoozabadi, A. Nucleation of methane hydrates at moderate subcooling by molecular dynamics simulations. J. Phys. Chem. C 118, 11310-11318 (2014).
Choudhary, N., Kushwaha, O. S., Bhattacharjee, G., Chakrabarty, S. & Kumar, R. Macro and Molecular Level Insights on Gas Hydrate Growth in the Presence of Hofmeister Salts. Ind. Eng. Chem. Res. 59, 20591-20600 (2020).
Lauricella, M., Meloni, S., English, N. J., Peters, B. & Ciccotti, G. Methane clathrate hydrate nucleation mechanism by advanced molecular simulations. J. Phys. Chem. C 118, 22847-22857 (2014).
Jiménez-Ángeles F. & Firoozabadi, A. Nucleation pathways of gas hydrates from molecular dynamics simulations. Proceedings of the 8th international conference on gas hydrates (ICGH8), 10.13140/2.1.1849.6009 Beijing, China, 28 July - 1 August (2014),
Barnes, B. C. & Sum, A. K. Advances in molecular simulations of clathrate hydrates. Curr. Opin. Chem. Eng. 2, 184-190 (2013).
Wu, J. Y., Chen, L. J., Chen, Y. P. & Lin, S. T. Molecular dynamics study on the growth mechanism of methane plus tetrahydrofuran mixed hydrates. J. Phys. Chem. C 119, 19883-19890 (2015).
Wu, J. Y., Chen, L. J., Chen, Y. P. & Lin, S. T. Molecular dynamics study on the nucleation of methane+ tetrahydrofuran mixed guest hydrate. Phys. Chem. Chem. Phys. 18, 9935-9947 (2016).
Gets, K. et al. Transformation of hydrogen bond network during CO₂ clathrate hydrate dissociation. Appl. Surf. Sci. 499, 143644; 10.1016/j.apsusc.2019.143644 (2020).
Belosludov, V. R. et al. Collective Effect of Transformation of a Hydrogen Bond Network at the Initial State of Growth of Methane Hydrate. JETP Lett. 115, 124-129 (2022).
Plimpton, S. Fast parallel algorithms for short-range molecular dynamics. J. Comput. Phys. 117, 1-19 (1995).
Nosé, S. A molecular dynamics method for simulations in the canonical ensemble. Mol. Phys. 52, 255-268 (1984)
Hoover, W. G. Canonical dynamics: equilibrium phase-space distributions. Phys. Rev. A 31, 1695-1697 (1985).
Hockney, R. W. & Eastwood, J. W. Computer simulation using particles. (McGraw-Hill, 1981)
Ryckaert, J. P., Ciccotti, G. & Berendsen, H. J. Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J. Comput. Phys. 23, 327-341 (1977).
Abascal, J. L. F., Sanz, E., Fernández, R. G. & Vega, C. A potential model for the study of ices and amorphous water: TIP4P/Ice. J. Chem. Phys. 122, 234511; 10.1063/1.1931662 (2005).
Martin, M. G. & Siepmann, J. I. Transferable Potentials for Phase Equilibria. 1. United-atom Description of n-Alkanes. J. Phys. Chem. B 102, 2569-2577 (1998).
Parui, S., & Jana, B. Factors promoting the formation of clathrate-like ordering of water in biomolecular structure at ambient temperature and pressure. J. Phys. Chem. B 123, 811-824 (2019).
Uchida, T. et al. Contribution of ultra-fine bubbles to promoting effect on propane hydrate formation. Front. Chem. 8, 480; 10.3389/fchem.2020.00480 (2020).
Belosludov, V. R. et al. Physical and chemical properties of gas hydrates: Theoretical aspects of energy storage application. Mater. Trans. 48, 704-710 (2007).
Makogon, T. Y. & Sloan, E. D. J. Phase equilibrium for methane hydrate from 190 to 262K. J. Chem. Eng. Data 39, 351-353 (1994).
Hirai, S., Okazaki, K., Tabe, Y. & Kawamura, K. CO₂ clathrate-hydrate formation and its mechanism by molecular dynamics simulation. Energy Convers. Manag. 38, S301-S306 (1997).
Kang, K. C., Hong, S. Y., Cho, S. J., Kim, D. H. & Lee, J. D. Evaluation of desalination by nanostructured hydrate formation and pellet production process. J. Nanosci. Nanotechnol. 17, 4059-4062 (2017).

No competing interests reported.

Download PDF

Version 1

posted

You are reading this latest preprint version

Molecular dynamics study of unexpected hydrate formation from supersaturated methane solution at low pressure

Status:

Version 1

Abstract

Figures

Introduction

Calculation Methods

Results And Discussion

Conclusions

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1