Effects of nanobubbles on methane hydrate dissociation A molecular simulation study

Hydrate dissociation is often accompanied by the formation of nanobubbles. Knowledge of the effects of nanobubbles on hydrate dissociation is essential for understanding the dynamic behavior of the hydrate phase change and improving the gas production efficiency. Here, molecular dynamics simulations were performed to study the methane hydrate dissociation kinetics with and without a pre-existing methane nanobubble. The results show that the hydrate cluster in the liquid phase dissociates layer-by-layer. This process is shown to be independent of the temperature and nanobubble presence at the simulation conditions. Hydrate dissociation does not always lead to nanobubble formation because the supersaturated methane solution can be stable for a long time. A steep methane concentration gradient was observed between the hydrate cluster surface and the methane nanobubble, which can enhance the directional migration of methane and effectively minimize the methane concentration in the liquid phase, thereby increasing the driving force for the hydrate dissociation. Our findings indicate that the presence of a nanobubble near the hydrate surface does not decrease the activation energy of hydrate dissociation, but it can increase the intrinsic decomposition rate. The average hydrate dissociation rate is linearly correlated with the mass flow rate towards the nanobubble. The mass flow rate is determined by the nanobubble size and hydrate-nanobubble distance. Our findings contribute to the fundamental understanding of the dissociation mechanism of gas hydrates in the liquid phase, which is crucial for the design and optimization of efficient gas hydrate production techniques


Introduction
Natural gas hydrates (NGHs) are non-stoichiometric, ice-like compounds in which guest molecules (primarily methane) are encapsulated in polyhedral cages formed by the hydrogen-bonded (H-bonded) concomitant host water molecules at low temperature and high pressure conditions [1].Three typical hydrate crystal structures (i.e, cubic sI, cubic sII, and hexagonal sH) have been identified.Methane hydrates often form cubic sI crystalline structures [2].NGHs have been extensively studied due to their unfavorable presence in offshore oil and gas pipelines, which can induce serious flow assurance problems in petroleum transmission operations [3,4].NGHs are considered as an unconventional energy resource due to their abundant reserves (roughly twice the energy storage than that of all other fossil fuels combined [5,6]) in marine sediments, which account for more than 90% of NGH reserves, as well as in permafrost regions [7][8][9].Therefore, the exploitation of NGHs can play an important role in alleviating the energy crisis caused by the rapidly increasing energy consumption and geopolitical instability.Depressurization and thermal stimulation are considered as two effective technologies to extract natural gas from hydrates occurring in marine deposits [10].Nevertheless, natural gas production from hydrates is a complex process that requires knowledge on thermophysics, geology, and mechanics.The volatile reservoir properties and reservoir physical property in poor create a major bottleneck for the commercially exploitation of natural gas from NGH reservoirs [11].Hydrates have great potential in other industrial applications, such as carbon dioxide (CO 2 ) capture and sequestration [12], hydrogen storage [13], seawater desalination [14], wastewater treatment [15], and gas transport [16].The basis for accelerating hydrate industrialization processes should be based on the knowledge of the kinetics of hydrate phase changes for several conditions, especially in the mining industry.
The methane to water ratio reaches 1:5.75 in the sI methane hydrate [1].The number of methane molecules released from the hydrate cage during mining is limited by the low solubility of methane in water.Consequently, most methane molecules aggregate and form nanobubbles after supersaturation to minimize water-methane hydrophobic interactions.Methane nanobubbles in water, formed due to the dissociation of methane hydrates, have been observed using transmission electron microscopy (TEM) by Uchida et al. [17].Daisuke et al. [18] visually observed the formation of small methane bubbles during methane hydrate dissociation in a glass micromodel.Yang et al. [19] observed that during xenon hydrate dissociation in a sedimentary matrix, xenon gas was metastably enriched in the water phase.Wu et al. [20] investigated the physicochemical characteristics and stability of nanobubbles generated by the dissociation of CO 2 hydrates.Such nanobubbles in the bulk solution have a long lifetime and are of small size (the diameter is usually less than 1 μm) [21,22].These experimental results directly or indirectly show that nanobubbles are formed during hydrate dissociation.Besides the small size, nanobubbles also have other characteristics, such as negligible rise speed in water, and extremely high internal pressure.Nanobubbles play an important role in the kinetics of hydrate formation and dissociation [23], and have attracted scientific attention for their potential applications in hydrate science [20].
Limited by the low spatial and temporal resolution, it is difficult to experimentally study the dynamic behavior of gas hydrate dissociation and the formation of nanobubbles in the liquid phase at the molecular level.Molecular Dynamics (MD) simulation is a complementary tool to experimental measurements, and can be used to study the structure and the dynamic behavior of nanobubbles during the hydrate dissociation process at the molecular level.In MD simulations of the methane hydrate dissociation process [23][24][25][26][27][28][29][30], nanobubble nucleation occurs during hydrate dissociation.Ripmeester et al. [24] performed MD simulation and showed that the formation of nanobubbles near the solid-liquid interface affects the mass transfer rate between the solid and liquid phases.Yagasaki et al. [25] reported that the hydrate dissociation process is slow before the formation of nanobubbles takes place, and the dissociation is enhanced after the nanobubble formation in the liquid phase.MD studies have shown that the assembly and stabilization of the nanobubbles may enable the absorption of the dissolved methane to decrease the concentration of methane in the liquid phase [25,27,31], and also alter the water structure near the gas/water interface and destabilize hydrate nuclei [32].These MD simulation results also show that the formation of nanobubbles can increase the decomposition of the residual hydrate.
Some studies show the formation of nanobubble has an opposite effect on hydrate dissociation.Uddin et al. [26] reported that the accumulation of nanobubbles at the hydrate interface could significantly decrease the dissociation rate due to the thermal resistance effect of nanobubbles.Enough heat should be provided to keep the hydrate dissociation rate constant, which means the heat transfer must be considered during the hydrate dissociation process [33].English and Phelan [34] proposed a simple model with coupled mass and heat transfer to describe the hydrate dissociation by fitting to the initial break-up rates computed from MD simulations.The negative effects of foams formation (containing a variety of sizes bubbles) on hydrate dissociation have been confirmed in macroscopic experiments.Wang et al. [35] observed that microbubbles can reduce the hydrate dissociation rate by using microfluidics, because the pressure propagation and heat transfer could be impeded by the large distribution of bubbles in the liquid.Uddin et al. [36] showed that the removal of newly formed bubbles could contribute to the intrinsic hydrate dissociation.
The combination of gas supersaturation and diffusion controls the formation of nanobubbles during hydrate dissociation [27,37], which strongly depends on temperature [23,27,38,39].Once nanobubbles are formed, they grow, move, and merge [23,40], and eventually evolve into two states, i.e., either the nanobubbles remain stable [30] or dissolve in water [23,30].Bagherzadeh et al. [23] found that dissolved methane could aggregate and develop into nanobubbles when the methane mole fraction reaches 0.044 in the homogeneous liquid phase.The gas concentration depends on the methane release rate (hydrate dissociation rate) and the gas diffusion in liquid water.Diffusion can be enhanced using a gas reservoir near the hydrate before the formation of nanobubbles [23,27,38].The rich gas source can promote the nucleation process of hydrate reformation at suitable conditions.A hypothesis of memory effects in hydrate formation was proposed based on nanobubbles [19,20,41].
In nature, hydrate sources often occur in systems with porous media.Solid and confined spaces affect the mass transfer and control nanobubble formation so that bulk and surface nanobubbles can exist in an actual hydrate mining environment.Both bulk and surface nanobubbles have been observed using TEM [17].Previous MD simulations showed that spherical cap-shaped nanobubbles are formed during hydrate dissociation within the porous structure.These nanobubbles play an important role in the further dissociation of the residual hydrate [40,[42][43][44].Thermodynamic inhibitors, such as methanol (CH 3 OH) and sodium chloride (NaCl), enhance the hydrate dissociation process by facilitating nanobubble formation by different mechanisms.Methanol stabilizes small clusters of methane molecules through its amphiphilic character [31,45].The presence of NaCl increases the hydrophobic interactions between dissolved methane and induced methane in the solution, thereby, enhancing nanobubble formation [30,31,45].CO 2 gas can also promote the formation of methane nanobubbles by inducing early nucleation [45,46].Static electric fields have been shown to enhance bulk nanobubble formation at massively elevated levels, and exceeding Henry's law solubilities [47].
During hydrate dissociation, the hydrate cluster acts as a large gas source in the liquid phase such that the water solution is in a nonequilibrium phase.A gas concentration gradient along the aqueous layer in a direction perpendicular to the hydrate dissociation surface is present [19,48].Nanobubbles can form near the hydrate surface randomly at varying distances and of any time [23][24][25][26]31,42].Therefore, the formation of nanobubbles can be quite unpredictable.It is very difficult to quantitatively study the nanobubble effects on the hydrate dissociation process.On this work, MD simulations are performed to investigate nanobubble and temperature effects on dynamic behavior and mass transfer of hydrate dissociation process.
This paper is organized as follows.In Section 2, we explain how to build two types of initial hydrate dissociation systems without and with a pre-existing nanobubble, respectively.This is followed by the details of the MD simulations.In Section 3, we present the kinetics of hydrate dissociation and nanobubble growth process, and calculate the hydrate dissociation rate per unit of hydrate surface area for different temperatures for the two types of systems.We also discuss the mechanism of nanobubble effects on hydrate dissociation.Finally, the conclusions of this study are summarized in Section 4. Our findings will help in understanding the hydrate dissociation process, as well as provide scientific reasoning and theoretical support for evaluating subsequent gas production during mining.

