Modeling and molecular simulation of natural gas hydrate stabilizers

Uncontrolled decomposition of natural gas hydrate may lead to serious marine geological disasters and air pollution. The model of natural gas hydrate and lecithin was established. The stability mechanism of lecithin to structure hydrate was studied by molecular dynamics simulation. The consistent valence force ﬁ eld (CVFF) and TIP3P potential models are used to de ﬁ ne the interaction between CH 4 -CH 4 and water-water species, respectively. The simulations are performed on a combination of a 2 × 2 × 4 unit cell of sI hydrate and a water liquid phase with lecithin. The results of the simulations indicate that lecithin molecules adsorb on the hydrate surface with their hydrocarbon chains crossing and forming a net structure, easily producing the hydrate memory e ﬀ ect, which will narrow the available space for hydrate methane and water movement. Compared to the pure water-hydrate model, the mean square displacement (MSD) values of hydrate methane and water molecules are much lower, indicating that the hydrate dissociates more slowly.


Introduction
Natural gas hydrate (NGH) has great potential as a future alternative fuel because of its large reserves, high energy density, and low environmental pollution (Sloan & Koh, 2007).Considering the fact that the vast majority of global NGHs exist in the shallow layer of deep water areas (Booth et al., 1999;Brown et al., 1996;Miles, 1995), the drilling of gas hydrate resources is the most direct and effective means for NGH exploration and development.However, hydrate is a safety hazard in drilling operations.When drilling in NGH sediments, the changes in temperature and pressure conditions as a due to in-situ formation disturbance can result in hydrate decomposition, thereby resulting in a series of drilling problems, such as borehole wall instability (Tan et al., 2005), gas-cutting trouble, and drilling fluid performance deterioration, which can result in submarine landslides and serious environmental pollution (Glasby, 2003;MacDonald, 1990;Paull et al., 1996).
In the cascade glacier fields of the Alaskan North Slope (Schofield et al., 1997) or the shallow strata of certain deep water areas, the structure of the hydrate layer is loose and the mud weight window is very narrow and has a low value.Therefore, the pressure control and temperature-pressure control methods, which are two types of widely accepted NGH layer drilling methods, cannot ensure rig safety in this type of formation.In recent years, BP has successfully solved this problem by adding a certain amount of lecithin to the drilling fluid when drilling a K-13 well (Schofield et al., 1997).The research results of Schofield, Jiang, and Ning generally expound that the inhibitory effect of lecithin on the hydrate decomposition was mainly caused by its adsorption onto the surface of the hydrates (Jiang et al., 2001;Schofield et al., 1997).Schofield et al. (Chen et al., 2007) proved that 0.3%wt of lecithin did not affect the hydrate formation thermodynamic conditions of the drilling fluid through a series of laboratory experimental studies.According to the results of these experimental studies, we can roughly determine that there is an optimal use concentration of lecithin, which also provides some practical support for our subsequent research on its inhibition mechanism.In spite of this, studies on hydrate chemical stabilizers are rare and most focus on experimental research.
Therefore, in this work, lecithin molecular adsorption processes on the sI hydrate surface and hydrate crystal dissociation are studied by the molecular dynamics simulation method.The aim of this work is to explain the hydrate deformation process and the mechanism of lecithin to stabilize hydrate crystals.

Simulation methodology
The molecular dynamics method simulates the behavior of lecithin adsorption on the methane hydrate surface as well as hydrate dissociation.The methane molecules were enclosed in both small and large types of cages.All of the molecular dynamics simulations with constant number, volume, and temperature (NVT) and constant number, pressure, and temperature (NPT) ensembles were performed using the Fotcite module of Accelrys Materials Studio (Alalwiat et al., 2005) software package.

