Competition between CH₄ hydrate formation and phase separation in a wetted metal–organic framework MIL-101 at moderate subcooling: molecular insights into CH₄ storage

Abstract

Adsorption-hydration hybrid technology has emerged as a promising technology to store CH₄ in porous materials, as it synergistically improves CH₄ storage capacity by combining CH₄ adsorption and hydrate formation. However, the fundamental mechanism involved in this technology remains elusive. Herein, we perform systematic molecular dynamics simulations to explore CH₄ hydrate formation in a metal–organic framework MIL-101 at moderate subcooling. Simulation results reveal that at moderate subcooling, CH₄ hydrate formation and phase separation of CH₄ to form nanobubbles occur simultaneously, and these two processes compete with each other for CH₄ molecules in the solution. The outcome of the competition is primarily governed by the relative stabilities of CH₄ hydrate solids and CH₄ nanobubbles, which are closely related to their sizes. It is revealed that CH₄ hydrate formation occurs exclusively in the outer space of MIL-101 cavities, whereas phase separation of CH₄ to form nanobubbles takes place in the MIL-101 cavities and their outer space simultaneously. The small nanobubbles in the MIL-101 cavities gradually shrink and finally disappear, as CH₄ molecules therein diffuse out and grow into large nanobubbles and large hydrate solids in the outer space. Moderate subcooling appears to facilitate the formation of large ordered CH₄ hydrate solids containing sI and sII domains. Additionally, it is found that lower subcooling and presence of MIL-101 both promote phase separation of CH₄. In the evolution of large and small nanobubbles during phase separation, coalescence of CH₄ nanobubbles and an interesting phenomenon similar to Ostwald ripening are observed. The molecular insights into the effects of the degree of subcooling on CH₄ hydrate formation in MIL-101 provide bottom-up guidance on optimizing pressure-temperature conditions for CH₄ storage in porous materials with adsorption-hydrate hybrid technology.

---

1. Introduction

With the rapid growth of the global economy and population, the demand for energy resources has increased consistently. Currently, natural gas (mainly CH₄) is a heavily consumed energy resource and could meet this demand due to its abundance worldwide, high energy density, cleanliness and low cost. Nevertheless, efficient methods for storing and transporting natural gas are required to maximize the usage of this clean energy resource. One promising method is utilizing gas physisorption in porous materials. With readily tunable structures and pore sizes, metal–organic frameworks (MOFs) are a very interesting class of porous materials for CH₄ storage. In the past decade, tremendous efforts have been made to improve CH₄ storage capacity in MOFs by optimizing their chemistry, structure and functionality. Recent experimental studies reported that pre-humidified MOFs could enhance CH₄ storage capacity,^1,2 especially at high pressure and moderate temperature.³ This interesting phenomenon is primarily attributed to the formation of CH₄ hydrate in the presence of H₂O in MOFs. Gas hydrates are one kind of ice-like solids with gas molecules entrapped into hydrogen-bonded water polyhedral cages. Similar to porous materials, CH₄ hydrate itself has excellent storage capacity, i.e., one volume hydrate solid can contain approximately 170 volumes of CH₄ at the standard temperature and pressure.⁴ Moreover, the unique self-preservation anomaly of gas hydrate⁵ can be leveraged for long-time preservation and transport of CH₄ hydrate in MOFs at ambient pressure and low temperatures. Therefore, researchers have proposed to combine adsorption and hydrate formation to efficiently store CH₄ in MOFs, i.e., the so-called adsorption-hydrate hybrid method.

Toward this end, there have been increasing experimental and theoretical investigations to explore CH₄ storage in MOFs via the adsorption-hydrate hybrid method.^6,7 Liu and co-workers explored CH₄ adsorption and hydrate formation in wet ZIF-8 with different water contents and observed that the storage capacity was raised by 56% with a water content of 35.1 wt%, which was conjectured to arise from hydrate formation in wet ZIF-8 pores.³ The experimental work by Casco et al. directly identified CH₄ hydrate formation in pre-humidified MIL-100(Fe) and ZIF-8, and revealed that the surface chemistry, pore structure and water–framework interactions would significantly affect the amount and nature of gas hydrates formed in MOFs.⁸ Nevertheless, it was speculated that CH₄ hydrate formation occurred primarily on the external surface or in the interparticle space due to the small cavities in most MOFs^9,10e.g., the cavity sizes in ZIF-8, MIL-53 and HKUST-1 < 1 nm. A recent experimental study on CH₄ hydrate formation in Cr-soc-MOF-1 and Y-shp-MOF-5 indicated that forming hydrate nanocrystals in MOF cavities would require a cavity size larger than the unit cell size of CH₄ hydrate (1.2 nm); moreover, encapsulating CH₄ hydrate nanocrystals inside Cr-soc-MOF-1 cavities could increase gas storage capacity by 50%.¹¹ Chen and co-workers investigated the morphology and growth kinetics of CH₄ hydrate in wetted ZIF-8 particles and quantified the effects of particle size, water content and pressure.^12,13 Additionally, the phase behavior of gas hydrate has been reported to change in porous materials. Liu et al. revealed a pore-diameter dependent effect on CH₄ hydrate equilibrium in wetted porous materials (ZIF-8, HKUST-1 and zeolite 13X) with a pore size ranging from 0.35 to 100 nm, i.e., strong inhibitory effect in nano-pores and marginal effect in macro-pores.¹⁴ The thermodynamic behavior of gas hydrates was observed to be relatively inhibited in the intergranular meso/macro pores of MIL-53.¹⁵ Certain MOFs have been reported to serve as effective promoters for CH₄ hydrate formation for enhancing natural gas storage. Koh and co-workers drastically improved CH₄ storage by over tenfold upon adding HKUST-1, ZIF-8 and ZIF-67, and these MOFs were found to significantly increase water-to-hydrate conversion and remarkably reduce the induction time for hydrate nucleation.^9,10,16,17 Recent theoretical studies have been mainly focused on the formation kinetics of CH₄ hydrate in MOFs. Our previous molecular dynamics (MD) simulation on CH₄ storage in mesoporous MIL-101 revealed that though the pre-adsorbed water in MIL-101 cavities substantially reduces CH₄ adsorption capacity, considerable CH₄ hydrate preferentially forms in the intergranular space and significantly enhances the overall storage capacity.¹⁸ A two-way transport of CH₄ and water in ZIF-8 was reported in MD simulation studies, and the synergist effect of physisorption and intergranular hydrates was found to promote CH₄ storage in ZIF-8.^19–21 Overall, the phase behavior and formation process of CH₄ hydrate in MOFs are far more complex than in a bulk phase.

