Does Confinement Enable Methane Hydrate Growth at Low Pressures? Insights from Molecular Dynamics Simulations

Natural methane hydrates are estimated to be the largest source of unexploited hydrocarbon fuel. The ideal conditions for methane hydrate formation are low temperatures and high pressures. On the other hand, recent experimental studies suggest that porous materials, thanks to their confinement effects, can enable methane hydrate formation at milder conditions, although there has not been a consensus on this. A number of studies have investigated methane hydrate growth in confinement by employing molecular simulations; however, these were carried out at either very high pressures or very low temperatures. Therefore, the effects of confinement on methane hydrate growth at milder conditions have not yet been elaborated by molecular simulations. In order to address this, we carried out a systematic study by performing molecular dynamics (MD) simulations of methane water systems. Using a direct phase coexistence approach, microsecond-scale MD simulations in the isobaric–isothermal (NPT) ensemble were performed in order to study the behavior of methane hydrates in the bulk and in confined nanospaces of hydroxylated silica pores at external pressures ranging from 1 to 100 bar and a simulation temperature corresponding to a 2 °C experimental temperature. We validated the combination of the TIP4P/ice water and TraPPE-UA methane models in order to correctly predict the behavior of methane hydrates in accordance to their phase equilibria. We also demonstrated that the dispersion corrections applied to short-range interactions lead to artificially induced hydrate growth. We observed that in the confinement of a hydroxylated silica pore, a convex-shaped methane nanobuble forms, and methane hydrate growth primarily takes place in the center of the pore rather than the surfaces where a thin water layer exists. Most importantly, our study showed that in the nanopores methane hydrate growth can indeed take place at pressures which would be too low for the growth of methane hydrates in the bulk.


