Collective Effect of Transformation of a Hydrogen Bond Network at the Initial State of Growth of Methane Hydrate

The molecular dynamics study of the rearrangement of the dynamic hydrogen bond network of liquid water to the crystal hydrogen bond network of methane hydrate in the process of their formation and growth is conducted. To reveal the nature of nucleation, the time dependences of the degree of “crystallinity” of the nearest environment of all water molecules, the rate of ordering of the hydrogen bond network, and the relation of these parameters to the rate of growth of methane hydrate are studied. The effect of the presence of sea salt ions and hydrate seed on these parameters is analyzed. Systems with a completely mixed gas, i.e., with the minimum induction time, are fabricated, and it is shown that gas hydrates can be nucleated simultaneously in the entire volume of the solution, which in turn indicates the collective formation of hydrates from liquid solution.


INTRODUCTION
The problem of access to pure water has become more relevant in many regions of the Earth in the last several decades. The only modern industrial solution to this problem is seawater desalination using various methods among which the hydrate method is outstanding, which is very close to the classical freezing desalination of seawater [1]. This method is based on the fact that salt ions remain in the aqueous solution when the crystal phase is formed. However, the problem of separation of the hydrate phase and seawater after the formation of hydrate and the problem of slow kinetics of the formation are still incompletely solved [2].
Gas hydrates are inclusion compounds, which are formed by displacing guest molecules in the cavity of the crystalline skeleton formed by water molecules [3]. The type of formed hydrate depends on the gas or the mixture of gases involved in formation. The most widespread structures of hydrates are cubic structure I (sI), cubic structure II (sII), and hexagonal structure III (sH).
Salt ions Na + , Cl -, Mg 2+ , Ca 2+ , and K + present in seawater are usually used as inhibitors of formation of hydrates and, consequently, affect the kinetics of their formation. To develop simple methods of syn-thesis of hydrates, it is of key importance to understand hydrate-formation mechanisms, particularly in the presence of various inhibiting impurities. The experimental study of kinetics is very difficult because nano-or even femtosecond resolution is required to directly observe the formation of hydrates [3]. At the same time, the molecular dynamics method is the most appropriate method for the study of these processes, although it requires large computational resources [4].
To understand the formation of hydrates and the dynamics of their growth, it is necessary to reveal the mechanism of nucleation, which is still incompletely clear because the process is stochastic. At present, there are two main theories of nucleation: homogeneous and heterogeneous nucleation. In the case of homogeneous nucleation, the new phase is formed directly from the initial phase after the formation of the so-called critical nucleus; after that, the growth of the hydrate phase occurs. It is noteworthy that this process is stochastic. In real systems, it is very difficult to observe heterogeneous nucleation because of the presence of impurities and numerous secondary factors, which are responsible for heterogeneous nucleation [5]. Heterogeneous nucleation occurs in the presence of the third phase, e.g., impurity or phase interface. The heterogeneous way for the formation of

CONDENSED MATTER
hydrate is energetically more favorable than the formation from the pure phase, which significantly increases the rate of nucleation [6]. Unfortunately, this mechanism is poorly studied because it is complex [7]. Studies of nucleation of hydrates continue because any commonly accepted opinion in favor of one of these two theories is absent [8][9][10][11]. During last decades, numerous qualitive and quantitative molecular dynamics studies the initial stages of the formation of hydrate at high excess pressure of low temperature upon supercooling were performed [12,13]. A key component of these studies is the calculation of the order parameters [14,15], which are a quantitative measure of the degree of ordering in the system and are used to estimate the dynamics and mechanisms of formation of hydrates and other key properties. The processes of formation of hydrates are associated with the ordering of both guest molecules and water molecules [16], which form the hydrogen bond network. Its role is important to understand the formation of hydrates because the hydrogen bond network determines the local properties of water, e.g., structural inhomogeneities [17], macroscopic properties of water, e.g., viscosity [18], thermodynamic properties of hydrates [19 and reference therein], and proton conductivity [20,21].
In this work, we implement the heterogeneous way of nucleation. The main aim is to study the transition of the hydrogen bond network of liquid water to the stable hydrogen bond network of hydrate. To this end, we use the molecular dynamics method and a set of our programs for the postprocessing and analysis of the resulting data. The study is performed for the water + methane, water + methane + hydrate seed, water + methane + sea salt ions, and water + methane + sea salt ions + hydrate seed model systems.