A fundamental understanding of CH₄ hydrate formation in porous material is indispensable to manipulating natural gas hydrate, which exists abundantly in deep-sea porous sediments and is considered a new energy resource.^22,23 Conceptually, CH₄ hydrate formation is inhibited in marine sediments due to the reduced water activity and capillary effect of porous sediments.²⁴ Several MD simulation studies explored the nucleation and growth of CH₄ hydrate in the slit-nanopores of silica or clay to reveal the effects of geological factors on hydrate formation in marine sediments, including sediment solid surfaces, seawater ions, and organic matter.^25–28 Recent experimental studies inferred that the crystal structure and stability of CH₄ hydrate confined in sediments might be mainly affected by the hydrophilic/hydrophobic properties of sediment surfaces.²⁹

Despite many achievements, the microscopic mechanism of hydrate formation in porous materials is still not fully elucidated. Low temperature near ice-points is often adopted in experimental and theoretical studies on CH₄ storage in MOFs using adsorption-hydrate hybrid technology.^3,8,11,18,19 Furthermore, high subcooling is commonly implemented in MD simulation studies on gas hydrate formation in a bulk phase.^30,31 Subcooling is referred to as the difference between the melting temperature of gas hydrate and an applied temperature. It is worth mentioning that CH₄ hydrate occurs at a moderate temperature in hydrate reservoirs due to the presence of a geothermal gradient therein. To the best of our knowledge, how moderate subcooling affects CH₄ adsorption and hydrate formation in porous material has not been explored thus far. To address this issue, MD simulations have been conducted in this study to investigate CH₄ hydrate formation at moderate subcooling in a mesoporous MOF namely MIL-101. The molecular insights obtained could provide a fundamental understanding of the formation process of natural gas hydrate in marine sediments, and bottom-up guidance on optimizing pressure-temperature conditions for CH₄ storage in porous materials with adsorption-hydrate hybrid technology.

2. Simulation models and methods

In this study, MIL-101 was adopted as a porous material for CH₄ hydrate formation due to its large pore volume. MIL-101 has a three-dimensional zeotype cubic structure with a lattice constant of 8.9 nm.³² Two kinds of spherical-like mesoporous cavities exist in MIL-101: a large cavity and a small cavity with a free diameter of 3.4 nm and 2.9 nm, respectively, which are sufficiently large to accommodate a unit cell of CH₄ hydrate. MIL-101 is assembled from two building blocks of 1,4-benzenedicarboxylic acid (1,4-BDC) and the octahedral Cr₃O trimer, which form corner-sharing supertetrahedra. Such a superhedron has 4 vertices and 6 edges and contains a free aperture of 0.86 nm, which is large enough for water and CH₄ molecules to cross. The large cavity contains 28 supertetrahedra accessible through hexagonal and pentagonal windows with 1.6 × 1.5 nm² aperture, while the small cavity consists of 20 supertetrahedra accessible through pentagonal windows with 1.2 × 1.2 nm² aperture.

The structure of MIL-101 was constructed based on experimental crystallographic data, and then a large cavity and a small cavity were cleaved with the cleaved bonds of Cr₃O trimers terminated by methyl groups, following the methods in our previous study.³³ The large and small cavities were placed into a simulation box with a size of 16.3 × 7.4 × 7.3 nm³, and the distance between the surfaces of the two cavity structures was 4.4 nm (Fig. 1). We used such a simplified MIL-101 model for two considerations. First, MOF samples in experimental studies often include intergranular mesopores among MOF particles, in addition to intrinsic cavities like the large and small MIL-101 cavities, thus, the outer space of the cavities can be regarded as intergranular mesopores. Second, the unit cell of MIL-101 is as large as 8.9 nm and direct simulation using a unit cell would be time-consuming. Thereafter, 2263 CH₄ and 26 014 H₂O molecules were randomly placed in the simulation box to form a homogeneous mixture with CH₄ mole fraction x_CH4 = 0.08. Such a homogeneous aqueous solution of CH₄ containing MIL-101 cavities was taken as the initial configuration for simulation up to 900 ns to explore CH₄ hydrate formation. It should be noted that the CH₄ solution was moderately supersaturated with x_CH4 = 0.08, much lower than that in sI crystalline hydrate (x_CH4 = 0.148).

Fig. 1 Starting configuration of the simulation system. A large cavity and a small cavity of MIL-101 separated by 4.4 nm were placed in a box filled with a homogeneous aqueous solution of CH₄. The cavities of MIL-101, CH₄ and water are shown as magenta sticks, green balls and light blue lines, respectively.

H₂O and CH₄ molecules were represented by TIP4P/Ice³⁴ and OPLS-UA models,³⁵ respectively. The Lennard–Jones parameters for MIL-101 were taken from the Universal Force Field (UFF),³⁶ which can precisely describe the diffusion and adsorption behavior of gas molecules in MOFs. The cross-interaction parameters between unlike species were calculated using the Lorentz–Berthelot combination rule. The initial configurations of large and small cavities of MIL-101 were structurally optimized using density-functional theory (DFT) in DMol3, and the atomic charges of MIL-101 were taken from our previous study.¹⁸ Electrostatic interactions were calculated using the particle mesh Ewald method.³⁷ Short-range nonbonded interactions were estimated with a cutoff of 1.0 nm. The pressure and temperature for the system were controlled with the Parrinello-Rahman barostat³⁸ and the Nosé–Hoover thermostat,³⁹ respectively. Trajectories were integrated using the leapfrog method with a time step of 2 fs and coordinates were saved every 40 ps. The MIL-101 cavities were position restrained with a moderate force constant of 1000 kJ mol⁻¹ nm⁻² during the simulation, thus, the cavities showed partial flexibility. The temperature and pressure for the system were set to 285 K and 50 MPa, corresponding to moderate subcooling of about 19 K. This is different from the high subcooling of 54 K (250 K and 50 MPa) usually adopted in previous MD simulation studies on CH₄ hydrate formation from a gas/water two-phase system or a homogeneous CH₄ solution system.^30,31,40 The choice of moderate subcooling is based on the fact that there is relatively low subcooling in actual CH₄ hydrate reservoirs under the sea floor, e.g., 283.15 K and 15 MPa for hydrate reservoirs in the South Sea of China.⁴¹ In addition, the conditions of 285 K and 15 MPa were also simulated to further examine the effect of lower subcooling of about 7 K on CH₄ hydrate formation. To reveal the role of MIL-101 in CH₄ hydrate formation and evolution, a control system without MIL-101 was simulated. Three repeated runs (Run1, Run2 and Run3) were performed using GROMACS-5.0.7 (ref. 42) for each condition, i.e., at 285 K and 50 MPa in the presence and absence of MIL-101, and at 285 K and 15 MPa in the presence of MIL-101. The corresponding subcooling temperatures were estimated based on the three-phase (water-hydrate-methane) coexistence curves of CH₄ hydrate,⁴³ which was simulated with the same potential models as used in this work (see Fig. S1†).

