As a dynamic research method for molecular systems, molecular dynamic (MD) simulation can represent physical phenomena that cannot be realized by experimental means and discuss the microscopic reaction mechanism of things from the molecular level. In this paper, the previous research results were reviewed. First, the MD simulation process was briefly described, then, the applicability of different molecular force fields in the natural gas hydrate (NGH) system was discussed, and finally, the application of MD simulation in the formation and decomposition law of NGH was summarized from the perspective of NGH mining. The results show that the selection of water molecular force field has a great influence on the simulation results, and the evaluation of water model applicable to the simulation of NGH under different thermodynamic states is still an open research field that needs to be paid attention to. The effect of surface properties of porous media (such as crystallinity and hydrophilicity) on hydrate needs to be further studied. Compared with thermodynamic inhibitors, kinetic inhibitors (such as amino acids) have more promising research prospects, and further research can be carried out in the screening of efficient kinetic inhibitors in the future.

1. Introduction

NGH is an ice-like solid substance formed by the combination of natural gas (mainly methane) with water under high pressure and low temperature, generally distributed in deep-sea sediments or permafrost [1, 2]. As a new energy source in the future, NGH has the advantages of large reserves, wide distribution, low pollution, and high energy density. Kvenvolden (1988) estimated that the global reserves of organic carbon in hydrate were , which was twice the total amount of coal, oil, and natural gas resources discovered on the earth so far [3]. In the standard state, 1 m3 NGH can be decomposed to produce 160~180 m3 of natural gas [1]. According to the size of guest molecules and the type and number of water molecules’ cavities, hydrate crystal structures can be divided into three categories, namely, type I [4], type II [5], and type H [6]. The structural differences of different types of hydrate are shown in Table 1 and Figure 1. As can be seen from the table, type I hydrate is a volumetric cubic structure, with a unit cell containing 46 water molecules, forming 2 small cavities and 6 large cavities. Among them, the small cavity is a pentagonal dodecahedron (512), while the large cavity is a tetrahedron (51262) composed of 12 pentagons and 2 hexagons [4]. Type II hydrate is composed of 136 water molecules, including 16 small cavities (512) and 8 large cavities (51264) [5]. Type H hydrate is a hexagonal crystal structure, and its unit cell contains three small cavities (512), two hollow cavities (435663), and one large cavity (51268) composed of 34 water molecules [6]. The most common methane hydrate in nature is generally of type I structure. In all types of hydrates, each cavity can contain only one guest molecule at most. In extreme cases (such as hydrogen hydrates under high pressure), multiple guest molecules may occupy a single cavity together [7].

The extraction methods of NGH mainly include depressuring, thermal stimulation, chemical reagent injection, and gas swapping [9]. Among them, the depressuring method is the most promising hydrate decomposition technology [1013], which mainly promotes hydrate decomposition by reducing the reservoir pressure below the hydrate equilibrium pressure. The depressuring method is economically feasible, but its action is slow. Moreover, in the process of mining, due to the heat absorption caused by hydrate decomposition, the reservoir temperature drops, and the free water will form ice, which hinders the further decomposition of hydrate [14, 15]. Thermal stimulation method refers to heating the reservoir in various ways, so that the reservoir temperature reaches above the equilibrium temperature, thus, promoting the decomposition of hydrate. Common heating methods include two types: one is to directly inject heat sources such as hot water, steam, or heating liquid into the formation; the other is to indirectly heat the formation through sound waves or electric tools [16]. There is a great heat loss in the process of heat injection, and the energy cannot be fully utilized. The method of chemical reagent injection is to inject chemical substances into the reservoir to change the thermodynamic equilibrium conditions for the stable existence of NGH. Generally, the chemical reagents are of high cost and will cause certain pollution to the formation [17]. The gas swapping method is to inject another guest molecule (such as CO2) into the reservoir. Compared with CH4, this guest molecule is easier to combine with water to form hydrate, and then the methane gas in the hydrate is replaced [18]. This method can always maintain the stability of the reservoir and can achieve the burial of CO2 gas, but there is a problem of low replacement efficiency [19]. The above four methods are all to convert solid hydrate into gas and water under certain temperature and pressure conditions, and then to extract gas. However, the hydrate reservoir is dense with extremely low permeability, and secondary hydrate is formed due to the change of temperature. Therefore, with the progress of exploitation, the saturation of each phase, effective porosity, and permeability of the reservoir all change dynamically [20], simple macroscopic physical experiment cannot accurately simulate the hydrate mining process.