Initial MD simulation configurations
To systematically study the effects of the presence of nanobubbles on gas hydrate dissociation in a liquid water environment at the molecular level, we build two different systems.The first system (denoted as the HW system) is composed of a hydrate cluster in a liquid water phase.The B. Fang et al. second system (denoted as the HW-nb system) contains a cylindrical methane nanobubble close to the undissociated hydrate cluster.Both systems are shown in Fig. 1.To better study the effects of relative location and size of the nanobubble on hydrate dissociation, we varied the HW-nb system as follows: First, the HW-nb system was divided into four areas along the X-axis (Areas 1-4) as shown in Fig. 1b.The solid hydrate phase is placed in Area 1, and a large methane nanobubble is placed in Area 3 (HW-nb-L system).The distance between the large nanobubble center and the hydrate cluster center is about 21.2 nm.To evaluate the effect of the bubble size on hydrate dissociation, we vary the nanobubble size by decreasing the number of methane molecules within the nanobubble.This results in two new configurations with medium and small nanobubbles; these systems are labeled as HW-nb-M and HW-nb-S, respectively.For the final system, we shifted the largesized nanobubble from Area 3 to Area 2 (HW-nb-LC).To establish a relation between the hydrate dissociation rate and the hydrate-nanobubble distance, the distance from hydrate center to nanobubble center was shortened to 10.2 nm, which is closer to the hydrate surface.The detailed characteristics of these systems are summarized in Table 1.
The dimensions of the initial simulation boxes are 43.08×21.54×3.61nm 3 , with a hydrate cluster located in Area 1.A cylindrical gas bubble is added to the HW-nb system, and the remaining space was occupied by the water molecules, i.e., ca.400,000 atoms are used in each system.The cylindrical nanobubble is used to ensure that the box size in the z-direction can be small, so that the number of water molecules can be relatively small.The methane hydrate cluster consists of 9×9×3 cubic unit cells.A unit cell of the sI hydrate is constructed as follows: the positions of oxygen atoms in the host water molecule are determined based on X-ray crystallography [49], and the hydrogen atoms of the host water molecule are adjusted to minimize the potential energy and net unit cell dipole moment based on the Bernal-Fowler rule [50].The guest molecules (methane) are placed at the center of the polyhedral cages constructed using the host water molecules and hydrogen-bond networks (methane mole fraction x CH v = 0.148).The lattice parameter of the sI hydrate unit cell is 12.03 Å [50].