Intermolecular potential models
The well-known TIP3P model is used for the waterwater interactions.The TIP3P model is a three-point model that has superior performance over the SPC water model when describing the energy and density of liquid water (Jorgensen et al., 1983).The intermolecular potential for the interaction of a pair of molecules is represented in the following form: (1) where i and j represent the respective charge sites, qi and qj are the point charges, ro-o is the distance between two oxygen atoms, and A and C are the potential parameters.All of the parameters of the TIP3P water model are presented in Table 1.
The consistent valence force field (CVFF) (Dauber-Osguthorpe et al., 1988) was employed to describe the water-methane, water-lecithin, and methane-lecithin interactions.Within the CVFF model, the non-bond interactions, including the van der Waals and electrostatic interactions, were represented by the Lennard-Jones (L-J) function and Coulombic representation, respectively.The equation of these interactions has a similar form as Equation (1), as presented in Equation (2) as follows: where m and n are the monomers in the simulation model, and Aij and Bij are the pair parameters.
Compared to the L-J function, A ii = 4ε i σ 12 and B ii = 4ε i σ 6 .Berger et al. (1997) adopted a simplified force field model, and simulates the DPPC (a phospholipid molecule, and it is very similar to lecithin molecule) the bilayer characteristics.The results agree well with the experimental results, which proves that this set of parameters can effectively describe the interactions between molecules.All electronic and pair parameters of the same atoms are listed in Table 2.
The Lorentz-Berthelot combination rules (Reid et al., 1987), as presented in Equation (3), are employed to calculate the interactions between different monomers.

Simulation method
Our molecular dynamics study involved two molecular dynamics simulation models.To study the effect of lecithin enhancement on hydrate stability, the first simulation model, which was later denoted as Model I, was built as a hydrate-liquid, water-lecithin model comprised of a 2 × 2 × 4 structure I hydrate unit cell filled with methane molecules (contains 736 TIP3P water molecules and 128 CVFF methane molecules), and a liquid phase composed of 810 TIP3P water and 3 lecithin molecules (Figure 1).Two of the lecithin molecules, hereafter referred to as Lecithin A and B, were placed near the right and left surfaces of the liquid phase, respectively, and the third molecule (referred to as C) was placed between Lecithin A and B in the liquid phase.The initial arrangement of the liquid water molecules was determined from the conditions of the minimum energy.The hydrate layer was placed on the left side of the liquid layer.The hydrate crystals were built according to the guidelines (Storr et al., 2004).In the hydrate layer, the initial positions of the oxygen atoms were determined from an X-ray single crystal diffraction experiment (Hollander & Jeffrey, 1977;Kirchner et al., 2004), and the hydrogen atoms of the water molecules were arranged disorderly.Every hydrate cage contained a methane molecule at its center.Simulation model II was the same as Model I except that it contained no lecithin molecules in its liquid layer.Periodic boundary conditions (PBC) were applied in all three directions.Long-range Coulombe forces

R E T R A C T E D
were treated by the Ewald summation method .The temperature was controlled according to the thermostats procedure of Nosé and the Nosé-Hoover dynamics (Fincham, 1992;Hoover, 1985;Shuichi, 1991) at 277 K, and the control parameter Q ratio was set at 1.0.The cut-off radius is 10 angstrom.The simulation began from this initial model and 10,000 time steps (step size 0.5 fs) were employed in the NPT simulation to obtain a stable start configuration for the NVT simulation.The pressure was controlled at 2.5 MPa with a Berendsen barostat (Berendsen et al., 1984) and the decay constant was set at 0.1 ps.
The final conformation of the NPT simulation was used as the initial conformation for the NVT simulations.To study the adsorption characteristics of the lecithin molecules on the hydrate surface, all atoms in the Model I hydrate crystal were constrained and the NVT simulation executed 1,600,000 time steps (step size of 0.5 fs).In its final conformation, another set of NPT simulations, wherein all of the hydrate atoms were unconstrained, was performed for a total time of 1000 ps (2,000,000 steps, step size of 0.5 fs) for the hydrate dissociation.For the simulations of Model II, which was taken as a matched group, a set of NPT simulations was performed and all of the computational details were the same as Model I.
Prior to initiating the dissociation simulations presented in this paper, we built a 2 × 2 × 2 sI hydrate crystal model using the same parameters as Models I and II.Four sets of simulations were performed to find the equilibrium pressures under different temperatures.The equilibrium pressures of the pure guest (methane) were described according to the experimental data (Frost & Deaton, 1946).To demonstrate the consistency of the obtained data with experimental data of each temperature, a semilogarithmic plot for the pressure versus temperature is plotted in Figure 2. As shown in the figure, the simulation results are in good agreement with the experimental data.