As an effective tool for the analysis of the behavior of complex systems, molecular dynamic (MD) simulation technique can explore the microscopic mechanism of various influencing factors in the simulation system from the molecular level. The technology has a certain research prospect in the simulation of hydrate phase behavior [2123], nucleation and growth [2426], transport properties [27, 28], and chemical inhibition [29, 30]. It can realize the physical phenomena and processes that cannot be achieved by the current experimental means, thus, providing a powerful theoretical guidance for the experimental research and thus improving the efficiency of research and development. Figure 2 shows the main flow of MD simulation. It can be seen from the flow chart that basic information of the system, such as initial temperature, number of particles, and type of atoms, is first read in. Then, the position and velocity are assigned initial values to calculate the potential energy of the system and the forces acting on all the particles. Finally, the Newtonian equation of motion is solved iteratively. The calculation of the potential energy and the solution of the equation of motion are repeated until the calculation system reaches a specified number of time steps. After the cycle is completed, the mean value of the measured quantity is calculated, and the simulation ends. In classical mechanics, forces are defined as follows: where is the force exerted on the particle; is the potential energy of the particle; is particle mass; is particle position; is particle velocity; is particle acceleration.

In Equation (1), the forces on the particles include intramolecular interactions (bond interactions) and intermolecular interactions (nonbond interactions). Among them, intramolecular interaction includes molecular bond length stretching term, bond angle bending term, dihedral angle torsion term, and out-of-plane vibration term, and intermolecular interaction includes electrostatic interaction and van der Waals interaction. Molecular potential energy can be simply expressed as follows [16]: where is the molecular potential energy; , , , and are bond length stretching term, bond angle bending term, dihedral angle torsion term, and out-of-plane vibration term, respectively. The and are van der Waals and electrostatic potential energy, respectively.

At present, there are many algorithms to solve the Newtonian equation of motion. The three most commonly used algorithms are the Verlet algorithm [31], leapfrog algorithm [32], and velocity-Verlet algorithm [33]. Where Verlet algorithm is derived from Taylor’s formula. By expanding the Taylor series of coordinates of a particle at time , position and velocity expressions can be obtained, respectively:

Verlet algorithm has time reversibility, but the denominator of its velocity expression contains the time step term, which makes the calculation of velocity have a large error. The leapfrog algorithm defines the particle velocity at half a time step. Its velocity and position are not defined at the same time, so it is unable to calculate the kinetic energy and potential energy of the particle at a certain time simultaneously. The velocity-Verlet algorithm is proposed on the basis of the Verlet algorithm. It has high precision, and the position and velocity are defined synchronously, but the calculation speed is slightly lower than that of the frog hop algorithm.

After the end of the cycle, different thermodynamic and kinetic properties, such as density, viscosity coefficient, mechanical parameters, potential and energy terms, etc., are calculated by statistical mechanics [8]. Existing molecular dynamic computer simulation programs mainly include DL_POLY [34], GROMACS [35], and LAMMPS [36].

Based on the above analysis, this paper will summarize the selection of the molecular field of the NGH system, the nucleation mechanism of NGH, and the application of MD in hydrate mining (such as depressuring, thermal stimulation, gas swapping, and chemical injection), and finally give some suggestions.

2. Force Field Model

2.1. CH4 Molecular Force Field Model

In general, intermolecular forces include van der Waals interactions, Coulombic interactions, and other interactions. CH4 molecule is a simple nonpolar spherically symmetric molecule, for nonpolar molecules, van der Waals force can usually be described by dispersion force, and the commonly used Lennard-Jones (LJ) potential function is described in this way, which is expressed as: where is the distance between atomic pairs; is the potential well depth, which reflects the strength of the interaction between two atoms. is the distance between atoms when the action potential is equal to zero; the 12 times the power represents repulsion, and the 6 times the power represents attraction. In methane hydrate simulation, single point Lennard-Jones (LJ) potential function [37], united-atom LJ potential function [38], and OPLS-AA (potentials for liquid calm-all atom) potential function [39] are commonly used to characterize methane molecular force field.

2.2. Carbon Dioxide Field Model