Force fields and simulation details
The TIP4P/2005 model was used to represent water molecules [51].The "Optimized Potentials for Liquid Simulations-United Atom" (OPLS-UA) force field was used to model methane molecules [52].This combination of force fields has been shown to be accurate for investigating phase change kinetic behavior of methane hydrates [25,53].For the nonbonded intermolecular interactions, Lennard-Jones and Coulombic potentials were used as per the original force field [51,52].The cutoff radius for the nonbonded interactions was set to 1.2 nm.
Long-range electrostatic interactions were calculated using the particle-mesh Ewald summation method [54].The leapfrog algorithm was used to integrate the Newton's equations of motion with a time step of 2 fs [55].The energies of the initial MD simulation configurations were first minimized using a conjugate gradient algorithm [56].This was followed by a long isothermal − isobaric (NPT) 200 ns simulation (in NPT-MD simulations, the total number of particles N, pressure P and temperature T are controlled and remain constant), which was The initial HW-nb configuration was divided into four areas along the X-axis.In the HW-nb-L system, a cylindrical pre-existing nanobubble was placed in Area 3 with a 5.24-nm radius.In the HW-nb-M system, a medium-sized nanobubble with a 3.49-nm radius was present in area 3.In the HW-nb-S system, a smallsized nanobubble with a 2.44-nm radius was placed in area 3 in the bulk phase.In the HW-nb-LC system, we decreased the nanobubble-hydrate distance based on the configuration of the HW-nb-L system.Methane molecules in pre-existing nanobubble and hydrate are colored by cyan and purple, respectively.For the unit sI hydrate cell, the 5 12 cage is shown in blue and the 5 12 6 2 cage is shown in green.

Table 1
Initial configurations of the HW and HW-nb systems with different nanobubble sizes and locations.All MD simulations were performed at 0.1 MPa.conducted to study the hydrate dissociation process and the effects of nanobubbles on this process.The velocity rescaling method and the Parrinello-Rahman extended-ensemble method were used for temperature and pressure coupling, respectively, the thermostat constant was 0.1 ps and the barostat constant was 1.0 ps [57,58].To evaluate the effect of the temperature on the dissociation process, MD simulations were performed for a series of temperatures (292, 302, 312, and 322 K), while maintaining the pressure at 0.1 MPa.Periodic boundary conditions were used along the three directions.All simulations were performed using the Gromacs 5.0.7 package [59,60].

Criteria to distinguish phase
During hydrate dissociation process, methane molecules will be released into the liquid water.Both the hydrate cluster and liquid water phases are composed of water and methane molecules.Therefore, it is necessary to develop a criterion for determining which water and methane molecules belong to the solid hydrate phase (or the liquid phase) during the dynamic dissociation process.Here, we used the F 3 order parameter proposed by Baez and Clancy [61].This parameter is based on the three-body (three water oxygen atoms) configuration that quantifies the deviation between a standard tetrahedron structure and a tetrahedron comprising a central water oxygen atom.The water oxygen tetrahedron configuration is different in the liquid (F 3 = 0.1) and the solid (F 3 = 0 e.g., hydrate, ice) phases [62].The expression of F 3 equals: where θ jik denotes the angle formed by i, j, k oxygen atoms with the oxygen atom i in center.The distance between two other surrounding water oxygen atoms j and k is within the range of 3.5 Å from the central oxygen atom i.The bracket <…> indicates an average over all the F 3 of the three-oxygen structure for the central oxygen atom i.Hence, we assume that the water molecule is in a solid phase (hydrate) if F 3 is smaller than 0.05, otherwise, water molecule is in the liquid phase.We assume that a methane molecule is in the hydrate phase when the number of surrounding water molecules in solid phase within a range of 5.5 Å is larger than 15.This is distance of the first water shell from the methane molecule in the hydrate phase [25,63].The number of methane molecules in the solid phase is used to quantify the undissociated hydrate during dissociation.
To quantify the effect of nanobubbles on hydrate dissociation, variations in the nanobubble size with time should be determined.Due to the difference in the density of liquid water and methane gas bubbles, we calculated the density map of the systems in the X-Y plane by dividing the simulation boxes into bins of 1 Å × 1 Å.The boundary between the gas nanobubble and liquid bulk phase was defined as the locus of squares with a density of 0.5 g/cm 3 , which corresponds to half the density of liquid water at normal atmospheric pressure and temperature.In Fig. 2, we show a snapshot (Fig. 2a) for a HW-nb-L system at 2 ns at 292 K and a pressure of 0.1 MPa, as well as the density map (Fig. 2b) of the corresponding system.By comparation of snapshots and the density map, it was verified that the nanobubble phase boundary is well-determined.The radius of the nanobubble can be estimated by calculating the area of the nanobubble.