3. Results and discussion

3.1 Competition between CH₄ hydrate formation and phase separation of CH₄ in solution at moderate subcooling

To monitor the formation of CH₄ hydrate in the system at moderate subcooling of 19 K (at 285 K and 50 MPa), a four-body structural order parameter (F₄) for water molecules was analyzed. F₄ is defined as ensemble averaged <cos(3θ)> and describes the torsion angle θ between the oxygen atoms of two water molecules within 0.35 nm and the outermost hydrogen atoms. It is often adopted to characterize the phase change of water by tracking the evolution of average value: 0.7, −0.04 and −0.4 for hydrate, liquid water and ice, respectively.⁴⁴ The evolution of F₄ during three repeated runs is shown in Fig. 2A. Additionally, the algorithm proposed by Jacobson et al. was used to recognize cage geometries⁴⁵ and seven types of cages (5¹², 5¹²6², 5¹²6³, 5¹²6⁴, 4¹5¹⁰6², 4¹5¹⁰6³ and 4¹5¹⁰6⁴) were identified to be the most abundantly formed during CH₄ hydrate formation. Fig. 2B displays the number of water molecules N_w in hydrate form (i.e., forming these seven types of cages).

Fig. 2 Evolution of (A) F₄ order parameter and (B) number of water N_w in hydrate form during three repeated runs at 285 K and 50 MPa with subcooling of 19 K.

As shown in Fig. 2A, F₄ values in the three repeated runs increase from −0.04 at the beginning of the simulations, indicating that CH₄ hydrate can easily nucleate in the supersaturated CH₄ solution with x_CH4 = 0.08. Afterwards, F₄ values in Run1 and Run2 increase sharply to 0.16 and 0.05 within the initial 200 ns and 100 ns, and then increase gradually to 0.24 and 0.11 in the remaining simulation, respectively. This suggests two growth stages for CH₄ hydrate: initial rapid growth and subsequent slow growth. In contrast, a different phenomenon is observed in Run3. The F₄ value in Run3 increases quickly to a maximum of 0.05 at 75 ns and then drops slowly to −0.04 at 310 ns (Fig. 2A), indicating that the formed CH₄ hydrate in Run3 decomposes gradually after the rapid hydrate growth stage. Furthermore, during the initial 50 ns, the F₄ value in Run3 increases faster and is larger than those in Run1 and Run2, reflecting faster hydrate growth at the initial stage in Run3. The final F₄ value in Run1 at 900 ns is much larger than that in Run2 (Fig. 2A) and more water molecules in Run1 are converted into hydrates (Fig. 2B), i.e., more hydrate solids are formed in Run1. It is noted that the final F₄ values in Run1 and Run2 (0.24 and 0.11) are much smaller than that for hydrate (0.7), as the water in the system is excessive (the composition of CH₄ in the system x_CH4 = 0.08 is lower than x_CH4 = 0.148 in sI crystalline hydrate) and there are still many water molecules in a liquid state. Such different formation behavior of CH₄ hydrate in the three repeated runs is also reflected in the evolution of N_w in hydrate form (Fig. 2B). In Run1 and Run2, N_w increases quickly to 7500 and 3100 within 200 ns and 100 ns, and then gradually to 10 500 and 5500 in the remaining simulation. In contrast, N_w in Run3 increases to 2500 at 75 ns and then decreases slowly to 0 at 310 ns. Further analysis will be performed below to elucidate the underlying mechanism in the formation kinetics of CH₄ hydrate.

In addition to hydrate formation, phase separation of CH₄ is also observed to occur in the supersaturated CH₄ solution with initial x_CH4 = 0.08 at moderate subcooling of 19 K (at 285 K and 50 MPa). It is noted that phase separation of CH₄ reduces gas concentration in the solution, which has been demonstrated to play an important role in hydrate formation as hydrate grows faster at a higher gas concentration.³¹ Thus, the CH₄ mole fraction in water (x_CH4) and the number of CH₄ molecules in nanobubbles (N_CH4) in the small and large cavities of MIL-101 and their outer space were analyzed. As shown in Fig. 3D, in Run1, N_CH4 in the outer space and in the large cavity of MIL-101 increase to a maximum of 230 and 260 at 50 ns and 120 ns, respectively, indicating the phase separation of CH₄ in these regions. As a consequence, x_CH4 in the outer space and in the large cavity decreases to 0.07 and 0.005 (Fig. 3A). In the subsequent simulation, with the sustained growth of CH₄ hydrate, x_CH4 in the outer space gradually increases to 0.085 and N_CH4 in nanobubbles decreases to 0, i.e., the nanobubbles in the outer space completely decompose. In Run2, similar phase separation of CH₄ in solution occurs. As shown in Fig. 3E, large nanobubbles containing up to 1200 CH₄ molecules form in the outer space, while small nanobubbles with 120 CH₄ molecules form in the large cavity of MIL-101. With the growth of CH₄ hydrate, the small nanobubbles in the large cavity quickly decompose and the large nanobubbles gradually shrink to smaller ones with 850 CH₄ molecules at 900 ns. Accordingly, x_CH4 in the outer space decreases sharply to 0.04 and then increases slowly to 0.052 (Fig. 3B). In Run3, phase separation of CH₄ occurs in the outer space of MIL-101 cavities and no obvious phase separation is observed in MIL-101 cavities. N_CH4 in nanobubbles in the large and small cavities is close to 0, while most of the CH₄ molecules (up to 2000) in the outer space gradually convert into nanobubbles (Fig. 3F). Thus, x_CH4 in the outer space decreases continuously to 0.01 during the whole simulation (Fig. 3C).

Fig. 3 Evolution of (A–C) CH₄ mole fraction in water (x_CH4) and (D–F) a number of CH₄ molecules in nanobubbles (N_CH4) in the absence of MIL-101 and in the presence of MIL-101 (in the small and large cavities and in their outer space) during three repeated runs at 285 K and 50 MPa with subcooling of 19 K.