The CO2 molecule is a simple nonpolar linear molecule. Although it has no polarity itself, it contains C=O polar bond. Therefore, C atom and O atom will have an obvious partial charge, and constraints on the C=O bond length and O=C=O bond angle are required. At this point, the use of a single LJ potential function to describe the CO2 molecular force field will bring some errors. Based on the above analysis, the LJ potential function combined with electrostatic interaction, that is, the common EPM (Elementary Physical Model) series Model [40] can be used to represent the CO2 molecular force field. The expression of the EPM model is as follows:

In the above equation, the first term is the constraint on the bond length, representing the rigid bond length; the second term is the constraint on the bond angle, which represents the resonant bond angle. The last represents electrostatic interaction, which works by introducing three charge points into the three atomic centers of a CO2 molecule. The gas-liquid coexistence curve and critical properties predicted by the EPM model are very close to the experimental values. In addition, the EPM2 model can be obtained by recalculating the potential parameters of the EPM model, which can effectively describe the gas-liquid behavior near the critical point of CO2 [38].

2.3. H2O Force Field Model

Water molecules are polar molecules with obvious partial charge on the atoms. Therefore, the LJ potential function combined with electrostatic interaction is usually used to express the water molecular force field. When using the MD method to simulate gas hydrate, there are a variety of water molecular potential functions to choose from, among which simple point charge (SPC) series and TIPnP series are most frequently used. The SPC potential function assumes that water molecules are rigid molecules. The interaction between molecules includes the LJ potential energy of oxygen atoms between different molecules and the Coulomb potential generated by the charge between oxygen atoms and hydrogen atoms in each molecule. The potential function expression can be written as follows:

The TIPnP series builds on the SPC series by introducing virtual atoms with only charge but no mass to describe hydrogen bonds, where can be 3 [41, 42], 4 [43, 44], and 5 [45, 46]. Because the TIP3P potential function is difficult to accurately characterize the properties of liquid water, and the TIP5P potential function has the disadvantage of complex calculation, the TIP4P potential function is often used in the actual hydrate simulation [47]. The TIP4P series of potential functions includes TIP4P/Ice [43], TIP4P/2005 [44], and TIP4P-EW [48]. However, it is known that the physical properties predicted by different water models tend to produce different results [22]. Li et al. used SPC/E, TIP4P, and other water models to study the microdecomposition characteristics of hydrate under the NVE ensemble. The results showed that the decomposition rate predicted by the SPC/E model was higher than that predicted by the TIP4P model under the same temperature and pressure conditions [49]. Choudhary et al. calculated the melting points of methane hydrate using TIP4P/ICE, TIP4P/2005, TIP4P, and SPC/E water models, respectively. The results showed that the melting points for methane hydrate at 10 MPa calculated by TIP4P/ICE and TIP4P/2005 water models were 325 K and 295 K, respectively, higher than that of the experimental results, while the melting points calculated by TIP4P and SPC/E models were between 270 and 265 K, lower than that of the experimental output [50]. Therefore, it is important to evaluate a more suitable water model under given thermodynamic conditions. In general, different water molecular models are needed to simulate different temperature states, such as the supercooling required for hydrate crystallization and the overheating required for decomposition [5155]. Studies have shown that TIP4P/Ice [22], TIP4P/2005 [56], and TIP4P-EW [57] are generally applicable to describe hydrate formation, while TIP4P and SPC/E are more applicable to hydrate decomposition [53, 58, 59]. However, the predictive nature of these water models still varies. Compared with the TIP4P/Ice model, the TIP4P/Ice model can predict the hydrate formation at higher temperature and lower pressure [50]. Therefore, it is still an open research field to explore the water model applicable to the study of NGH in different thermodynamic states.

3. Regularity of Hydrate Nucleation

Cage hydrate is a solid crystal structure formed by nucleation of mixed solution consisting mainly of water and guest molecules under certain thermodynamic conditions. Understanding the nucleation mechanism of hydrate is of great significance for understanding the basic chemical principle and application of these complex structural compounds. Because the size of the critical nucleus and the nucleation process are all nanoscale, it is impossible to observe directly at the molecular level by experimental method. Therefore, MD simulation is an ideal tool to study the nucleation at the molecular level.