■ INTRODUCTION
Gas hydrates, also known as clathrate hydrates, are a subset of nonstoichiometric crystalline inclusion compounds that are formed when the self-assembly of water molecules into a 3D hydrogen-bonded framework of cavities enclathrates small gas molecules. The ideal conditions for gas hydrate formation are usually low temperatures (<300 K) and high pressures (>6 bar), and their structure is stabilized by van der Waals forces. To date, more than 130 molecules (or hydrate formers) have been identified that form gas hydrates. 1 There are three common crystalline structures of gas hydrates, namely structure I (sI), structure II (sII), and structure H (sH). The type of structure they adopt is determined by a range of factors, i.e., the formation conditions and the type and the size of the guest molecules that are enclathrated. One of the key properties of gas hydrates is their extraordinary gas storage capacity. At full occupancy (i.e., all cages are fully occupied), 1 m 3 of a gas hydrate can store up to 173 m 3 (STP) of gas. 2 In the past few decades, there has been a surge of interest in gas hydrate research due to their relevance to flow assurance, 3,4 global warming, 5−7 and marine geohazards. 8,9 Gas hydratebased technologies have been proposed in various fields, including but not limited to gas mixture separation, 10 energy recovery, 11 and gas storage and transportation. 12 The biggest challenges in exploiting gas hydrate-based technologies are their slow formation and dissociation kinetics and a poor understanding of the formation and dissociation mechanisms of gas hydrates. There is a considerable amount of literature on gas hydrate formation and dissociation in the bulk phase, with or without the presence of impurities such as hydrate promoters and inhibitors, via experiments 13−19 and molecular simulations. 20−28 Consequently, the phase equilibria and other thermophysical properties of bulk gas hydrates are well established. On the other hand, the nature of hydrate formation and dissociation in a confined nanospace, which is of crucial relevance to understanding the appearance of natural gas hydrates in complex porous environments, is a matter of ongoing scientific discussion.
Methane hydrates are the most abundant form of gas hydrates and are typically found naturally in permafrost regions and marine sediments. Estimates suggest that the amount of energy stored as natural methane hydrates is at least twice that of all other hydrocarbon-based fuels combined, making it the largest source of unexploited fuel. 4,29 They are of the cubic sI crystal type, where each unit cell is comprised of two "small" dodecahedron cages (denoted as 5 12 ) and six "large" tetrakaidecahedron cages (denoted as 6 2 ) with coordination numbers of 20 and 24, respectively. An average cavity radius of 3.95 and 4.33 Å for the small and large cages, respectively, is sufficient to accommodate one methane molecule per cage. Hence, a nominal stoichiometry of 1CH 4 :5.75H 2 O is obtained for methane hydrates at full occupancy.
Recent studies have focused on the provision of solid surfaces in gas hydrate studies, both experimentally and computationally. An improved understanding of methane hydrate formation in confinement is essential for a wide range of scientific and industrial purposes, especially in the advancement of hydrate-based technologies. To date, there has been little agreement about the precise nature of the confinement effects. Some studies 30,31 have suggested that the effect of inhibition is stronger than that of promotion in a confined space due to the reduced water activity. However, the promoting or inhibiting effects of the porous materials are still largely unknown, owing to the experimental challenges associated with quantifying hydrate formation comprised of nucleation and growth in such systems where the length and time scales involved are often inaccessible experimentally. In order to address this, researchers have employed computational approaches, such as molecular dynamics (MD) and Monte Carlo (MC) simulations, together with experimental studies in order to study gas hydrate systems.
Cha et al. 32 found that methane hydrate formation is promoted in the presence of bentonite, for which they suggested that the clay mineral surface could act as a nucleation site by forming an ordered layer of adsorbed water on the surface. Conversely, Handa and Stupin 33 reported the thermodynamic inhibition of confined methane and propane hydrates in silica gel with an average pore diameter of 15 nm as a result of reduced water activity. Similar observations can also be noted in the work of Uchida et al., 34 who used porous glasses with a pore size distribution ranging from 10 to 50 nm. Several studies 35−37 have utilized molecular simulations in order to study the promoting effects of Na montmorillonite in regard to methane hydrate formation, which showed that clathrate-like water structures enclosing the methane molecules near the clay surfaces could stabilize the confined methane hydrates. One study by Bai et al. 38 reported the nucleation mechanism of carbon dioxide hydrates in the presence of hydroxylated silica surfaces, which is a three-stage process; ice-like structures trigger the formation of intermediate hydrate structures, which subsequently act as nucleation seeds for hydrate growth. Moreover, they also demonstrated that the pore size has little effect on the formation mechanisms, aside from narrow pores where no carbon dioxide hydrate can be formed due to steric constraints. Bagherzadeh et al. 39,40 conducted MD simulations with rigid silica surfaces and demonstrated that confinement affects the interfacial curvature of the confined fluids as well as their distribution. Surprisingly, faster dissociation rates were observed in confined systems and for systems without a buffer water layer adjacent to the silica surface. The latter finding suggests that the surfaces destabilize the hydrate phase such that there exists a layer of disordered water close to the surfaces, and as such hydrate formation does not occur from the surfaces. Based on the same surface structure of silica, Liang et al. 41,42 identified that methane hydrate growth progresses preferentially toward the bulk solution and that the surfaces can act as a reservoir of methane molecules. In addition, they also noticed the formation of incomplete cages near the surfaces that could initiate methane hydrate nucleation. In an analysis of methane hydrate formation in clays of different pore sizes, Yan et al. 43 noted that the surface water exhibits certain cage-like structures. The presence of intercalated cations causes distortions and defects of the intercalated methane hydrates due to disruptions in hydrogen bonding. In a study comparing hydrophilic and hydrophobic surfaces, Nguyen et al. 44 found that the latter enhances local gas density and promotes the structural ordering of water near the solid surface, which consequently aids in gas hydrate formation. While most studies utilizing solid surfaces treated them as rigid and nonflexible structures, He et al. 45 explored the use of flexible structures of graphite and silica surfaces in order to study methane hydrate formation with various initial configurations of methane nanobubbles in slit-shaped nanopores. They demonstrated that a methane nanobubble of negative curvature (concave shaped) adsorbed on the hydrophobic graphite surface resulted in a lower Young− Laplace pressure at the fluid interface, which effectively lowers the bulk methane concentration; hence, there was no observable methane hydrate formation. On the other hand, the hydroxylated silica surface promotes the formation of a convex-shaped methane nanobubble, which consequently results in methane hydrate formation. Furthermore, they also showed that, in a homogenized solution, methane hydrate preferentially forms in the bulk away from the surfaces regardless of the hydrophobicity of the surfaces. Overall, these studies highlight the complexity of studying gas hydrate formation in porous media.
Recent experimental studies by Silvestre-Albero and coworkers suggested that porous materials can promote methane hydrate formation at conditions milder than those required in the bulk, thanks to the confinement effects. 46−50 They found that wet activated carbon, which contains well-defined slit pores, modifies the methane adsorption isotherms such that methane adsorption is greatly enhanced compared to that of the dry samples. Through synchrotron X-ray powder diffraction, they further confirmed that methane hydrate formation took place in the confined nanospace and that the storage capacity increased with increasing pore sizes up to an optimum size of 25 nm.
As mentioned above, there have been a number of studies that investigated methane hydrate growth in confinement by employing molecular simulations. However, these studies were carried out at temperatures and pressures where hydrate growth is certain in the bulk, i.e., at either very high pressures or very low temperatures. 38,[41][42][43]51 Therefore, it has so far not been possible to confirm whether the confinement effects can facilitate methane hydrate growth at milder conditions compared to the bulk. For instance, the work of He et al. 45 shows that while methane hydrate growth takes place in hydroxylated silica slit pores, no growth is observed in graphene slit pores. Although their work provides valuable The Journal of Physical Chemistry C pubs.acs.org/JPCC Article insights about the effects of surface chemistry, it cannot show whether confinement promotes methane hydrate growth because their simulations were carried out at 400 bar, at which point methane hydrate growth occurs in the bulk phase, i.e., without the effect of confinement.
In order to answer whether the confinement effects can facilitate methane hydrate growth at milder conditions compared to the bulk, we carried out a systematic study by performing molecular dynamics simulations of methane water systems in the bulk and in hydroxylated silica pores at various pressures. Our study shows that, indeed, methane growth can take place in narrow pores at pressures lower than those required for the growth of methane hydrate in the bulk.