It is observed that phase separation of CH₄ in solution is a process of attracting and consuming aqueous CH₄ molecules. During phase separation of CH₄ in the large cavity of MIL-101 in Run1 (0–120 ns) and Run2 (0–80 ns), aqueous CH₄ molecules therein are gradually depleted due to the formation of nanobubbles (Fig. 3D and E), and x_CH4 in the large cavity is much lower than x_CH4 in the outer space (Fig. 3A and B). Thus, a gas concentration gradient is established between these two regions and drives CH₄ molecules to diffuse into the large cavity to further promote phase separation therein, as evidenced by the increase in the number of CH₄ in the large cavity during 0–120 ns in Run1 and during 0–80 ns in Run2 (Fig. 4). By contrast, such CH₄ transport from the outer space into MIL-101 cavities is not observed in the large cavity in Run3 and in the small cavity during the three repeated runs, as no obvious phase separation occurs therein. Interestingly, as shown in Fig. 4, the number of CH₄ molecules in the large cavity decreases from a maximum during 120–900 ns in Run1 and during 80–300 ns in Run2, while the number of CH₄ molecules in the large cavity in Run3 and in the small cavity during the three repeated runs decreases constantly from the beginning of the simulation. This indicates that CH₄ molecules diffuse out from MIL-101 cavities into the outer space, probably to facilitate the formation of large nanobubbles and large hydrates in the outer space of cavities, as discussed below. Such two-way migration of CH₄ molecules between the cavities and outer space was also observed in a recent experimental study on CH₄ hydrate formation in activated carbons.⁴⁶ As a result, in the late stage of simulation, x_CH4 in the cavities of MIL-101 is low and no obvious nanobubble forms therein (Fig. 3A–C). In contrast, the finally formed large CH₄ nanobubbles are always located in the outer space of MIL-101 cavities.

Fig. 4 Evolution of the number of CH₄ molecules in the large and small cavities of MIL-101 during three repeated runs at 285 K and 50 MPa with subcooling of 19 K.

In order to directly inspect the processes of CH₄ hydrate formation and phase separation of CH₄ in the system at moderate subcooling of 19 K, we prepared the simulation snapshots of Run1, Run2 and Run3 (Fig. 5), which show the evolution of seven types of hydrate cages. For easy visualization, the periodic boundary of the simulation box was modified following the procedure in Fig. S2† to ensure the nucleation and growth of hydrate within the box (MIL-101 and CH₄ nanobubbles may cross the box boundary). In the homogeneous CH₄ solution in the presence of MIL-101 (Fig. 5A(I)–C(I)), hydrate cages uniformly form and generate several cage clusters (Fig. 5A(II)–C(II)). Most of these clusters grow to large hydrate (Fig. 5A(II)–(VII), B(II)–(VII), and C(II)–(IV)), while a few decompose (e.g., the two clusters on the left side in Run2, see Fig. 5B(II)–(IV)). As shown in Fig. 5A(II)–(VII), with the hydrate growth in Run1, CH₄ molecules are gradually converted into hydrate and finally depleted in the solution (Fig. 5A(VII)). This indicates that CH₄ hydrate formation is also a process of attracting and consuming many aqueous CH₄ molecules. From Fig. 5, it is found that hydrate formation always occurs in the outer space, not in the MIL-101 cavities due to the low x_CH4 therein (Fig. 3A–C). Under the moderate subcooling conditions, phase separation of CH₄ to form nanobubbles occurs in the supersaturated CH₄ solution and competes with the formation of CH₄ hydrate for aqueous CH₄ molecules. In Run1, phase separation of CH₄ takes place in the large cavity of MIL-101 and forms a small nanobubble, which decomposes completely with the growth of large hydrate in the outer space (Fig. 5A(III)–(VI)). In Run2, during the phase separation of CH₄, one small nanobubble and one large nanobubble form in the large cavity of MIL-101 and in the outer space (Fig. 5B(III)), respectively. Later, with the sustained growth of hydrate, the small nanobubble in the large cavity of MIL-101 completely disappears, while the large nanobubble in the outer space shrinks slowly (Fig. 5B(III)–(VII)). In Run3, one large nanobubble forms in the outer space of MIL-101 cavities during phase separation of CH₄, while two cage clusters grow to two middle-sized hydrates (Fig. 5C(II)–(IV)). Interestingly, with the sustained growth of the large CH₄ nanobubble, both middle-sized hydrates gradually decompose due to the accumulation of CH₄ molecules in the nanobubble (Fig. 5C(IV)–(VII)). Such competition between hydrate formation and CH₄ phase separation is also visualized in Videos S1–S3† for Run1, Run2 and Run3, respectively. The results for such competition can be summarized as follows: in Run1, almost all CH₄ molecules in the system are converted into two hydrates (Fig. 5A(VII)); in Run2, about half of CH₄ molecules in the system form a large nanobubble and the other CH₄ molecules form a large hydrate (Fig. 5B(VII)); in Run3, almost all CH₄ molecules in the system form a huge nanobubble (Fig. 5C(VII)). Additionally, the stability of CH₄ nanobubbles and CH₄ hydrate appears to be related to their sizes, i.e., the larger the size, the higher the stability. Thus, during the growth of large hydrate in Run2, the small nanobubble in the large cavity of MIL-101 decomposes more easily than the large nanobubble in the outer space (Fig. 5B(III)–(V)), whereas two middle-sized hydrates in Run3 (Fig. 5C(IV)) decompose more easily than the large hydrate in Run2 (Fig. 5B(IV)). Thus, the outcome of the competition between hydrate formation and phase separation is mainly determined by the relative stabilities of CH₄ hydrate solids and CH₄ nanobubbles.

Fig. 5 CH₄ hydrate formation and phase separation of CH₄ to form nanobubbles during Run1 (A(I)–(VII)), Run2 (B(I)–(VII)) and Run3 (C(I)–(VII)) at 285 K and 50 MPa with subcooling of 19 K. Hydrate cages are displayed as sticks with different colors (green for 5¹², blue for 5¹²6², red for 5¹²6³, orange for 5¹²6⁴, cyan for 4¹5¹⁰6², purple for 4¹5¹⁰6³ and pink for 4¹5¹⁰6⁴). MIL-101 cavities, CH₄ and water molecules are shown as magenta sticks, green balls and light blue lines, respectively. The black circles in the panels indicate the CH₄ nanobubbles formed in MIL-101 cavities.