At present, there are mainly two theories about hydrate nucleation. The first theory is called the “unstable cluster aggregation hypothesis” [60], which holds that methane molecules in solution are surrounded by clathrate-shaped dissolving shells that have the characteristics of hydrates, and that hydrates are formed from these clathrate clusters. The second theory, known as the “local structure hypothesis” [61], is concerned with the local ordering of guest molecules: random fluctuations cause a certain number of methane molecules to appear in a hydrate-like arrangement, which in turn induces water molecules to form hydrate-specific hydrogen bond structures. At the basic research level, the difference between the two theories is whether the order of water is driven by guest molecules, or whether the order of guest molecules is driven by water molecules. At the practical level, the unstable clustering method relies on particle-cluster aggregation to grow, while the local structure model requires the guest molecules to be more clustered and move over a longer distance [62].

According to the nucleation theory, it is necessary to overcome the Gibbs free energy barrier to form stable methane hydrate from unstable states [63]. The size of the barrier is determined by the energies required to form crystal-fluid interface and the crystal itself. The driving force to overcome the free energy barrier is the difference in chemical potential between the new and the old phases [64]. Under the conditions of low temperature, high pressure, and supersaturation (excessive methane in solution), the driving force is higher, which is beneficial to hydrate formation. In the theory of nucleation, nucleation time is strongly correlated with supersaturation, and with the increase of supersaturation, the nucleation time decreases gradually.

In order to understand the mechanism of hydrate nucleation at the molecular level, many scholars have tried to use the MD method to study the nucleation and growth of hydrate. The hydrate nucleation simulation was initially carried out in the bulk phase, but this approach used an unreasonably high methane concentration and assumed that the movement of hydrophobic guest molecules in the liquid was unrestricted [55, 65]. To avoid this problem, common hydrate nucleation studies were carried out at the gas-water interface. Moon et al. (2003) simulated the nucleation process of methane hydrate at the methane/water interface using the MD method, and the results showed that under moderate supercooling conditions, the growth rate of hydrate was accelerated, and a hydrate cluster of 280 water molecules could be formed within 10 ns [66]. Zhang et al. (2008) carried out MD simulation in water/methane interface under the temperature of 220, 235, 240, 245, and 250 K, respectively. The results showed that the nucleation of the hydrate has a significant dependence on temperature. The methane hydrate nucleates rapidly at lower temperature, while the hydrate structure disappears completely at higher temperature; 240 K is the transition temperature between the growth and the melting of the hydrate [67]. Hawtin et al. (2008) carried out MD simulation of spontaneous nucleation of methane hydrate at the methane/water interface. The research results showed that the initial nucleation of hydrate was not consistent with any common bulk crystal structure, but containing the structural units of these crystals. In addition, the formation of water cages is closely related to the collective arrangement of methane molecules [62]. Walsh et al. (2009) constructed a methane-water equilibrium interface at 425 K and provided snapshots of methane hydrate nucleation at different time steps, as shown in Figure 3. The results showed that the cages formed in the early stage were all small cages, and some cages shared a surface, which was in accordance with type II structure. However, due to spatial constraints and a thermodynamic preference for structure Type I, larger cages subsequently emerged, nucleating and growing to result in a combination of two main type (type I and type II) hydrate crystals. The two types of hydrate are connected by an uncommon 51263 cage, which promotes the coexistence of the two types of hydrate structures [68]. Zhang et al. (2015) carried out 6 precision NVE MD simulation in a methane-water system under 80 MPa and 250 K. The results showed that methane hydrate nucleation through various channels, including the formation of metastable amorphous phase intermediates, as well as direct forming hydrate crystalline phase, which proved that it is possible to directly form stable crystalline phases during the process of hydrate nucleation [69]. In existing studies, the formation of hydrate structures usually requires very high pressure, dynamic changes in temperature, or preexisting hydrate structures. Felipe and Abbas simulated supersaturated methane-water mixture at a series of temperatures. They found that as the amount of supercooling reduced, the order of the crystal structure increased. Crystal structures I and II form and coexist at moderate temperatures. Crystallization began with spontaneous formation of amorphous clusters in which structure I, II, and other ordered defects occurred [55]. In addition to supercooling, the formation of hydrates is also triggered when the concentration of methane in aqueous solutions exceeds a certain value. Walsh et al. quantified the mole fraction of methane dissolved in water at the moment of nucleation. They found that the mole fraction was strongly correlated with temperature, pressure, and interface geometry. In the induction period before nucleation, the methane mole fraction was 0.02~0.04 when the temperature was 245~250 K and the pressure was in the range of 50~4000 bar [55]. Guo and Rodger found that the critical value of methane mole fraction was around 0.05, beyond which the hydrate structure would spontaneously form [70].