Results and discussion
Gas hydrate dissociation is an endothermic process since heat is required for the breaking of the H-bonded network.The variation in the potential energy with time for all systems at different temperatures is shown in Fig. 3.The potential energy initially increases (until the As shown in Fig. 3, the time required for hydrate dissociation of each simulation is marked.This suggests that the duration of hydrate dissociation varies, although the size of the hydrate cluster remains the same in each system.As shown in Fig. 3(a and b), hydrate cluster dissociation is faster as temperature increases.In the presence of a methane nanobubble, it takes less time for hydrate dissociation to finish at the same temperature and pressure.As shown in Fig. 3(c), the nanobubble size and hydrate-nanobubble distance affect the hydrate dissociation process.In practice, depressurization technology is the most practical method for gas hydrate exploitation.As shown in Fig. 3, the potential energy increases during the hydrate dissociation as this process demands heat that is supplied via the thermostat.If the heat supply is insufficient during hydrate dissociation, the increase of the potential energy of the system should be supplied from the kinetic energy of the system, which leads to a large temperature drop of the system with hydrate dissociation [40,42,43].The temperature drop is very unfavorable for gas production.Therefore, heat should be added during gas production using depressurization.

Kinetics of hydrate dissociation in the presence of a nanobubble
Typical simulation snapshots showing the hydrate dissociation process for the system with and without the nanobubble at 292 K and 312 K are shown in Fig. 4. Initially, the undecomposed hydrate phases are in the shape of a square in the X-Y plane, due to the initial simulation system construction settings (Fig. 1).At this phase, no methane molecules have escaped from the H-bonded cages of the hydrate crystal into the bulk liquid phase.After the initial hydrate dissociation, the shape of the hydrate cluster changed from a square into a roughly circular crosssection, which can be attributed to the Gibbs-Thomson effect [64].The corners of the hydrate cluster having small curvature dissociate faster than the planar parts, and the spherical-shaped residual hydrate shrinks in a layer-by-layer manner with a curved decomposition front until the dissociation is completed.In our simulations, the temperatures were set to 292 K, 302 K, 312 K, and 322 K.For the hydrate cluster, the temperature of the inner hydrate was high enough to induce a phase change.However, the inner hydrate could not decompose prior to hydrate surface decomposition.The reason behind this phenomenon is that unreleased methane molecules can stabilize the hydrate cage formed by the H-bonds of water molecules.Therefore, due to mass transfer limitations of inner gas molecules the hydrate is maintained in a crystalline state when the temperature is higher than phase equilibration conditions [65][66][67].For the hydrate surface cage, if a few water molecules escaped from the complete H-bonded cage network, gas molecules can easily diffuse out of the incomplete cage and facilitate the dissociation process [67,68].Therefore, the hydrate dissociation process is heterogeneous.Hydrates dissociate from the outer to the inner surface in sequence.This is similar to what is reported in previous studies [40,61,69].Regarding hydrate dissociation with a pre-existing nanobubble in the bulk phase, the dissociation process did not change and decomposed from the outer to the inner surface at the same simulation conditions, as shown in Fig. 4. In hydrate dissociation dynamics, this heterogeneity is not unique.If the surface mass transfer is limited and the system is heated rapidly, the hydrate dissociation process can occur in a homogenous manner [70,71].However, this type of homogeneity was not observed in our simulations, even for a wide temperature range (292-322 K).
Our simulation results show no new nanobubble (methane cluster) formation in systems without (Fig. 4a -c) and with (Fig. 4d -f) a preexisting nanobubble at 292 K during hydrate cluster dissociation.The phenomenon of new nanobubble formation is shown to be common in previous simulations and experiments involving hydrate dissociation [23,26,27,38].With hydrate dissociation, methane molecules escape from the cages in hydrate surface and diffuse into the liquid bulk phase (Fig. 4).If the number of released methane molecules in a specific area is larger than the solubility limit, methane molecules can form a cluster and evolve into a bubble if the process of methane release continues, such as in the hydrate dissociation process [38].In systems without a pre-existing nanobubble, the slow dissociation rate facilitates sufficient diffusion of the released methane molecules and prevents aggregation into bubbles.The average mole fraction of methane in water after hydrate dissociation is 0.018 (there were 1944 methane molecules and 11,178 water molecules in the initial hydrate cluster, and 96,910 water molecules in the bulk phase), which is much larger than the solubility of methane in water (mole fraction of ca.7×10 − 6 ) [72] and less than the critical methane fraction in a homogeneous liquid phase (mole fraction of ca.4.4×10 − 2 ) [23].This indicates that the methane concentration was supersaturated but did not reach the barrier for bubble formation.When we increase the dissociation temperature (312 K), hydrate dissociation was accelerated and new nanobubbles appeared around the residual hydrate cluster in the HW system as can be seen in Fig. 4(h and i).The results for the HW-nb-L system were comparable to those of the HW system.There was no new bubble formation at 292 K, while several newly formed bubbles were observed at a high temperature (312 K) in the bulk phase (Fig. 4(gl)).These results indicate that there is no connection between hydrate dissociation and nanobubble formation, because the methane molecules can be dissolved in a supersaturated liquid at the timescale of the MD simulation.
As shown in Fig. 4(df) and (j-i), the pre-existing nanobubble grows as the number of methane molecules decrease in the bulk phase, particularly in the area around the nanobubble.The nanobubble moves similar to the Brownian motion, and more released methane molecules can diffuse into the nanobubbles.In Section 3.3, we have already discussed the effects of nanobubbles on the hydrate dissociation rate by a quantitative analysis.Based on the simulations, we found that nanobubbles do not necessarily form in the liquid phase during hydrate dissociation.When temperature increases, the probability of nanobubble formation increases because more methane molecules are released in a short time.The rapid increase in the number of methane molecules without sufficient diffusion leads to nanobubble formation near the residue hydrate cluster.This suggests that nanobubble formation strongly depends on the local methane density in water.The preexisting nanobubble in water can grow and move during the hydrate dissociation process, and it appears to play a minor role in hydrate dissociation manner and new nanobubble formation.When the temperature is further increased to 322 K, several new nanobubbles form in a shorter time compared to those formed at 312 K in the system without (Fig. 5a) and with (Fig. 5b) a pre-existing nanobubble.Only a small number of released methane molecules can be found in the initially existing nanobubble, indicating that new methane bubbles form quickly before methane diffuses into the existing nanobubble.
For systems without and with the nanobubble, the number of methane molecules in the hydrate cluster as a function of time is shown in Fig. 6(a).The exponential-like decay indicates that the hydrate dissociation rate (-dn/dt) is not constant (as shown in Fig. S1 of Supporting Information) and may be influenced by the decrease in the area of hydrate dissociation front (shown in Fig. 4).The hydrate dissociation time in both systems is shown in Fig. 6(b) as a function of temperature.The dissociation time decreases exponentially with increasing temperature.This clearly shows that temperature plays an important role in hydrate dissociation within a specific temperature range.As shown in Fig. 6(b), the hydrate dissociation time in the system with a pre-existing nanobubble is shorter than that of the system without a nanobubble at same temperature and pressure.The dissociation time difference decreases with an increase in the temperature, indicating that the existence of a nanobubble around the hydrate cluster can facilitate the decomposition process at the same temperature and pressure.This effect become less pronounced as the temperature increases.The trend in the variation of the number of methane molecules in the hydrate cluster with time (Fig. 6a) is the opposite to the variation in the potential energy of the system with time (Fig. 3).For each simulation, the duration of hydrate dissociation shown in Fig. 6(b) is close agreement with the data in Fig. 3.These results also show that the distinguished criteria given for the methane molecule in hydrate phase is reasonable.

Dissociation rate per unit area of hydrate surface in the presence of a nanobubble
Because the residual hydrate clusters shrink during the dissociation process, and the surface area of dissociation front decreases nonlinearly against time, the hydrate dissociation rate is not constant.To remove the effect of surface area change, the gradient of the number of methane molecules in hydrate cluster can be calculated from [25] where N(t) is the number of methane molecules in the hydrate phase, S (t) is the dissociation front surface area of the hydrate cluster, and k(t) is the dissociation rate per unit area.As discussed in Section 3.1, hydrates dissociate with an irregular dissociation front surface.During this process, the residual hydrate cluster rotates and translates.As shown in Fig. 4, hydrate clusters have a cylindrical shape except for the initial stage.Therefore, we assumed that where r(t) is the radius of the hydrate cluster, l z is the length of the simulation box along the Z-axis, and ρ is the density of methane in the hydrate crystal.When ρ is considered a constant, we obtain The computed dissociation rates k(t) are shown in Fig. 7 for the system without a pre-existing nanobubble at 292 K (other computed rates for the system without and with a nanobubble at different temperature are shown in Fig. S2 in Supporting information).The shapes of k(t) curves are similar: first, the dissociation rates are high and decrease rapidly, second, the dissociation rates of the hydrate remain constant with a small fluctuation and speed up the dissociation process at the final stage.In the initial stage, the initial configuration of the hydrate phase was prepared as a cubic and not a spherical structure (Fig. 1), and methane molecules were present in an incomplete cage on the hydrate surface.These methane molecules can escape easily from the open cages.This, together with the fast dissociation rate of the hydrate cube corners, results in the rapid dissociation of the hydrate in the initial stages of dissociation.Subsequently, the hydrate dissociation rates per Fig. 5. Simulation snapshots of (a) HW and (b) HW-nb-L systems at t = 5 ns at 322 K in the X-Y plane.After 5 ns, the hydrate dissociation process was completed and several new nanobubbles were formed (represented in purple).The color scheme is the same as in Fig. 4.  surface area fluctuated around a constant value.The decomposition rate becomes stable, thus, we obtained the hydrate dissociation rates per surface area at different temperatures by averaging k(t) in this period.In the final stage, the hydrate dissociation rates rapidly increased because the stability of the hydrate weakened due to the small residual hydrate cluster, leading to a sudden collapse of the hydrate cluster.Snapshots for a 5 12 6 2 cage decomposition in the final stage are shown in Fig. 7.The collapse of this cage begins with multiple sides of cage and spreads to the whole cage rapidly without the support of host-host interactions, in sharp contrast to the two-step pathway in hydrate dissociation at other dissociation stages [68].The dissociation rate increases with the increasing of temperature in both systems as shown in Fig. S2.
Based on the average hydrate dissociation rate in the stabilization phase, the hydrate dissociation rates per surface area in the two systems were calculated at different temperatures according to the logarithmic form of the Arrhenius equation [73] where E a is the activation energy and k 0 is the intrinsic decomposition rate (pre-exponential factor).The Arrhenius plot of k(T) in Fig. 8 shows a linear relationship between k(T) and reciprocal of temperature.By performing linear regression, the activation energy computed for the hydrate dissociation process in the liquid phase was calculated as 96.9 ±2.8 kJ/mol for the HW system without the pre-existing nanobubble.This value is nearly equal to the MD simulation result reported by Yagasaki (96 kJ/mol [25]) but higher than the experimentally measured value of 81 kJ/mol [74].This discrepancy may be due to the methane occupancy in the hydrate cage.The typical experimental sample is approximately 95% occupied, while in our simulations, the hydrate cage is fully occupied.The difference in occupancy induced a scaling factor (0.89) in the calculation of the activation energy based on the MD simulation results by English et al. [34].In the presence of a nanobubble (HW-nb-L system), the activation energy required for the dissociation process was calculated as 97.8±4.0 kJ/mol, which is slightly higher than that required in a system without a nanobubble.This indicates that the nanobubbles play a small role in the activation energy of the hydrate dissociation process.However, the intrinsic decomposition rate k 0 was different for the two systems.The calculated values of ln k 0 are 45.17 and 45.73 for the two systems, respectively.The weak temperature dependence of the pre-exponential factor (ln k 0 ) is negligible due to the narrow temperature range (292-322 K).This indicates that the pre-existing nanobubble can change the intrinsic properties of phase transition and plays an important role in hydrate dissociation.It is important to note that the hydrate dissociation process is also determined by the pressure of the system [75].Due to the low compressibility of liquid water, it is difficult to study the pressure effects on hydrate dissociation in liquid water by MD simulations, therefore, the pressures of the systems are kept at atmospheric pressure in our simulations.
Our results show that nanobubbles can increase the hydrate dissociation rate by approximately 3% in the HW-nb-L system at the simulation conditions.The driving force for hydrate dissociation is proportional to the difference (△μ) between the chemical potentials of the solid (μ Hydrate ) and methane solution (μ CH4 + μ water ) [76].With hydrate dissociation, methane is released from the incomplete H-bonded cages, and the number of methane molecules increase in the liquid phase around the residual hydrate cluster.The chemical potential of methane in an ideal solution equals where μ θ CH4 (T, p) is chemical potential in the pure state, and c CH4 is mole fraction of methane.Prior to nanobubble formation, the chemical potential of methane increases with the methane concentration due to the hydrate dissociation process along with a decrease in the driving force for phase change.In Fig. 9, we show the average methane density and concentration along the X-axis as a function of the simulation time for the systems without and with the nanobubble.In the HW system (Fig. 9a), in the initial stage, the density of methane near the hydrate was larger than that of the region far from the hydrate cluster, because methane molecules are released from the hydrate and diffuse into the liquid water phase.The methane density tends to balance along the Xaxis in the liquid phase because of methane transfer.The overall methane density in the bulk phase increased due to the continuing methane emission during the dissociation process.In the HW-nb-L system (with the large-sized nanobubble) (Fig. 9(b)), a strong density (concentration) gradient of methane in the liquid phase between the residual hydrate and nanobubble was found.Similar results were also observed experimentally [19,48].In Fig. 10, we present the gradients of methane concentration from the hydrate cluster surface to the middle of the liquid phase and the pre-existing nanobubble surface in the HW and HW-nb-L systems as a function of simulation time.Due to the nanobubble, the negative derivative of the methane concentration near the residual hydrate surface significantly increased in the HW-nb-L system compared to the HW system.The methane diffusion process in water obeys Fick's law: J = -D▽c CH4 , (J is the diffusion flux, D is the diffusion coefficient).The large negative derivative in the HW-nb-L system indicates that a number of released methane molecules diffuse into the pre-existing nanobubble.The methane concentration in the liquid phase near the residual hydrate decreases, thereby decreasing the chemical potential of methane compared with that observed in the HW system and increasing the driving force.
As shown in Fig. 4, the residual hydrates remaining in the spherical shape as hydrate clusters dissociated changed the curvature of the hydrate surface in the X-Y plane.According to the Gibbs-Thomson equation [25,77], the difference in the chemical potential of a spherical crystal cluster of one component system is expressed as: where μ ∞ is the chemical potential with a flat surface and γ is the specific surface-free energy.Therefore, the decrease in the curvature of the solid phase may increase the chemical potential difference and the driving force with hydrate dissociation.Plots of k(t) as a function of 1/N 1/2 as shown in Fig. 11.k(t) does not increase linearly with 1/N 1/2 , indicating that the change in the curvature of the hydrate surface may play a minor role in governing the hydrate dissociation rate.This result may be attributed to the inhibition induced by an increase in the methane  concentration in the liquid phase.To some extent, these two functions are interdependent and interactive, suggesting that the dissociation rate remains virtually unchanged.

Effect of the nanobubble size and distance on hydrate dissociation
During hydrate dissociation, nanobubbles of different sizes form near the hydrate surface at different distances [23][24][25][26]31,42].This phenomenon was also observed in our simulations (Fig. 4 and Fig. 5).To investigate the effect of nanobubble size and hydrate-nanobubble distance on hydrate dissociation, we simulate the hydrate dissociation process with nanobubbles of different sizes (HW-nb-L, HW-nb-M, and HW-nb-S) and at different hydrate-nanobubble distances (HW-nb-LC system).Simulation snapshots showing hydrate dissociation for the remaining HW-nb systems at different times are shown in the Supporting Information (Fig. S3).The hydrate clusters dissociate from the hydrate surface to the inner in a layer-by-layer manner.The variation in the number of methane molecules in the hydrate cluster with time, dissociation completion time, and the average dissociation rate of all four systems calculated using Equation (4), as well as those of HW system, are shown in Fig. 12.The presence of nanobubbles can enhance the hydrate dissociation process in different degrees at the same temperature and pressure, regardless of the size and distance of the pre-existing nanobubble (Fig. 12(a)).When the initial nanobubble is close to the hydrate cluster, it can significantly increase the average hydrate dissociation rate as shown in Fig. 12(b).Regarding the effect of the nanobubble size, the simulations do not indicate a clear trend.The hydrate cluster dissociates the fastest in the system with the medium-sized nanobubble, followed by the system with the large-sized nanobubble, and the dissociation is the slowest in the system with the small-sized nanobubble (Fig. 12).Our results show that the effect of the small nanobubble on the hydrate dissociation has the characteristic of a time lag.The variation in the number of methane molecules in hydrate cluster with time is almost the same in the HW-nb-S system with a small nanobubble and HW for the first 60 ns (marked in Fig. 12(a) by a circle), indicating that this smallsized nanobubble may have a minor effect on methane diffusion at first due to the large distance.As the released methane molecules get closer to the nanobubble, the effect of the nanobubble is enhanced.In the simulations with the large and the medium pre-existing nanobubbles, the effect of the nanobubble was observed in the initial dissociation stage.As shown in Fig. 4 and Fig. S3, both the nanobubble and hydrate c c Fig. 9. Methane density along the X-axis at different simulation times for HW (a) and HW-nb-L (b) systems at 292 K.The insets show methane concentration along the X-axis near the hydrate surface.A time interval of 1 ns is used to calculate the average value of the methane densities and concentration distributions.c x Fig. 10.Negative derivatives of the methane concentration in the liquid phase along the X-axis.The distance (x) is taken from the residual hydrate surface to the central position of the bulk phase and from the center of area 1 to the bulkphase center after the hydrate dissociated in the HW system.The calculation distance is from the hydrate surface to the pre-existing surface and from the center of area 1 to the bubble surface after the hydrate dissociated.The lines are third order polynomials.The coefficients of determination R 2 are equal to 0.9405 and 0.9488 for the HW and HW-nb-L systems, respectively.cluster in the liquid phase can move during hydrate dissociation due to intermolecular interactions.The distance between the residual hydrate and nanobubbles also change with time.For this reason, the mechanism behind the effect of the nanobubble size on hydrate dissociation is still unclear.It is necessary to study the dynamic changes of the nanobubble in the liquid phase during the hydrate dissociation process.
During hydrate dissociation, methane molecules on the hydrate surface are released, which diffuse into the liquid phase, resulting in a rapid increase in the methane concentration near the hydrate cluster.If there is a nanobubble in the liquid phase, mass transfer from the liquid phase to the gas bulk is driven by the concentration gradient.The variation in the number of methane molecules as a function of time in the gas phase for varying nanobubble radii in the four Hw-nb systems is shown in Fig. 13.Methane molecules in the nanobubble come from two sources: the hydrate cluster and the pre-existing nanobubble.As shown in Fig. 13, at first the number of gas molecules and nanobubble radius was decreased, which may be attributed to the fact that initially there is no methane in the liquid phase.Therefore, the mass transfer occurred from the nanobubble to the liquid bulk phase, and methane molecules were partially solvated into the liquid phase.The time required for this process varied with approximately 20 ns for the HW-nb-L (Fig. 13(a)) and the HW-nb-M (Fig. 13(b)) systems and approximately 40 ns for the HW-nb-S system (Fig. 13(c)).There is little time for the solvation process in the HW-nb-LC system (Fig. 13(d)).In this stage, the released methane molecules from hydrate rarely diffuse into the nanobubble, and therefore, the nanobubble had a negligible effect on the methane concentration around the hydrate cluster and not influenced the hydrate dissociation process.As an abundant amount of methane was released from the hydrate cluster, the direction of mass transfer changed and the number of methane molecules originated from the pre-existing nanobubbles maintained the balance (black line in Fig. 13(ad)).The released methane molecules from the hydrate cluster diffused into the nanobubble from the liquid phase due to the supersaturated state of the liquid (red line in Fig. 13(ad)), leading to an increase in the amount of methane in the nanobubble and the nanobubble size until the hydrate dissociation finished.As discussed above, large amounts of gas molecules diffusing into the nanobubble can decease the methane concentration in the liquid phase, which can increase the dissociation rate by enhancing the driving force.
As discussed above, the hydrate dissociation rate highly depends on the methane concentration around the residual hydrate cluster at the same temperature and pressure conditions.The mass transfer transits across the gas-liquid interface in the liquid phase.The mass flux of methane can be described by where J CH4 is the mass flux of methane molecules, K CH4 is the mass transfer coefficient, and C l and C eq are the methane concentrations in  water and in water at the equilibrium condition, respectively.For the dissociation system with a nanobubble, the decrease in the methane concentration in water was related to an increase in the mass rate across the nanobubble surface.
dm/dt = J CH4 ⋅A (9) where A is the surface area of the bubble.As shown in Fig. 13, the size of the nanobubble increases with the dissociation process, and A is not a constant.Therefore, we attempted to establish a relationship between the molecular flow rate and the dissociation rate for the four systems with the pre-existing nanobubble at 292 K.This is shown in Fig. 14.We feel that the effect of the nanobubble size and hydrate-nanobubble distance on the hydrate dissociation process is caused by the decreasing methane concentrations in the liquid phase.The hydrate dissociation rates are linearly related to the rate of mass transfer into the nanobubble in the surrounding liquid phase.