3.2 Structures and crystallinity of incipient CH₄ hydrate formed at moderate subcooling

It is found that moderate subcooling may facilitate the formation of crystalline CH₄ hydrate. There are two types of crystalline forms for CH₄ hydrate, namely, cubic structure I (sI), and cubic structure II (sII). sI hydrate is constructed from 5¹²6² and 5¹² cages with a ratio of 6 : 2, while 5¹²6⁴ and 5¹² cages constitute sII hydrate with a ratio of 8 : 16, i.e., the unit cell of sII hydrate has eight times more 5¹² cages than in sI hydrate.⁴⁷ Moreover, 5¹² cages are arranged in different ways in sI and sII hydrates: in sI hydrate, 5¹² cages fill the interstices of 5¹²6² cages but never share faces, whereas many 5¹² cages share faces in sII hydrate. It is observed that a large ordered incipient CH₄ hydrate containing sI and sII domains forms on the right side of the system in Run1 (Fig. 5A) and in the middle of the system in Run2 (Fig. 5B), and the detailed formation process is displayed in Fig. 6A and B, respectively. A sII motif containing several 5¹² and 5¹²6⁴ cages forms in the outer space of MIL-101 cavities (Fig. 6A(I) and B(I)), which originated from an amorphous cluster of cages. This sII motif grows into a medium-sized sII domain, and a sI motif is observed to emerge on the lateral side of the sII domain (see the circles in Fig. 6A(II) and B(II)), which serves as a nucleus of the sI domain. Subsequently, the sII domain and the sI motif in the same hydrate co-grow simultaneously and rapidly (Fig. 6A(III) and B(III)). At the end of the simulation, a large hydrate with long-range order forms, which occupies almost the whole space between the large and small cavities of MIL-101. Moreover, crystalline sI and sII domains coexist in the ordered hydrate by connecting to each other with the linking cage 5¹²6³ (Fig. 6A(IV) and B(IV)). The dynamic formation process of the large ordered CH₄ hydrate in Run1 and Run2 is also visualized in Videos S4 and S5.† The co-formation of sI and sII CH₄ hydrate observed in this work is considered a kinetically favored process. It is well known that CH₄ hydrate adopts sI as a thermodynamically stable phase. Nevertheless, at the nucleation and early growth stages of the ordered CH₄ hydrate solids, the 5¹² cage is the most abundant and its growth is obviously earlier and faster than other cage types (Fig. 7C and D), and these abundantly formed face-sharing 5¹² cages (Fig. 6A(I) and B(I)) could facilitate the formation of sII hydrate. The reason why the 5¹² cage is kinetically favored over other cage types is because its size is most fitting for CH₄, and more importantly, the angle in its pentagonal ring is quite close to the angle in water hydrogen bonding. Interestingly, such coexistence of sI and sII was also reported in a previous experimental study on CH₄ hydrate formation at moderate subcooling.⁴⁸ It is noted that stochasticity plays a role in the formation of these ordered incipient hydrates. The incipient large hydrate formed on the left side of the system in Run1 (Fig. 5A(VII)) and the two middle-sized hydrates formed in Run3 (Fig. 5C(IV)) are relatively amorphous with short-range order, only containing certain small sI or sII motifs (Fig. S3†).

Fig. 6 Detailed formation process of a large ordered CH₄ hydrate containing sI and sII domains during Run1 (A(I)–(IV)) and Run2 (B(I)–(IV)) at 285 K and 50 MPa with subcooling of 19 K. Cages are in different colors (green for 5¹², blue for 5¹²6², red for 5¹²6³, and orange for 5¹²6⁴). MIL-101 cavities and CH₄ are shown as magenta sticks and green balls, respectively. The systems are rotated to clearly show the crystalline hydrate structure.

Fig. 7 (A and B) Crystallinity and (C and D) number of seven types of cages for the large ordered hydrate formed on the right side of the system during Run1 (see Fig. 5A) and in the middle of the system during Run2 (see Fig. 5B) at 285 K and 50 MPa with subcooling of 19 K.

To quantify the crystallinity degree of the large ordered CH₄ hydrates formed in Run1 and Run2, crystallinity was calculated from the ratio of the number of characteristic cage links of a crystalline hydrate to the total number of cage links in the system.³¹ Crystallinity with a value 0 indicates a completely amorphous hydrate and a value of 1 means a perfect crystalline hydrate. Fig. 7A and B display the crystallinities of the ordered hydrates formed on the right side in Run1 and in the middle in Run2, respectively. Fig. 7C and D illustrate the corresponding evolution of seven types of cages. For the hydrate in Run 1 at 900 ns, its sI and sII crystallinities are about 0.4 and 0.5 (Fig. 7A), suggesting that the fractions of sI and sII structures are comparable. At the end of simulation (900 ns), as shown in Fig. 7B, the sI and sII crystallinities of the ordered hydrate in Run2 are about 0.2 and 0.55, respectively, indicating that most structures of the hydrate belong to sII and a small portion belongs to sI. Additionally, it is found that the evolution of sI and sII crystallinities is closely related to the growth of sI and sII domains in the ordered hydrate, i.e., the evolution of seven types of cages therein. For example, the number of sI-specific 5¹²6² cages in the ordered hydrate in Run1 increases from 10 at 100 ns to 75 at 450 ns (Fig. 7C), indicating the growth of the sI domain, and makes the corresponding sI crystallinity in Run1 increase from 0.15 at 100 ns to 0.45 at 450 ns (Fig. 7A). Similarly, with increasing the number of sII-specific 5¹²6⁴ cages from 500 to 900 ns (Fig. 7C), the sII crystallinity in Run 1 is enhanced from 0.35 to 0.5 at 900 ns (Fig. 7A). Due to the coexistence of large sI and sII domains in the ordered hydrates in Run1 and Run2, 5¹², 5¹²6² and 5¹²6⁴ cages are the most abundant cage type followed by the 5¹²6³ linking cage (Fig. 7C and D). Several meta-stable 4¹5¹⁰6², 4¹5¹⁰6³ and 4¹5¹⁰6⁴ cages are observed in the incipient ordered hydrates, and most of them exist on the hydrate surfaces and convert into regular cages during hydrate growth.

3.3 Effects of lower subcooling and presence of MIL-101 on CH₄ hydrate formation