In nature, NGH can form in different geological environments, such as in sedimentary rocks or loose clays. Therefore, the effect of solid surface properties (such as crystallinity and hydrophilicity) on hydrate formation should be considered in the study of hydrate nucleation. The influence of solid surface on hydrate formation may be caused by changing the local structure of water molecules near the liquid-solid interface, thus, changing the hydrate formation pathway. For example, experimental measurements show that methane hydrates form faster under the surface of bentonite than in solution [71]. Existing molecular simulation is mostly focused on the dynamics research of hydrate formation and decomposition in the gas-liquid two-phase system [26, 66, 68, 7274] and limited gas-liquid-solid three-phase system [75, 76]. However, the properties of water near the hydrophilic silica surface are different from those of water in the bulk phase, and the surface of silica can promote the growth of methane hydrate [75, 76]. Bai et al. studied the nucleation process of CO2 hydrate caused by hydroxylation of silicon surface in two-phase [77] and three-phase [78] systems and found that the presence of silicon surface could accelerate the nucleation of hydrate, partly because the surface of silicon stabilized the structure of adsorbed water. On the basis of the above studies, Bai et al. further studied how the hydrophilicity and crystallinity of solid surfaces regulate the local structure of adjacent molecules and the nucleation of CO2 hydrate [79]. The results showed that the hydrophilicity of solid surface changes the local structure of water molecules and the distribution of gas phase near the liquid-solid interface, thus, changing the mechanism and dynamics of gas hydrate nucleation. Hydrate nucleation is more likely to occur on a relatively hydrophilic surface. Different from surface hydrophilicity, surface crystallinity has little influence on the local structure of adjacent water molecules and the nucleation of gas hydrate. During the initial phase of gas hydrate growth, the molecular structure induced by the crystal surface is more orderly than that induced by an amorphous solid surface. The hydrate formation process under different solid surfaces is shown in Figure 4. In the figure, the systems on different silicon surfaces are labeled as crys-n or amor-n, where “crys” and “amor” correspond to crystalline and amorphous solid surfaces, respectively, and the amount of represents the percentage of -OH groups on the surface.

4. Law of Hydrate Decomposition

At present, the kinetics of hydrate decomposition is becoming more and more important. On the one hand, only by fully understanding the decomposition kinetics of hydrate can a feasible hydrate mining technology be developed. On the other hand, a good knowledge of the kinetics of hydrate decomposition is very important for the design of kinetic inhibitors or anticondensate inhibitors. The exploitation methods of gas hydrate mainly include depressuring, thermal stimulation, gas swapping, and chemical reagent injection. The principle of the former two methods is to make the reservoir conditions cross the solid-liquid equilibrium curve, and the principle of the latter two is to make the phase line move to the left and up. The schematic diagram of the four mining methods is shown in Figure 5. Currently, relevant laboratory experiments and field-scale simulations have been carried out, but the decomposition information of methane hydrate at the microscopic level is still limited. It is of great significance to explore the hydrate decomposition mechanism on the molecular scale for improving the feasibility of hydrate large-scale exploitation and solving the hydrate plugging problem in natural gas pipelines.

4.1. Depressuring and Thermal Stimulation