Adsorption characteristics of lecithin on the hydrate surface
To study the adsorption effects of lecithin on the hydrate surface, the hydrate crystal of Model I was constrained and NVT simulations were performed.The initial conformation was derived from the first NPT simulations.The dimensions of the model were 20.77Å × 20.77 Å × 75 Å.The simulation results showed that the morphology of lecithin was changing constantly.Figure 3(a-c) represent the snapshots of Model I at the beginning, middle, (400 ps), and end of the simulation at 800 ps, respectively.To facilitate the following analysis, the three lecithin molecules placed in the left, right, and middle of the water phase will be referred to as lecithin A, B, and C, respectively.
The aliphatic oxygen atoms of lecithin A formed a hydrogen bond with hydrate water, whereas the oxygen and hydrogen atoms of the hydroxyl groups (phosphate and choline groups) formed hydrogen bonds with the liquid water.Both of the hydrocarbon chains of the fatty acids laid flat on the hydrate crystal surface and exhibited molecular conformational changes.Similar to lecithin A, hydrogen bonds were formed between the aliphatic oxygen atoms of lecithin B and the water on the hydrate surface.At the beginning of the simulations, the

R E T R A C T E D
In summary, the lecithin molecules adsorbed onto the hydrate surface and their phosphate and choline groups attracted other nearby lecithin molecules.The distortion of the hydrocarbon chains may have resulted in a net structure.

Gas hydrate dissociation characteristics
In this section, we calculated the radial distribution function (RDF) profiles, average atomic mass densities on the z-axis direction, MSD profiles, and diffusion coefficients of both models.The hydrate dissociation characteristics under the influence of the lecithin molecules were examined by comparing Models I and II.

Snapshot structures and the radial distribution function (RDF)
The radial distribution function (RDF) (Hansen & McDonald, 1990)  Figure 4 presents the initial conformation of Models I and II.The hydrate crystals of both models in the simulation, which presented the initial moment, maintained the clathrate structure.
After 6 ps, the hydrate in Model I remained structurally intact and exhibited a small deformation near the hydrate-liquid interface.Meanwhile, the morphological changes of the three lecithin molecules were small and their positions remained unchanged.The deformation of the hydrate crystals in Model II was significantly larger than those in Model I, where the first hydrate crystal layer of Model II dissociated by half.This result indicates that the clathrate structure had been opened and the methane molecules could easily diffuse into the liquid phase.In Figure 5(b), the shapes of the RDF profiles for the hydrate O-O and C-C atom pairs in Model I were nearly identical to the profile at 0 ps except for the peak height.The reduced height of the peaks compared to the profile at 0 ps was likely attributed to the slightly larger thermal motion of the water molecules near the interface.The RDF profiles of Model II exhibited the same characteristics as the original curve except for the lower peak height.
When the simulation evolved to 100 ps, the first hydrate crystal layer in Model I was almost completely dissociated, whereas the remainder of the crystal maintained a clear clathrate structure.The three lecithin molecules exhibited constantly changing conformations during the simulations but remained on the hydrate surface.At the same time, the second hydrate layer was opened.The dissociation rate differences may be attributed to the two lecithin layers on the hydrate surface, which attracting the water molecules dissociated from the hydrate crystal, thereby reducing the amount of space available for movement for the hydrate water and the methane molecules.In Model I, the O-O pair RDF profile (Figure 6(b)) represented the same characteristic as the profile at 0ps except for the lower heights of the first, second, and third peaks and the disappearance of the remaining peaks.In Model II, the characteristic peak heights of the O-O and C-C pair RDF profiles exhibited a continuous decrease.In addition, some small fluctuations were observed in front of the first peak in both C-C pair RDF profiles, thereby indicating the gathering trend of the free methane molecules.
At 1000 ps, t he second hydrate layer in Model I exhibited serious deformation.In addition, the clathrate structure was opened but not completely dissociated, whereas the third hydrate layer in Model II dissociated completely.As compared to the O-O pair In Model I, the hydrate dissociation rate of the second period (100 ps to 1000 ps) was much lower than that of the first period (0 ps to 100 ps).In the first period, the first hydrate layer dissociated completely.In the second period, only the clathrate structure of the second layer was opened.The lecithin molecules were absorbed onto the interface and reduced the amount of space available for free water and methane molecule movement, thereby reducing the dissociation rate.In Model II, no obvious difference was observed in the dissociation rates.