To reveal the effect of low subcooling, such as in CH₄ hydrate reservoirs in the South Sea of China, on CH₄ hydrate formation, the condition of 285 K and 15 MPa with a corresponding lower subcooling of 7 K was chosen for further simulation. Fig. 8 shows the evolution of the simulation system under this condition during three repeated runs. It is found that at the low subcooling of 7 K, the process of phase separation of CH₄ to form nanobubbles becomes dominant in a homogeneous CH₄ solution with x_CH4 = 0.08, while it becomes more difficult for CH₄ hydrate formation. Only at the initial simulation stage (Fig. 8A(I)–(II) for Run1, Fig. 8B(I)–(III) for Run2, and Fig. 8C(I)–(IV) for Run3), several hydrate cages (no more than 10, see Fig. S4†) form and soon decompose. In contrast, phase separation of CH₄ to form nanobubbles takes place easily in the solution. As a result, x_CH4 in the small and large cavities of MIL-101 and their outer space decreases sharply to around 0.01 (Fig. S5A–C†), while N_CH4 in nanobubbles in the outer space of MIL-101 increases rapidly (Fig. S5D–F†). Nevertheless, no obvious phase separation occurs in MIL-101 cavities, as N_CH4 therein is close to 0. Additionally, migration of CH₄ molecules from the MIL-101 cavities to the outer space is also observed (Fig. S6†), similar to the situation at 285 K and 50 MPa (Fig. 4).

Fig. 8 Phase separation of CH₄ to form nanobubbles during Run1 (A(I)–(VII)), Run2 (B(I)–(VII)) and Run3 (C(I)–(VII)) at 285 K and 15 MPa with subcooling of 7 K. Hydrate cages are colored as in Fig. 5. MIL-101 cavities, CH₄ and water molecules are shown as magenta sticks, green balls and light blue lines, respectively. The red arrows in the panels indicate the coalescence of CH₄ nanobubbles.

An interesting phenomenon is observed in the evolution of CH₄ nanobubbles at a low subcooling of 7 K (at 285 K and 15 MPa). With the rapid phase separation of CH₄ in solution, several small nanobubbles emerge in the system in the three repeated runs (Fig. 8A(III), B(II) and C(II)), and two small nanobubbles are observed to merge into a large nanobubble (Fig. 8A(III), (IV), B(II), B(III), C(II) and C(III)), and then only form two stable nanobubbles (Fig. 8A(V), B(IV) and C(IV)). Interestingly, as shown in Fig. 8A(V)–(VI) for Run1, and in Fig. 8B(IV)–B(V) for Run2, the large nanobubble gradually grows into a larger size, while the small nanobubble shrinks slowly until it completely disappears (Fig. 8A(VI)) or directly merges to the large nanobubble (Fig. 8B(VI)). Finally, the large nanobubble elongates in one direction to form an ellipsoid shape and develops into a stable cylindrical nanobuble crossing the simulation box (Fig. 8A(VII), B(VII) and S7†). Such evolution of CH₄ nanobubbles also occurs in Run3 (Fig. 8C(IV)–(VII)), though it takes a much longer time as both nanobubbles have similar sizes and are very stable. The evolution of large and small nanobubbles in solution is similar to Ostwald ripening and is to minimize the surface energy of nanobubbles. Phase separation of CH₄ in solution and coalescence of CH₄ nanobubbles in the three repeated runs at 285 K and 15 MPa are visualized in Videos S6–S8.†

To unravel the role of MIL-101 in CH₄ hydrate formation and evolution, we compare the results from three repeated runs in the presence and absence of MIL-101, respectively, under the same condition (285 K and 50 MPa with subcooling of 19 K). The main effect of MIL-101 is that its presence may promote phase separation of CH₄ in solution to form nanobubbles, and then reduce the aqueous CH₄ concentration and affect CH₄ hydrate formation. In addition, the large and small cavity structures of MIL-101 provide two nanospaces with a diameter of 3.4 and 2.9 nm in solution. In the absence of MIL-101, x_CH4 slightly decreases to 0.075 only at the initial simulation stage (Fig. 3), because a few CH₄ molecules form small clusters in solution, and no phase separation and CH₄ nanobubble formation are clearly observed in all three repeated runs (Fig. 9). In contrast, in the presence of MIL-101, phase separation of CH₄ and formation of large nanobubbles occur in the large cavity of MIL-101 and the outer space (Fig. 3 and 5). The underlying mechanism for MIL-101 to promote phase separation of CH₄ in solution may arise from two factors. The first factor is that the atoms of polar groups in MIL-101, including Cr and O in the Cr₃O trimer, O in the carboxylate groups of 1,4-BDC and F, bind surrounding water molecules (as confirmed by the high peaks between these atoms and water oxygen in the g(r) curve, see Fig. 10A), thus weakening their ability to solvate CH₄ molecules in solution. The second factor is that the methyl and benzene groups in MIL-101 adsorb CH₄ molecules, as evidenced by the peaks at around 0.4 nm in the g(r) curve for CH₄ and carbon atoms in these groups of MIL-101 (Fig. 10B), which may facilitate the formation of CH₄ nanobubbles. It is noted that the presence of MIL-101 does not change the mode for hydrate formation in homogeneous CH₄ solution, i.e., hydrate cages uniformly form at several sites in the system and develop into small clusters, which merge into one or two large hydrate clusters and then sustainedly grow to form large hydrate solids, as shown in Fig. 5 and 9. It is worth mentioning that CH₄ hydrate solids formed in the absence of MIL-101 (Fig. 9A(VI), B(VI) and C(VI)) are quite amorphous with only certain sI and sII motifs (Fig. S3C and D†), further highlighting the important role of stochasticity in the formation of crystalline hydrate from a homogeneous CH₄ solution at moderate subcooling of 19 K (at 285 K and 50 MPa). CH₄ hydrate formation in the three repeated runs in the absence of MIL-101 is visualized in Videos S9–S11.†

Fig. 9 CH₄ hydrate formation in CH₄ solution in the absence of MIL-101 during Run1 (A(I)–(VI)), Run2 (B(I)–(VI)) and Run3 (C(I)–(VI)) at 285 K and 50 MPa with subcooling of 19 K. Hydrate cages are colored as in Fig. 5. CH₄ and water molecules are shown as green balls and light blue lines, respectively.

Fig. 10 Radial distribution functions g(r) of (A) water oxygen atom (O_w) and (B) methane around specific atoms of MIL-101, which are labeled in a fragmental cluster of MIL-101 in panel (B). The g(r) curves were calculated at 285 K and 50 MPa with subcooling of 19 K.

3.4 Implications for CH₄ storage in porous materials via forming hydrates