DETAILS OF CALCULATIONS
In this work, we combine the molecular dynamics method implemented in the LAMMPS package [22] and a set of our programs for the structural analysis of the hydrogen bond network of water molecules.
The simulation was performed for several systems of various sizes and compositions. The first group of systems consisted of 1200 water molecules and 104 methane molecules; this amount of methane makes it possible to transfer half of the water molecules into hydrate. We produced variants of this system with seawater salt ions Na + and Cl -(the other ions in seawater were disregarded in the model because their concentrations were very low) with a concentration of one pair of ions per 100 water molecules, as well as systems consisting of the large and small cavities of the sI structure joined together, and without hydrate seed.
Water molecules were described by the four-point TIP4P/Ice potential [23], which is appropriate to describe phase transitions between different phases of water. Methane molecules were treated as spherical neutral particles with van der Waals potential [24]. The parameters of the potential for Na + and Clwere taken from [25,26].
The initial structures were created by a random distribution of the required numbers of corresponding molecules; after that, the structure was optimized in order to eliminate overlapping of molecules. In the system with hydrate seed, seed was considered as a solid; i.e., the relative position of molecules was fixed. The temperature and pressure in the simulation were controlled by a thermostat and a Nosé-Hoover barostat, respectively [27,28]. Systems were simulated at a temperature of 270 K, pressures of 50 and 1000 bar, and a step of 1 fs.
The molecular structure of solutions, the coordinates of atoms and molecules, and the mutual orientation of water molecules dynamically vary in time. In this work, the hydrogen bond network is considered as dynamic [17]; i.e., it joins fluctuating hydrogen bonds between neighboring water molecules, whose formation, existence, and decay at each time are determined by the geometric criteria: O-O distance Å and mutual orientation angle ([29], the number of bonds with is insignificant). These criteria were used to analyze instantaneous configurations and to search for small and large cavities of sI hydrate. The number of cavities and average values of the characteristics of the hydrogen bond network described below were calculated for individual molecular configurations obtained every 100 ps. It is important that only bonds formed between water molecules are taken into account when analyzing hydrogen bond network.
The parameters F 3 [14] and F 4 [15] were used to analyze the structure of water. The parameter F 3 is calculated as the value average over all water molecules , where α is the O-O-O angle between molecules nearest to the considered water molecule in the center and summation is performed over six independent angles that make the molecule with four neighboring molecules when the ice rule is valid. This parameter characterizes the coefficient of difference of the local structure of the water molecule from the ideal tetrahedral structure and is 0 and ≈0.1 for the perfect crystal and liquid water, respectively. The ideal tetrahedral angle is 109.47°. The parameter F 4 = , where ϕ is the torsion angle between nearest molecules, is the average value of the torsion angles between all neighboring water molecules. The characteristic F 4 values for hydrates with liquid water and ice are approximately 0.7 and -0.3, respectively.
RESULTS AND DISCUSSION Figure 1 shows the time dependences of the numbers of (a) small 5 12 and (b) large 5 12 6 2 cavities filled with sI hydrate methane molecules formed in the system. In the systems without hydrate seeds, small cavities are formed earlier than large ones because the former are energetically favorable than the latter. Large cavities are formed only in the immediate vicinity of small cavities because faces of small cavities are involved in their formation. The hydrate seed promotes the formation of large cavities. The presence of salt ions insignificantly affects the rate of growth or the number of cavities. The ratio of the numbers of small and large cavities in the sI hydrate is from 2 to 6, but this ratio in our case is different because of a small number of methane molecules in the system. It is seen that the number of cavities in all systems is saturated in 150-200 ns.
The method chosen to form initial structures allowed one to reduce the induction time, which reaches several microseconds in other molecular simulation studies, to insignificant values [10]. In this work, a configuration with the content of methane in the solution sufficient for the beginning of the growth of hydrate is formed, which is a necessary condition of hydrate formation. According to the preliminary data for Freon 14, propane, or isobutane, this effect can also be manifested in other hydrate systems.
The convenient parameters F 3 and F 4 were used to analyze the spatial orientation of water molecules and their distribution. The time dependence of the parameter F 3 is shown in Fig. 2a, where it is seen that the structures of solutions are ordered with time, according to a decrease in the parameter F 3 with time. Ordering in the tetrahedral surrounding of water molecules at 50 bar occurs more rapidly than at 1000 bar, which is due to a higher mobility of water molecules at a lower pressure. Similarly, it is seen that the presence of salt in water reduces the nearest tetrahedral surrounding of water molecules in the hydrogen bond network, which is confirmed by X-ray scattering data and simulation in [30], where a decrease in the number of hydrogen bonds in the short-range order in the Na + + Clsolution compared to pure water was observed. The dependences approach a constant value and oscillate around this value because almost all methane molecules initially present in the solution pass to the hydrate phase.
The time dependence of the parameter F 4 is shown in Fig. 2b, where it is seen that the presence of hydrate seed leads to a more pronounced intermolecular ordering independent of the pressure, which indicates the gradual formation of hydrate. The comparison of data for the water + methane and water + methane + sea salt ions systems without seed shows that the presence of salt makes the mutual orientation of water molecules more inherent in hydrate, which is explained by a decrease in the mobility of water molecules interacting with Clions [30]. Similar to the parameter F 3 , hydrate grows, involving almost all methane molecules initially existing in the solution. The resulting value F 4 ≈ 0.2-0.25 differs from a value of 0.7 for hydrate because only half of water molecules can be involved in hydrate and water molecules that are on the surface of hydrate and form hydrate make a significant contribution.
It is seen that the parameters F 3 and F 4 are saturated in ~100 ns, which significantly differs from a similar time for the dependences of the numbers of small and large cavities in hydrate. This indicates that the ordering of the hydrogen bond network precedes a visible growth of hydrate itself. In other words, the growth of hydrate involves a larger volume of the system than the volume occupied by formed cavities. (Color online) Numbers of (а) small 5 12 and (b) large 5 12 6 2 cavities containing methane versus the time of simulation. Data for the water + methane, water + methane + hydrate seed, water + methane + sea salt ions, and water + methane + sea salt ions + hydrate seed systems are shown in black, orange, blue, and green, respectively. The solid and dotted lines show data at pressures of 50 and 1000 bar, respectively. Figure 3 shows the number of molecules with the parameter F 3 < 0.025 (which are the closest to the crystal phase in water phases) divided by the total number of water molecules, as well as the number of hydrogen bonds connecting these molecules. An increase in the number of such molecules indicates the general ordering of the hydrogen bond network. It is seen that the number of tetrahedrally ordered molecules approaches a constant value in the first 100 ns, as for the parameters F 3 and F 4 . To the same time, the strongest binding of these molecules by hydrogen bonds is observed; the fraction of these molecules reaches 40-45 and ~35% in systems with and without seeds, respectively. The degree of binding in the first 50 ns is rather low, particularly, in systems without seed, which indicates the delocalization of such molecules and the absence of a single nucleation center. A high degree of binding shows that most of these molecules form a solid phase. Thus, it is seen that hydrate continues to be formed when the hydrogen bond network is transformed and corresponds to the solid phase. In this case, an "amorphous" hydrate structure can be formed, transferring in time to the crystal phase.