Conclusions
In this study, we investigated the effects of the presence of nanobubbles on hydrate dissociation kinetics in the liquid phase.A hydrate cluster was placed in a liquid bulk phase with and without a pre-existing nanobubble near the hydrate cluster.MD simulations were performed at different temperatures and for varying nanobubble sizes and hydratenanobubble distances.The simulation results showed that the hydrate dissociation process does not necessarily result in the formation of a new nanobubble because methane molecules can be dispersed in a methane supersaturated water solution.For higher temperatures, new nanobubbles were formed, even if there was a pre-existing bubble in the liquid bulk phase.The pre-existing nanobubble does not change the hydrate dissociation mechanism, regardless of the nanobubble size and distance.In all cases, the hydrate dissociates layer-by-layer in a shrinking core manner.By constructing a numerical model, the hydrate dissociation rate per surface area was calculated for each system, which showed that the dissociation rates of the hydrate clusters are enhanced rapidly along with the increase of the temperature.Temperature plays an important role in hydrate dissociation.In the presence of a preexisting nanobubble near the hydrate cluster, a methane concentration gradient was observed between the residual hydrate and the nanobubble surface.For mass transfer, the hydrate cluster acted as a methane source and the nanobubble acted as a gas reservoir.This concentration gradient of methane enhanced the methane transport process of released methane, thereby, decreasing the methane concentration in the liquid phase compared to the system without a nanobubble at the same  temperature and pressure.This phenomenon could increase the difference in the chemical potential between the solid and liquid phases and increase the dissociation driving force, thereby facilitating the hydrate dissociation process.The nanobubble does not decrease the activation energy of hydrate dissociation, but it can change the intrinsic decomposition rate of hydrate.Finally, we investigated the effects of the nanobubble size and distance on hydrate dissociation.The simulation results showed that the hydrate dissociation rates were linearly related to the rate of methane mass flow from the liquid to the nanobubble phase.
Formation of nanobubbles is a commonly occurring phenomenon during hydrate dissociation.Our investigation of the effects of nanobubbles on hydrate dissociation will help improving the understanding of the hydrate dissociation process.The existing nanobubble near the hydrate cluster could significantly facilitate the dissociation process at the same temperature and pressure.If large amount of nano or micro bubbles form during hydrate dissociation process, mass transfer can become more complex, and the interaction between these bubbles cannot be ignored, further studies are necessary to investigate; On the other hand, heat transfer can be impeded, resulting in a decrease in temperature, which may reduce the rate of hydrate dissociation [35].Although known physical and chemical technologies can control the formation of nanobubbles in water [78], heat losses should be taken into account and compensated for during hydrate dissociation.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.Initial configuration of the hydrate dissociation (a) without (HW) and (b) with (HW-nb) nanobubbles of different sizes (large, medium, and small).The initial HW-nb configuration was divided into four areas along the X-axis.In the HW-nb-L system, a cylindrical pre-existing nanobubble was placed in Area 3 with a 5.24-nm radius.In the HW-nb-M system, a medium-sized nanobubble with a 3.49-nm radius was present in area 3.In the HW-nb-S system, a smallsized nanobubble with a 2.44-nm radius was placed in area 3 in the bulk phase.In the HW-nb-LC system, we decreased the nanobubble-hydrate distance based on the configuration of the HW-nb-L system.Methane molecules in pre-existing nanobubble and hydrate are colored by cyan and purple, respectively.For the unit sI hydrate cell, the 512 cage is shown in blue and the 512 6 2 cage is shown in green.