When the adsorption-hydrate hybrid method is used for CH₄ storage, it is expected to adsorb CH₄ in the intrinsic nanopores of porous materials and form CH₄ hydrates in the intergranular space, leading to a maximum gas storage capacity. In this study, a moderate temperature of 285 K and a pressure of 50 MPa (i.e., moderate subcooling of about 19 K) are implemented for CH₄ hydrate formation in MIL-101. The main finding is that at such moderate subcooling, CH₄ hydrate formation and phase separation of CH₄ to form nanobubbles occur simultaneously in a homogeneous CH₄ solution, and these two processes compete with each other. Nevertheless, the finally formed large CH₄ nanobubbles always stay in the outer space rather than in the confined space of MIL-101 cavities. When the subcooling is lowered to 7 K (285 K and 15 MPa), the process of phase separation of CH₄ to form nanobubbles becomes dominant over CH₄ hydrate formation. In contrast, at higher subcooling of 29 K (275 K and 50 MPa), our previous study¹⁸ has shown that phase separation of CH₄ to form nanobubbles is suppressed (Fig. S8†), while CH₄ hydrate formation becomes dominant. These observations suggest that proper low subcooling (i.e., pressure-temperature conditions) should be adopted to facilitate CH₄ hydrate formation and meanwhile suppress phase separation of CH₄. Recent experimental studies also highlight the importance of proper pressure-temperature conditions for CH₄ storage in porous materials with the adsorption-hydrate hybrid method. Chen et al. observed CH₄ adsorption and hydrate formation in wetted ZIF-8 at 275.15 K and 7.0 MPa, however, no hydrate was formed at a higher temperature of 283.15 K.¹³ Koh et al. determined the dissociation temperatures of CH₄ hydrate in HKUST-1, ZIF-8 and ZIF-67 at 8.0 MPa were 283.45–283.85 K, 281.15–281.65 K and 281.75–282.05 K, respectively.^9,10 Natural gas hydrates in actual marine sediments may also form under moderate subcooling, e.g., 283.15 K and 15 MPa for hydrate reservoirs in the South Sea of China. Therefore, it is crucial to properly control operating conditions to achieve high-performance CH₄ storage in porous materials. It is noted that to ensure CH₄ hydrate formation within acceptable simulation time, this study has simulated a moderately supersaturated homogeneous CH₄ solution with x_CH4 = 0.08, i.e., CH₄ molecules were directly dispersed in water; while in previous experimental studies,^8,9,11 the cavities and intergranular space of pre-humidified MOFs were filled with water molecules. Thus, future MD studies should simulate the pre-humidified MOFs exposed to the CH₄ gas phase, to elucidate how subcooling affects the diffusion and dissolution of CH₄ into water and the subsequent nanobubble formation and hydrate formation processes in MOFs.

4. Conclusions

Aiming to reveal the effect of moderate subcooling on CH₄ hydrate formation, MD simulations have been performed to investigate CH₄ hydrate formation from a homogeneous CH₄ solution with x_CH4 = 0.08 in the presence of MIL-101 at 285 K and 50 MPa (i.e., at subcooling of 19 K). The results show that at this moderate subcooling, CH₄ hydrate formation and phase separation of CH₄ to form nanobubbles occur simultaneously, and compete with each other for aqueous CH₄ molecules. The outcome of the competition is primarily determined by the relative stabilities of CH₄ hydrate solids and CH₄ nanobubbles, which are closely related to their sizes. CH₄ hydrate formation occurs exclusively in the outer space of MIL-101 cavities, whereas phase separation of CH₄ to form nanobubbles takes place simultaneously in the MIL-101 cavities and their outer space. The small nanobubbles in the MIL-101 cavities gradually shrink and finally disappear, as CH₄ molecules therein diffuse out to facilitate the growth of large nanobubbles and large hydrates in the outer space. Finally, the formed large CH₄ nanobubbles and large CH₄ hydrates always stay in the outer space rather than in the confined space of MIL-101 cavities. Moderate subcooling appears to facilitate the formation of large ordered CH₄ hydrates containing sI and sII domains. Additionally, it is revealed that a lower subcooling of 7 K (at 285 K and 15 MPa) promotes the phase separation of CH₄ to form nanobubbles and inhibits hydrate formation in a homogeneous CH₄ solution. In the evolution of CH₄ nanobubbles during phase separation, coalescence of CH₄ nanobubbles to form large nanobubbles and an interesting phenomenon similar to Ostwald ripening are observed, i.e., the large nanobubble sustainedly grows to a larger size, while the small nanobubble shrinks gradually until completely disappears. The presence of MIL-101 can facilitate the phase separation of CH₄ in solution. These molecular insights into the effects of the degree of subcooling on CH₄ hydrate formation in MIL-101 could have important implications on optimizing pressure-temperature conditions for CH₄ storage in porous materials with adsorption-hydrate hybrid technology.

Author contributions

Zhongjin He: investigation, methodology, software, validation, visualization, formal analysis, data curation, writing – original draft, writing – review & editing, funding acquisition. Jianwen Jiang: conceptualization, methodology, supervision, writing – review and editing. Guosheng Jiang: resources, writing – review and editing. Fulong Ning: resources, writing – review and editing, project administration, funding acquisition.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (No. 42376202 and 41976203), and National Science Foundation for Distinguished Young Scholars (42225207).

References