N N
The bulk transition of the hydrogen bond network of water to the hydrate network in the case of the uniform initial distribution of the gas can be observed as the formation of small and large cavities in the hydrate structure over the entire model volume, rather than only in a certain region. Figure 4 shows the time dependence of the ratio of the numbers of molecules forming small and large cavities to the total number of cavities. It is seen that cavities at small times are not connected to each other (have a smaller number of   Fig. 3. (Color online) Time dependences of (a) the fraction of molecules with F 3 < 0.025 of the total number of molecules and (b) the number of hydrogen bonds coupling these molecules divided by the number of molecules with F 3 < 0.025. Data for the water + methane, water + methane + hydrate seed, water + methane + sea salt ions, and water + methane + sea salt ions + hydrate seed systems are shown in black, orange, blue, and green, respectively. The solid and dotted lines show data at pressures of 50 and 1000 bar, respectively.
No. 3 2022 common molecules than those in a single pentagonal or hexagonal face). With the time, they begin to be connected and fill a large fraction of the volume of the model cell, which reduces the ratio, which is equal to 5.75 in the case of ideal hydrate. In the absence of the seed, larger 5 12 6 2 cavities consisting of 24 water molecules are more often formed at initial times, whereas the formation of small cavities consisting of 20 water molecules occurs at a higher rate around the seed. CONCLUSIONS To summarize, the water + methane and water + methane + sea salt ions systems at various pressures have been studied by the molecular dynamics method. It has been shown that the creation of gas-saturated solution can significantly reduce the induction time and provide a rapid growth of the hydrate phase from the prepared solution. Growth has occurred over the entire volume simultaneously without the localization of nucleation centers. Furthermore, it has been shown that the rearrangement of the hydrogen bond network begins and ends earlier than the alignment of water molecules into the crystal structure of hydrate is visually observed. The presence of seed promotes a higher rate of hydrate formation even under conditions of the complete mixing of methane in water. The presence of seawater ions insignificantly affects the kinetics of hydrate formation but affects the geometry of the short-range order of water molecules. Data for the water + methane, water + methane + hydrate seed, water + methane + sea salt ions, and water + methane + sea salt ions + hydrate seed systems are shown in black, orange, blue, and green, respectively. The solid and dotted lines show data at pressures of 50 and 1000 bar, respectively. N N