■ METHODOLOGY
The direct phase coexistence methodology developed by Ladd and Woodcock 52 has been shown to predict the phase equilibria of gas hydrate systems in the bulk successfully. 53−58 In this method, a three-phase configuration consisting of a hydrate crystal, liquid water, and vapor slabs coexists. We used this method in order to study the gas hydrate systems in confinement.
Water, Methane, and Hydroxylated Silica Interaction Potentials. Water and methane molecules were represented by the TIP4P/ice water model, which is a four-site model with a massless particle for the oxygen charge, 59 and the TraPPE United Atom (TraPPE-UA) model, 60 respectively. 53,55,56 The hydroxylated silica surfaces were described by a CHARMMbased force field developed by Lopes et al., 61 which allows for the flexibility of the silica surfaces.
Initial Configuration of Bulk Systems. The preparation of the initial configuration of the bulk systems closely follows the procedure set out by Conde and Vega 53 and can be briefly summarized as follows. In order to generate the three-phase configuration, we first obtained a single-crystal structure of methane hydrate, which was determined from XRD measurements as reported by Kirchner et al. 62 While retaining the positions of carbon and oxygen atoms at their crystallographically determined positions, we added hydrogen atoms as well as dummy atoms to the oxygen atoms according to the geometry of the TIP4P/ice model. Hydrogen atoms were not added to the carbons because we used a united atom methane model. The resulting unit cell structure was then replicated in order to create a supercell containing 3 × 3 × 3 (27) sI methane hydrate unit cell. Finally, three slabs consisting of (i) 3 × 3 × 3 (27) sI methane hydrate unit cells, (ii) 1242 water molecules, and (iii) 216 methane molecules were equilibrated separately using short NVT-MD simulations. The y-and zdimensions of the water and methane slabs were kept equal to the hydrate slab (i.e., 3.563 nm × 3.563 nm). Then, they were brought together with a 0.1 nm buffer distance between them. Finally, energy minimization via steepest descent was carried out in order to avoid bad contacts. The dimensions of the initial configuration of the bulk systems ( Figure 1) were approximately 8.409 × 3.563 × 3.563 nm 3 .
Initial Configuration of Confined Systems. The initial configuration of the confined systems used in this study consisted of a three-phase configuration (i.e., solid hydrate crystal, liquid water, and methane gas) placed within a slitshaped silica nanopore as shown in Figure 2. In order to generate the silica slit pore, a silica block was first created by replicating the α-cristobalite unit cell 21 × 12 × 10 times. Then, a portion of the supercell that corresponded to the thickness of four sI methane hydrate unit cells (≈4.75 nm) was removed from the center of the silica block, leaving behind a slit pore. Finally, non-bridging oxygen atoms were subsequently protonated, yielding silanol groups on the surfaces.
The three-phase configuration in the slit pore was created from 7 × 5 × 4 (140) sI methane hydrate unit cells, which consisted of 1225 methane molecules and 6650 water molecules. In the crystalline form, this would occupy about 80% of the slit-pore volume shown in Figure 2. For the crystal slab, 3 × 5 × 4 (60) sI methane hydrate unit cells were used. The 700 methane molecules and 3845 water molecules in the remaining 4 × 5 × 4 (80) unit cells were used in order to create methane and water slabs, respectively, by equilibrating them separately using short NVT-MD simulations. The three slabs were then brought together with a 0.1 nm buffer distance between them in order to create a three-phase configuration. Finally, the three-phase configuration was placed in the slitshaped pore of the hydroxylated silica. In order to avoid bad contacts, energy minimization via steepest descent was performed. The dimensions of the initial configuration of the confined systems were 10.45 × 5.97 × 6.95 nm 3 . A flowchart illustrating the preparation of the initial configurations is provided in Figure S1.
Simulation Settings. All MD simulations were performed using GROMACS 2019.4 molecular dynamics simulations software. 63 The leapfrog algorithm with a time step of 2 fs was used for integrating Newton's equations of motion. Water molecules were held rigid by the SETTLE algorithm. 64 Periodic boundary conditions were applied in all directions. The short-range non-bonded interactions were modeled with the Lennard-Jones (LJ) potential. The Lorentz−Berthelot combining rules were used to calculate the cross-interaction LJ parameters of unlike atoms. The long-range Coulumbic interactions were calculated using a smooth particle mesh ewald (PME) method of a fourth order polynomial with a mesh width of 0.12 nm. 65 A cutoff distance of 11 Å was used   ■ RESULTS AND DISCUSSION Model Validation. Before investigating whether confinement effects promote methane hydrate growth, we tested the molecular models employed for methane and water for the bulk systems. For this, we conducted MD simulations in the NPT ensemble for the bulk systems at the following five pressures: 1, 10, 24, 40, and 100 bar. Previous studies have associated the decrease and increase in potential energy of the system with phase transitions, i.e., hydrate growth and dissociation, respectively. 27,53,55 Figure 3 shows the potential energy of the bulk systems as a function of the simulation time. For the systems simulated at 1, 10, and 24 bar, i.e., below the equilibrium pressure at 2°C, a characteristic spike signifying hydrate dissociation 53−55 followed by an immediate dip in the potential energy was observed, after which the simulations became unstable. In these simulations, a rapid expansion of the simulation box was observed as a result of the dissociation of the full hydrate slab and the consequent release of the enclathrated methane molecules to the gas phase. The hydrate dissociation rate was fastest at the lowest pressure. On the other hand, for the bulk system at 40 bar, the potential energy increased slightly but then stabilized after a microsecond, i.e., fluctuated around a mean value. The initial increase in the potential is associated with the partial dissociation of the hydrate phase, leading to the expansion of the simulation box. However, the subsequent stabilization of the potential energy shows that, although hydrate growth was not observed within the 6 μs simulations time, a state of three-phase coexistence was established ( Figure  S2) and growth may be observed at longer simulation times. Finally, in the case of 100 bar, the steady decrease in the potential energy over the course of the simulation indicated hydrate growth, which can be corroborated by the snapshots presented in Figure 4 where shrinkage of the simulation box can be noted. The characteristic two-step decrease in the potential energy, where a "steady" decrease in the potential energy is followed by a "rapid" decrease in the potential energy that was reported in previous studies for hydrate growth, 27,53,55 was not observed for the bulk system we simulated at 100 bar. This implies that the hydrate growth was "gradual". Overall, the results from the bulk systems demonstrate that the molecular models employed in them were able to correctly predict the behavior of the methane hydrates in bulk such that no growth was observed when the external pressure applied was below the bulk equilibrium pressure of 32 bar; above that pressure, either the growth or the onset of growth of the methane hydrates was observed.
We also looked into how dispersion corrections affect the simulation of methane hydrate growth in the bulk. For this, the simulation in the bulk at 24 bar was repeated with the dispersion corrections applied to both the energy and the pressure. In stark contrast to the experimental data and the simulation at 24 bar, in which dispersion corrections were not applied, turning the dispersion corrections on led, incorrectly, to complete crystallization within about 2.5 μs (Figures S2 and  S3). It is apparent that the dispersion corrections can yield artificial effects such that they result in an increased driving force for methane hydrate growth, which can lead to crystallization at conditions where such a result cannot be observed in the bulk. This is in line with the findings of Michalis et al., which highlighted that dispersion corrections will cause gas compression and that such effects are more pronounced at low pressures.
Effect of Confinement. In order to understand the effect of confinement on the growth of methane hydrates, MD simulations were conducted in the NPT ensemble where threephase configuration methane hydrates were subjected to the presence of a silica slit pore (Figure 2) at the following five pressures: 1, 10, 24, 40, and 100 bar. It is apparent from the results shown in Figure 5 that hydrate growth was observed in all simulations. The higher the pressure, the faster the growth rate of the methane hydrates became in confinement. The width of the silica pores, i.e., the distance between the two surfaces in the z-direction, decreased with increasing pressure ( Figure S5). This means that the confinement effect is the smallest in the 1 bar case; however, methane hydrate growth was still observed. Interestingly, only the simulation at 100 bar displayed the commonly reported characteristic of a two-step decrease in the potential energy. The results shown in Figure 5 can be used to interpret the progression of the simulation, i.e., a "gradual" growth step occurs in the first 1.75 μs, and a "spontaneous" growth step follows consecutively. On the other hand, all other simulations exhibit only "gradual" growth. Besides that, the growth rate was observed to be slower at 40 bar than that at lower pressures. This result is somewhat counterintuitive and may be explained by the fact that hydrate Figure 3. Variation of the potential energy of the bulk systems as a function of the NPT-MD simulation time. At 1, 10, and 24 bar, the methane hydrates are not stable and dissociate, corresponding to an increase in the potential energy. On the other hand, at 40 bar a threephase coexistence is observed, corresponding to a plateau in the potential energy. Finally, at 100 bar methane hydrate growth is observed, corresponding to a decrease in the potential energy.
The Journal of Physical Chemistry C pubs.acs.org/JPCC Article formation is stochastic, meaning that such an effect is more pronounced when the driving force is low.
In the snapshots displayed in Figure 6 for the simulation conducted at 100 bar, at first a "gradual" growth of methane hydrates from the hydrate−water interface ( Figure 6A) and the formation of a convex-shaped methane slab can be observed. Subsequently, a "spontaneous" lateral growth of the hydrate slab can be seen due to the formation of a methane nanobubble ( Figure 6B−D), which was also observed by He et al. 45 In addition, incomplete and distorted hydrate cages can be observed near the silica surfaces. Although complete cages were initially present adjacent to the silica surfaces, these were replaced by a thin water layer during the course of the simulation without causing the dissociation of the remaining cages. This also aligns with earlier experimental and theoretical observations, 39,40,45,68 which showed that the surfaces are coated with a thin layer of ordered water due to strong hydrogen bonding with the silanol groups on silica surface. It has been hypothesized that the stabilizing effect of silica on the cage-like water structures may promote hydrate formation 45 and that the formation of cage-like water structures is a necessary step for gas hydrate formation. 69 In the "gradual" growth step, hydrate growth is a mass transfer-limited process, as the methane molecules must first dissolve and diffuse through the water slab to the hydrate− water interface in order to be enclathrated; hence, the potential energy decreases steadily. However, in the "spontaneous" growth step, when the growing hydrate slab came in contact with the supersaturated methane solution, the growth was accelerated and consequently accompanied by a steep decrease in the potential energy of the system. 55 Prior studies have noted the effect of surfaces with different hydrophobicities on the wetting angle of the fluid coming in contact with the surfaces 45,70,71 such that it affects the Young−Laplace pressure (P YL ) that arises from the interfacial curvature. Hence, this influences the concentration of methane in the aqueous phase, which plays an important role in hydrate formation. 24,72−74 Consistent with the findings of He et al., 45 our results show that hydrophilic hydroxylated silica surfaces promote convexshaped (positive curvature) methane nanobubbles where the direction of P YL is pointing outward. In addition, they also pointed out that P YL pointing outward increases the aqueous concentration of methane and consequently a higher hydrate growth rate can be observed.
Comparing the energy profiles of the bulk ( Figure 3) and confined ( Figure 5) systems, it is clear that in the bulk systems the growth of the methane hydrate crystals was dependent on the external pressure applied and that the simulation outcome reproduced the experimental methane−water phase equilibria. On the other hand, in the confined systems, hydrate growth can be observed regardless of the external pressure applied, even at pressures lower than the minimum pressure required for the formation of methane hydrates in the bulk. This result could be attributed to an increased tangential pressure as well as an increased probability of formation for the clathrate-like water structures in confinement. 75 The latter is essential for the spontaneous onset of hydrate formation according to the cage adsorption hypothesis by Guo et al. 69