Atomic and molecular density profile analysis
To obtain the density distribution profiles, we first isometrically divided both Models into 600 along the (0 0 1) direction (from hydrate to liquid phase).Then, we calculated the relative concentration of all types of atoms with the Forcite module in the MS package.The relative concentration is defined as follows: where C R,i is the relative concentration of atom i; N i and N i,total are the number of atoms i in the slab and the total model, respectively; and V slab and V total are the volume of the slab and total model, respectively.
To observe the density distribution characteristics of all types of atoms and molecules, we converted the relative concentrations to density data.The conversion formula is presented as follows:

R E T R A C
T E D where M i is the relative atomic weight of atom i; ρ i , ρ m,i , and ρ total are the densities of atom i, molecule and total slab, respectively; ρ j and ρ k are the densities of the atoms in the molecule and the slab, respectively.Corresponding to Figure 4, the atomic density profiles of Model I (Figure 8(a)) and Model II (Figure 8(b)) represent the sI methane hydrate characteristic structure.The curves in both Figure 8(a-b) were in the range of 0 ~47.5 Å (z-direction) and indicate the hydrate characteristic density profile, whereas the remainder represents the liquid phase.In the hydrate, the methane profiles (black lines) and oxygen profiles (blue lines) exhibited an obvious corresponding relationship.
The high-density methane peaks coincided with the low valleys of the oxygen density, whereas the lowdensity peaks corresponded to the high valleys between the high-density peaks of the oxygen profiles.The regular pattern of the four hydrogen peaks (red line) also reflected the characteristics of the hydrate crystal structure: the two coincided with the methane peaks, whereas the other two were located at the high valley position of the oxygen profile.The highly structured hydrogen pattern was not observed in the liquid phase.Thus, this characteristic can be used to determine the complete dissociation and restructuring of the hydrate (if necessary).Although the lecithin molecular density profiles (magenta line) exhibited no obvious structural features, they can be used to indicate the position of the lecithin molecules.
As presented in the configuration snapshots (Figures 4-7) and the density distribution profiles , the lecithin molecule at the left of the liquid phase exhibited constant deformed but did not move during the dissociation simulation process.The two lecithin molecules at the right side moved to the right (left interface of the hydrate) prior to 100 ps and then maintained its position for the remainder of the simulation.From 0 ps to 100 ps, the hydrate in Model I dissociated faster than in Model II, though only a small difference was observed between both models.In the next period (100-1000 ps),

R E T R A C T E D
restructured hydrate layer reduced the allowed movement space for the water and methane molecules that were newly dissociated from the hydrate, thereby significantly reducing the hydrate dissociation rate in Model I.This rapid restructuring was consistent with the distinguishing features observed in the hydrate experiments, or the so-called "hydrate memory" effect (Ohmura et al., 2003).This effect refers to the phenomenon in which water molecules that previously formed a hydrate more easily and rapidly reformed hydrate crystal following the fulfillment of the thermal equilibrium condition as compared to water molecules that never formed a hydrate structure.The phenomenon that the induction period for the formation of hydrate after the decomposition of hydrate is significantly reduced is called hydrate memory effect (Makogon, 1981).Observations of the hydrate memory effect increased the credibility of our simulation results.In contrast, neither hydrate conformation nor density profile characteristics were observed in the dissociated area in Model II, thereby indicating the absence of hydrate crystal reformation during that period and an undisturbed decomposition process.
Mean square displacement (MSD) and diffusion coefficients Figure 12(a) presents the MSD of the water molecules and methane molecules in Models I and II.In both models, the MSDs of the host water molecules and guest methane molecules increased in accordance with the simulations.With respect to the MSD amplitude,   where R i (t) is the MSD of the type i molecules at time t and D is the diffusion coefficient.As presented in Figure 12(b), the diffusion coefficients of the water and methane molecules of Model I were 2.4693 × 10 −10 m 2 /s and 1.1250 × 10 −10 m 2 /s, respectively, whereas those of Model II were much larger, reaching 4.3195 × 10-10 m2/ s (H 2 O) and 2.1053 × 10 −10 m 2 /s (CH 4 ).
Combining the results of the above analyses (snapshot structures, RDF profiles, atomic & molecular density profiles, MSDs, and diffusion coefficients), the dissociation process of hydrate with the nearby lecithin molecules were roughly established.At the start of the dissociation, the released water and methane molecules had reduced movement space due to the presence of the lecithin molecules, which resulted in slower decomposition.The hydrate memory effect was produced much earlier under the continued dissociation as compared to those without lecithin molecules due to the presence of more molecules that have previously formed hydrate in the area between the lecithin molecules and hydrate surface.The new hydrate layer between lecithin and hydrate further limited the movement of the water and methane molecules.In addition, the hydrate dissociation rate continued to decrease.From a macroscopic aspect, the hydrate crystal dissociation rate can be greatly reduced by lecithin in the liquid phase.The results of the laboratory evaluation (Chen & Kamath, 2006;Jianguo et al., 2010;Yan et al., 2012) and field application (Schofield et al., 1997;Uchida et al., 1999) also confirmed that lecithin is a very effective type of hydrate stabilizer (or promoter).Thus, the simulation results in this paper conformed with the objective laws.