1. A. S. Jalilov, Y. L. Li, C. Kittrell and J. M. Tour, Nat. Energy, 2017, 2, 932–938.
2. E. Soubeyrand-Lenoir, C. Vagner, J. W. Yoon, P. Bazin, F. Ragon, Y. K. Hwang, C. Serre, J. S. Chang and P. L. Llewellyn, J. Am. Chem. Soc., 2012, 134, 10174–10181.
3. L. Mu, B. Liu, H. Liu, Y. T. Yang, C. Y. Sun and G. J. Chen, J. Mater. Chem., 2012, 22, 12246–12252.
4. E. D. Sloan, Nature, 2003, 426, 353–359.
5. A. Falenty, W. F. Kuhs, M. Glockzin and G. Rehder, Energy Fuels, 2014, 28, 6275–6283.
6. L. Borchardt, M. E. Casco and J. Silvestre-Albero, Chemphyschem, 2018, 19, 1298–1314.
7. Z. Y. Wang, J. Duan, S. J. Chen, Y. Fu, X. F. Li, D. Wang, M. Zhang, Z. Q. Zhang, D. D. Liu and F. H. Wang, J. Energy Storage, 2022, 50, 104195.
8. M. E. Casco, F. Rey, J. L. Jordá, S. Rudic, F. Fauth, M. Martínez-Escandell, F. Rodríguez-Reinoso, E. V. Ramos-Fernández and J. Silvestre-Albero, Chem. Sci., 2016, 7, 3658–3666.
9. S. Denning, A. A. A. Majid, J. M. Lucero, J. M. Crawford, M. A. Carreon and C. A. Koh, ACS Appl. Mat. Interfaces, 2020, 12, 53510–53518.
10. S. Denning, A. A. A. Majid, J. M. Lucero, J. M. Crawford, M. A. Carreon and C. A. Koh, ACS Sustainable Chem. Eng., 2021, 9, 9001–9010.
11. C. Cuadrado-Collados, G. Mouchaham, L. Daemen, Y. Q. Cheng, A. Ramirez-Cuesta, H. Aggarwal, A. Missyul, M. Eddaoudi, Y. Belmabkhout and J. Silvestre-Albero, J. Am. Chem. Soc., 2020, 142, 13391–13397.
12. J. Duan, Q. C. Li, Y. Fu, S. J. Chen, Y. X. Zhang and D. D. Liu, Fuel, 2022, 322, 124266.
13. S. J. Chen, D. Wang, Z. Y. Wang, Y. Fu, Y. H. Xu and D. D. Liu, Energy, 2023, 283, 129013.
14. H. Liu, S. Y. Zhan, P. Guo, S. S. Fan and S. L. Zhang, Chem. Eng. J., 2018, 349, 775–781.
15. D. Kim, Y. H. Ahn and H. Lee, J. Chem. Eng. Data, 2015, 60, 2178–2185.
16. S. Denning, J. M. Lucero, A. A. A. Majid, J. M. Crawford, M. A. Carreon and C. A. Koh, J. Phys. Chem. C, 2021, 125, 20512–20521.
17. S. Denning, A. A. A. Majid, J. M. Crawford, J. D. Wells, M. A. Carreon and C. A. Koh, Fuel, 2022, 325, 124920.
18. Z. J. He, K. Zhang and J. W. Jiang, J. Phys. Chem. Lett., 2019, 10, 7002–7008.
19. Z. Y. Wang, J. Duan, S. J. Chen, Y. Fu, Y. X. Zhang, D. Wang, J. L. Pei and D. D. Liu, Chem. Eng. J., 2022, 430, 132901.
20. J. Duan, X. Y. Jiang, Y. Fu, S. J. Chen and M. C. Zi, Fuel, 2023, 337, 126851.
21. Z. Li, N. Li, J. Y. Kan, B. Liu and G. J. Chen, Fuel, 2023, 351, 129055.
22. R. Boswell, Science, 2009, 325, 957–958.
23. Z. R. Chong, S. H. B. Yang, P. Babu, P. Linga and X.-S. Li, Appl. Energy, 2016, 162, 1633–1652.
24. M. B. Clennell, P. Henry, M. Hovland, J. S. Booth and M. Thomas, Ann. N.Y. Acad. Sci., 1999, 912, 887.
25. Z. J. He, P. Linga and J. W. Jiang, Langmuir, 2017, 33, 11956–11967.
26. Z. He, F. Mi and F. Ning, Energy, 2021, 234, 121260.
27. Z. J. He, F. L. Ning, F. Y. Mi., B. Fang and G. S. Jiang, Energy Fuels, 2022, 36, 3628–3639.
28. Z. Zhang, P. G. Kusalik, C. Liu and N. Wu, Energy, 2023, 285, 129414.
29. S. Takeya, H. Fujihisa, Y. Gotoh, V. Istomin, E. Chuvilin, H. Sakagami and A. Hachikubo, J. Phys. Chem. C, 2013, 117, 7081–7085.
30. M. R. Walsh, C. A. Koh, E. D. Sloan, A. K. Sum and D. T. Wu, Science, 2009, 326, 1095–1098.
31. Z. J. He, K. M. Gupta, P. Linga and J. W. Jiang, J. Phys. Chem. C, 2016, 120, 25225–25236.
32. G. Ferey, C. Mellot-Draznieks, C. Serre, F. Millange, J. Dutour, S. Surble and I. Margiolaki, Science, 2005, 309, 2040–2042.
33. Y. F. Chen, R. Babarao, S. I. Sandler and J. W. Jiang, Langmuir, 2010, 26, 8743–8750.
34. J. L. F. Abascal, E. Sanz, R. G. Fernandez and C. Vega, J. Chem. Phys., 2005, 122, 234511.
35. W. L. Jorgensen, J. D. Madura and C. J. Swenson, J. Am. Chem. Soc., 1984, 106, 6638–6646.
36. A. K. Rappe, C. J. Casewit, K. S. Colwell, W. A. Goddard and W. M. Skiff, J. Am. Chem. Soc., 1992, 114, 10024–10035.
37. T. Darden, D. York and L. Pedersen, J. Chem. Phys., 1993, 98, 10089–10092.
38. M. Parrinello and A. Rahman, Phys. Rev. Lett., 1980, 45, 1196–1199.
39. S. Nose, Mol. Phys., 1984, 52, 255–268.
40. S. Sarupria and P. G. Debenedetti, J. Phys. Chem. Lett., 2012, 3, 2942–2947.
41. X. W. Qin, Q. Y. Liang, J. L. Ye, L. Yang, H. J. Qiu, W. W. Xie, J. Q. Liang, J. A. Lu, C. Lu, H. L. Lu, B. J. Ma, Z. G. Kuang, J. G. Wei, H. F. Lu and B. B. Kou, Appl. Energy, 2020, 278, 115649.
42. B. Hess, C. Kutzner, D. van der Spoel and E. Lindahl, J. Chem. Theory Comput., 2008, 4, 435–447.
43. M. M. Conde and C. Vega, J. Chem. Phys., 2010, 133.
44. L. A. Báez and P. Clancy, Ann. N.Y. Acad. Sci., 1994, 715, 177.
45. L. C. Jacobson, W. Hujo and V. Molinero, J. Phys. Chem. B, 2009, 113, 10298–10307.
46. G. D. Zhang, M. T. Sun, B. J. Liu and F. Wang, Chem. Eng. J., 2020, 396, 125256.
47. E. D. Sloan and C. A. Koh, Clathrate Hydrates of Natural Gases, CRC Press, Boca Raton, FL, 3rd edn, 2008.
48. J. M. Schicks and J. A. Ripmeester, Angew. Chem., Int. Ed., 2004, 43, 3310–3313.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d3ta06952d

---