This section briefly introduces the MD simulation process of gas hydrate decomposition caused by depressuring and thermal stimulation. Because it is difficult to reduce the pressure of the hydrate system to realize the hydrate decomposition process, there are few reports on the decomposition mechanism simulated by MD. Baez and Clanzy simulated the decomposition of hydrate clusters, and the results showed that in the process of dissolution, the crystal cavity located at the interface lost part of its structure, leaving only part of dodecahedron and tetrahedral cage. Before the complete decomposition, the remaining structure was a dodecahedron cavity, which could be maintained for a period of time [80]. English et al. simulated the decomposition process of methane hydrate crystals surrounded by liquid phase at 276.65 K and 68 bar. The liquid phase is composed of pure water or methane-water mixture. The composition of methane is 50% to 100% of the theoretical maximum value of the corresponding hydrate. Four different types of crystals are used in the two-phase system. These crystals include cavities containing and without methane, where the cavities containing methane contain 80 to 100 percent of the theoretical maximum. The simulation results showed that when the methane content is within the range of 80% ~100%, the methane content in the hydrate phase has little influence on the decomposition rate. In contrast, the decomposition rate of empty hydrate clusters is faster. In contrast, the decomposition rate of empty hydrate clusters is faster. In addition, the diffusion of methane molecules from the crystal surface to the surrounding liquid layer may be a major control factor for the decomposition rate of hydrate [72]. English and Clarke studied the thermal decomposition of CO2 hydrate interfaces in 300~320 K liquid water using equilibrium and nonequilibrium MD simulations. The results showed that the decomposition rate is closely related to temperature, and the higher the temperature is, the faster the decomposition rate will be. The mass and thermal coupling model is used to observe and distinguish two different decomposition mechanisms. The second well-defined region is independent of composition and temperature, in which the lattice framework of the remaining two-dimensional nanoscale systems is inherently unstable [81]. Ding et al. simulated the decomposition process of methane hydrate under the conditions of , , 320, and 325 K under the NPT ensemble. Figure 6 showed snapshots of the continuous change process of methane hydrate from crystal to liquid. As can be seen from the figure, the decomposition process of methane hydrate is composed of two continuous stages. The first stage is the process of molecular diffusion, which leads to the increase of cell size and the fracture of lattice structure. In the second stage, methane molecules escape from the broken hydrate cavity and gather together. Among them, the fracture of lattice structure is the key step of hydrate decomposition [82]. Myshakin et al. carried out simulation studies on the decomposition process of methane hydrate with initial methane cage occupancy of 100%, 95%, and 85% at 6.8 MPa. The results showed that the decomposition rate was closely related to the hydration number, and the decomposition rate of the cage showed the behavior of Arrhenius, which was consistent with its formation mechanism. It is worth mentioning that when the temperature is higher than the calculated hydrate decomposition temperature, part of the water cage formed in reverse around the methane molecule [83]. Iwai et al. studied the decomposition process of methane and carbon dioxide hydrate, respectively, under the NPT ensemble. TraPPE unit point and 5-point force field model were used for methane, and EPM2 (3-point) and SPC/E models were used for carbon dioxide and water molecules, respectively. The simulated temperature and pressure were 270 K and 5.0 MPa. The results showed that the methane molecular potential energy model has little influence on the decomposition process of methane hydrate. In the decomposition process, the water cage first breaks down, and then the gas molecules escape from the cage. Under calculation conditions, methane hydrate is more stable than carbon dioxide hydrate [84]. Smirnov and Stegailov conducted MD simulations of pressures up to 5000 bar for different models of water molecular force fields. They calculated the dynamic stability boundary of the superheated metastable type I structure and analyzed the effects of heating rate, system size, and cage occupancy on the stability of the system. The results showed that a perfect type I lattice can withstand great overheating, the dynamic stability boundary is generally dependent on the heating rate and system size, and the reduction of cage occupancy reduces the decomposition temperature systematically [85].

4.2. Gas Swapping

In recent years, the important prospect of carbon dioxide storage in clathrate hydrates has attracted much attention, and many researchers have proposed the use of carbon dioxide to replace methane in frozen soils or deep oceans. However, various technical and operational challenges remain in this process. Releasing large amounts of methane into the atmosphere, for example, can cause serious harm to the environment. In order to extract natural gas from hydrate reservoirs and realize CO2 burial, scholars proposed the original idea of replacing CH4 with CO2 in NGH [8688]. Smith et al. studied the possibility of substitution process [89], and Ohgaki et al. [86] showed that the equilibrium pressure of CO2 hydrate was lower than that of CH4 hydrate at the same temperature. In addition, the heat of formation of CO2 hydrate can provide heat for the decomposition of CH4 hydrate, which solves the problem of lack of heat source in the development of NGH [86]. Recent experimental studies have shown that CH4 can be extracted from CH4 hydrate by reacting with CO2 [9093]. The thermodynamic calculation results also confirm the feasibility of the hydrate exploitation by CO2 replacing CH4 gas under appropriate conditions [94].