Fig. 2 .
Fig. 2. Snapshots (a) and density map (b) of the HW-nb-L system at 2 ns.A time interval of 1 ns is used to calculate the average value of densities.

Fig. 3 .
Fig. 3.The potential energy of the system as a function of time in systems (a) without (HW) and (b) with (HW-nb-L) a nanobubble at different temperatures.(c) The potential energy as a function of time in HW and HW-nb systems for varying nanobubble sizes and nanobubble-hydrate distances at 292 K. Dashed lines represent the completion of the hydrate dissociation.

Fig. 4 .
Fig. 4. Typical simulation snapshots of the hydrate dissociation process for the HW and HW-nb-L systems at 292 K (first six figures a-f), and 312 K (next six pictures g-l).Hydrate cages are marked in red.The pre-existing nanobubble grew during the dissociation process, methane molecules inside the nanobubble are in cyan.A few new formed nanobubbles (methane molecules inside are marked in purple) were formed at a high dissociation temperature (312 K).

B
.Fang et al.

Fig. 6 .
Fig. 6.The number of methane molecules in the hydrate cluster in HW and HW-nb-L systems (a) at different temperatures.(b) Hydrate dissociation times of the two systems.

Fig. 7 .
Fig. 7. Dissociation rates of the per hydrate surface area in HW systems at 292 K.The initial and final hydrate dissociation parts in the HW system are indicated by the light red areas.Snapshots of the 5 12 6 2 cage dissociation process in the final hydrate dissociation stage are shown.

Fig. 8 .
Fig. 8. Arrhenius plots of ln k(T) as a function of 1000/T for the HW and HWnb-L systems.R 2 is the coefficient of determination in linear fitting.

Fig. 11 .
Fig. 11.Correlation between the dissociation rate and the radius of curvature (1/N 1/2 ) of the hydrate cluster in the two systems without (a) and with (b) the preexisting nanobubble at 292 K.The black line represents average values in each system.

Fig. 12 .
Fig. 12. Variations in the number of hydrate methane molecules in hydrate cluster with time (a).Hydrate dissociation time and the average dissociation rate (b) in the HW-nb system with a pre-existing nanobubbles of different sizes and at different distances in the liquid phase at 292 K.

Fig. 13 .
Fig. 13.Variation in the number of methane molecule in the nanobubble and nanobubble radius in the (a) HW-nb-L, (b) HW-nb-M, (c) HW-nb-S, and (d) HW-nb-LC systems.The methane molecules in the bubble can originate either from the hydrate phase or from the pre-existing nanobubble.

Fig. 14 .
Fig. 14.Average hydrate dissociation rate as a function of the molecular flow rate to the nanobubble from the liquid phase.