Conclusion
In this work, the molecular dynamics simulations (with CVFF) of a lecithin solution-hydrate model were divided into the following two processes: adsorption and dissociation.As compared to the MDS of the pure water-hydrate model, the influence of lecithin on the dissociation of the hydrate was examined.The results can be summarized as follows: (1) The lecithin molecules adsorbed onto the hydrate surface, and their phosphate and choline groups attracted other nearby lecithin molecules.The distortion of the hydrocarbon chains resulted in the formation of a net structure.
(2) As compared to the pure water-hydrate model (Model II), the hydrate dissociation rate was always lower in the presence of lecithin molecules near the interface.Initially, the movement space of the released water and methane molecules was reduced by the presence of the lecithin molecules near the water-hydrate interface, which resulted in slower hydrate dissociation.The hydrate memory effect was also observed much earlier and a new hydrate layer was formed because most of the water and methane molecules between the lecithin molecules and the hydrate surface were released from the hydrate.Hence, the space available for movement was further compressed and a continuously decreasing hydrate dissociation rate was observed.

R E T R A C T E D
Based on the above conclusions, subsequent research and development of natural gas hydrate stabilizer can further explore the interaction mechanism between lecithin and its similar structures and natural gas hydrate, with a view to developing a hydrate stabilizer with better performance.

Disclosure statement
No potential conflict of interest was reported by the authors.

R E T R A C
T E D

Figure 1 .
Figure 1.Molecular structural formula of the lecithin molecule.

Figure 2 .
Figure 2. Pressure versus temperature plot showing a comparison between the MD simulation results and the experimental data.

Figure 3 .
Figure 3. Evolution of Model I. Water molecules in liquid phase are hidden for the convenience of observing lecithin molecule formation and position changes.The water molecules in NGH are white and red sticks, the methane molecules are gray and white sticks, and the lecithin molecules are gray, red, pink, and blue.The light blue dotted line represents the hydrogen bonds.(a) Starting configuration.(b) Configuration of 400 ps.(c) Last configuration.
reflects the characteristics of the crystal microstructure and can serve as an index of hydrate stability.To quantitatively analyze the dissociation process, two major RDFs of O-O atom pairs and C-C pairs at a series of times are shown in Figures 5(b), 6(b), and 7(b).

Figure 4 .
Figure 4. Snapshot of Model I (the upper one) and Model II.Same conventions in Figure 3 apply.

Figure 6 .
Figure 6.Model conformation and RDF curves at 100 ps.(a) Model conformation; (b) RDF curves.R E T R A CT E D

Figure 8 .
Figure 8. Mass density profiles at 0ps.The black line corresponds to carbon; blue to oxygen; red to hydrogen; magenta to lecithin; and olive is the net density.

Figure 12 .
Figure 12.(a) MSDs and (b) diffusion coefficients of the host water molecules and guest methane.

Table 2 .
Pair-parameters in lecithin model.