■ CONCLUSIONS
In this study, by performing MD simulations, we investigated whether methane hydrate growth can take place in confinement at milder conditions compared to the bulk. This was largely motivated by recent experimental work that used porous materials for the synthetic growth of methane hydrates. Furthermore, previous molecular simulation studies of methane hydrate growth in pore systems were carried out at very high pressures, making it impossible to determine whether confinement could enable methane hydrate growth at pressures lower than those required in the bulk. Prior to this study, little was known about hydrate growth in confinement at low simulation pressures (≤100 bar). In order to gain an understanding of the effect of confinement on methane hydrate formation at low pressures, we employed the direct phase coexistence approach in our MD simulations. In the simulated bulk systems, the external pressure dictates whether methane hydrate growth was observed at constant temperature such that the behavior of the methane hydrates obeys the expected experimental methane hydrate phase equilibria. On the other hand, in confinement, methane hydrate growth was observed even at the pressures lower than the minimum pressure required for hydrate formation in the bulk regardless of the external pressure applied, suggesting, indeed, that confinement effects enable methane hydrate growth at relatively milder conditions. We further showed that applying dispersion corrections (i.e., tail corrections) could yield artificial effects that incorrectly lead to the crystallization of the hydrates. Overall, our findings provide insights into how confinement promotes hydrate formation and can help others develop accurate methodologies and simulation settings in order to study the behavior of gas hydrates in confinement. However, many points remain to be investigated. For instance, the pore width and shape can also be factors for the hydrate growth. Furthermore, it is not totally unambiguous how the external pressure applied in the simulations of the confined systems is representative of the pressure of the bulk methane used in the experiments in order to make synthetic methane hydrates in wet porous materials. In order to develop a full picture of the effects of confinement on hydrate growth in porous materials, additional modeling and experimental research will be needed.
Main procedures for MD simulations of the confined systems, snapshots from the bulk system simulation at 40 bar, the variation of the potential energy of the bulk systems as a function of the NPT-MD simulation time with and without dispersion corrections at 24 bar and 271.65 K, snapshots from the bulk system simulation with dispersion corrections applied at 24 bar, and the variation of the pore widths of the confined systems as a function of the NPT-MD simulation time at various pressures and 271.65 K (PDF) Files containing simulation input parameters (ZIP)

Notes
The authors declare no competing financial interest.

■ ACKNOWLEDGMENTS
We are grateful to the UK Materials and Molecular Modelling Hub for computational resources, which is partially funded by EPSRC (EP/P020194/1). K.B.Y. is grateful for the H. Walter Stern Ph.D. Studentship from the UCL Department of Chemical Engineering.