To investigate the stability of CH4, CO2, and CH4-CO2 mixed gas hydrate, Geng et al. carried out MD simulation and stability energy calculation. A type I CH4 hydrate, CO2 hydrate, and CH4-CO2 mixed hydrate simulation system with 100% cage occupancy was constructed in a unit simulation box with periodic boundary conditions. The simulation results showed that CO2 molecules are more suitable for large cavities than CH4 molecules, and CH4-CO2 mixed gas hydrate is the most stable, that is, CH4-CO2 mixed gas hydrate may be formed through the replacement of CH4 in the hydrate [95]. Tung et al. conducted MD simulation of CO2 displacement of CH4 gas under conditions of 280 K and 6 MPa, and displacement snapshots at different simulation times are shown in Figure 7. According to the distance between liquid water and solid hydrate, the replacement mechanism can be divided into two types, namely, direct exchange of methane and carbon dioxide and transient coexistence of methane and carbon dioxide in a cavity. During the displacement process, there are two types of cage openings through which carbon dioxide/methane molecules enter/leave the cage. The two openings are located at the junction of hydrogen bond, hexagon, and pentagon and the adjacent of the two pentagons, respectively (see Figure 8) [96]. Qi et al. conducted an interface MD simulation similar to Tung et al. The results showed that in CH4 hydrate, the area near the interface between gas phase and hydrate is the most likely place where CH4 can replace CO2 [97]. Bai et al. studied the effects of memory effect, mass transfer limitation, and chemical potential of guest molecules on replacement. Their research results showed that the amorphous layer greatly hindered the mass transfer of guest molecules, and the presence of the amorphous layer reduced the replacement rate between CO2 and CH4 and prevented the further decomposition of CH4 hydrate [98].

4.3. Inhibitor

Under certain temperature and pressure conditions, NGH forms in formation or oil and gas transportation pipelines, which causes clogging and has adverse effects on the economy and even the environment. In order to avoid this problem, many scholars began to carry out studies on NGH inhibitors [99102]. Inhibition can now be achieved in two ways: thermodynamic or kinetic. Thermodynamic inhibitors, such as methanol or ethylene glycol [103], move the three-phase equilibrium boundary to more extreme positions, reducing the possibility of hydrate formation. The technology is predictable and reliable, but often requires extensive use, is therefore costly, and may cause some environmental harm. Kinetic inhibitors can delay the nucleation process of hydrate or reduce the growth of hydrate crystals [104, 105], and their dosage is generally much smaller than that of thermodynamic inhibitors [106, 107], so they have a broader application prospect.

Many scholars have tried to screen suitable kinetic inhibitors, but the progress was slow because of the lack of understanding of the working mechanism of inhibitors. Of the several inhibitory mechanisms available, the most common is that kinetic inhibitors irreversibly bind to the surfaces of small hydrate crystals and critical clusters in a manner similar to antifreeze proteins [1]. In fact, this kinetic inhibitor probably promotes hydrate nucleation, because molecules that favor surface binding also create seed sites for nucleation, but then prevent the nanocrystals from growing to the observable size. These compounds have been classed as nucleating inhibitors because they experimentally delay the formation of hydrates, although their activity is actually produced by inhibiting crystal growth [108]. Another mechanism is reversible binding to the hydrate surface or to the hydrate/water interface region. In this case, changes in water structure were observed to destabilize small hydrate clusters, thus, increasing the time to induce nucleation [90]. Selective dissolution of guest molecules to inhibitors may also play a role. In fact, it is possible that multiple mechanisms are available, and the success of the synergistic inhibitor mixture is the result of the interaction of different inhibitors in the mixture with many different hydrate nucleation and growth mechanisms.

Amino acids are a new class of hydrate kinetics inhibitors, which are considered as biodegradable and nontoxic additives [109], and their inhibitory mechanism has been widely concerned. Amino acids contain two groups: a hydrophobic (side-chain) group and a hydrophilic (carboxylic, amino) group, which forms strong hydrogen bonds with water. Li et al. simulated the inhibition effect of PVP-A, PVP-E, and PVP on the growth of hydrate, and the results showed that the inhibition effect of PVP-A was the best. The inhibitor combined with liquid water molecules on the surface of hydrate nucleus to form hydrogen bonds, preventing the growth of hydrate nucleus from reaching the critical size [110]. Kuznetsova et al. carried out the polyethylene caprolactam (PVCap) impact on the formation and decomposition of methane hydrate research. The results showed that the PVCap behavior is closely related to its concentration in the water. When the concentration is low, PVCap tends to remain at the methane-water interface, and there is no interaction between them; when the content of PVCap is high, the interface tension between methane-water is significantly reduced, and separate methane bubbles are formed in the randomly dispersed methane phase [111]. Maddah et al. used MD to simulate the growth kinetics of methane hydrate at different concentrations (<1.5 wt.%) of glycine, alanine, proline, and serine. The results showed that all amino acids had inhibitory effects on methane hydrate, and the inhibitory effects were sequentially ordered as , and the inhibitory mechanism was the interaction between amino acids and solution that disturbed the hydrate structure [112]. Natural inhibitors, such as glue, chitosan, cassava powder, and antifreeze protein, also have certain inhibitory effects on hydrates. However, there are few studies on the inhibitory mechanism of these natural products. Pectic polysaccharide is a kind of nontoxic, natural products. Xu et al. studied the mechanism of pectin inhibited the growth of the methane hydrate by MD simulation and gave the pectin under different simulated times when the mass fraction of 2.46% and 3.62%, respectively, of methane hydrate generation process conformation (Figure 9). As can be seen from Figure 9, at the initial stage of simulation, methane gas accumulates around pectin, and the gas-liquid phase splits, forming a pectin layer between water molecules and guest molecules, increasing the mass transfer resistance, and resulting in difficult hydrate growth. At 20 ns, water molecules interact with pectin through hydrogen bonds, disrupting the hydrogen bond network between water molecules and hindering the further growth of gas hydrate (Figure 9(b), 20 ns enlarged figure) [113]. Zi et al. conducted MD simulation to explore the role of asphaltenes in the formation of methane hydrate, and the results showed that the asphaltenes located at the gas-water interface promoted the formation of hydrates, while the asphaltenes in water had a slight inhibitory effect on the growth of hydrates [114].

5. Conclusions

On the basis of previous studies, this paper briefly summarized the simulation process of molecular dynamics, discussed the applicability of each molecular force field in the hydrate system, and summarized the formation and decomposition rules of natural gas hydrate. The following are the main understandings:

(1) The force field model is a model to describe the force exerted on molecules. In the simulation of NGH, the choice of water molecule model has a great influence on the final simulation results. Generally speaking, it is difficult for a single water molecule model to describe the crystallization and decomposition process of hydrate well at the same time. Therefore, different water molecule models are needed to simulate separately. In the future research, it is still an open research field to explore the water model applicable to the study of NGH in different thermodynamic states

(2) The research on the nucleation mechanism of hydrate at the gas-water interface shows that under the conditions of supercooling, high pressure, and gas supersaturation, hydrates are more easily formed. However, the experimental results show that the surface properties of the sediments (such as crystallinity and hydrophilicity) may also influence the formation of hydrates. Therefore, it is of great significance to strengthen the characterization of the surface properties of hydrate reservoir and study the gas-liquid-solid interface effect for the application of MD simulation in the study of natural gas hydrate

(3) The extraction methods of NGH mainly include depressuring, thermal stimulation, gas swapping, and chemical reagent injection. Because it is difficult to reduce the pressure of the hydrate system to realize the hydrate decomposition process, there are few reports on the decomposition mechanism simulated by MD. The MD simulation results of hydrate decomposition at different temperatures show that the decomposition rate is closely related to temperature, and the higher the temperature is, the faster the decomposition rate will be. The gas swapping method can not only realize carbon dioxide storage but also hydrate mining without damaging the geological structure of the reservoir. However, its research is still in the theoretical stage, and further research should be carried out in terms of replacement location and replacement efficiency in the future

(4) Hydrate inhibitors are of great significance in preventing hydrate from blocking pipes, which can be divided into thermodynamic inhibitors and kinetic inhibitors. Kinetic inhibitors are of low dosage and have a broader application prospect. Amino acids are a new class of hydrate kinetics inhibitors, which are considered as biodegradable and nontoxic additives. The hydrophilic (carboxylic, amino) groups of amino acids can form strong hydrogen bonds with water and interfere with the growth of hydrate. The inhibitory effect of amino acids is closely related to their type and concentration. Generally speaking, the higher the concentration is, the more obvious the inhibitory effect is. Screening efficient and economical kinetic inhibitors and finding their mechanism of action are still the long-term research topics of molecular dynamics simulation

Conflicts of Interest

The authors declare that they have no conflicts of interest.


The authors are grateful to the National Natural Science Foundation of China (51991365), Key Special Project for Introduced Talents Team of Southern Marine Science and Engineering Guangdong Laboratory (Guangzhou) (GML2019ZD0105), Key Program of Marine Economy Development (Six Marine Industries) Special Foundation of Department of Natural Resources of Guangdong Province ([2020]047), and China Geological Survey Project (No. DD20211350).