Solubility of Methane in Water: Some Useful Results for Hydrate Nucleation

Molecular insights into methane hydrate dissociation: Role of methane nanobubble formation

Understanding the underlying physics of natural gas hydrate dissociation is necessary for efficient CH4 extraction and in the exploration of potential additives in the chemical injection method. Silica being "sand" is already present inside the reservoir, making the silica nanoparticle a potential green additive. Here, molecular dynamics (MD) simulations have been performed to investigate the dissociation of the CH4 hydrate in the presence and absence of ∼1, ∼2, and ∼3 nm diameter hydrophilic silica nanoparticles at 100 bar and 310 K. We find that the formation of a CH4 nanobubble has a strong influence on the dissociation rate. After the initial hydrate dissociation, the rate of dissociation slows down till the formation of a CH4 nanobubble. We find the critical concentration and size limit to form the CH4 nanobubble to be ∼0.04 mole fraction of CH4 and ∼40 to 50 CH4 molecules, respectively. The solubility of CH4 and the chemical potential of H2O and CH4 are determined via Gibbs ensemble Monte Carlo simulations. The liquid phase chemical potential of both H2O and CH4 in the presence and absence of the nanoparticle is nearly the same, indicating that the effect of this additive will not be significant. While the formation of the hydration shell around the nanoparticle via hydrogen bonding confirms the strength of interactions between the water molecules and the nanoparticle in our MD simulations, the contact of the nanoparticle with the interface is infrequent, leading to no explicit effect of the nanoparticle on the dynamics of methane hydrate dissociation.

Dissociation line and driving force for nucleation of the nitrogen hydrate from computer simulation. II. Effect of multiple occupancy

In this work, we determine the dissociation line of the nitrogen (N2) hydrate by computer simulation using the TIP4P/Ice model for water and the TraPPE force field for N2. This work is the natural extension of Paper I, in which the dissociation temperature of the N2 hydrate has been obtained at 500, 1000, and 1500 bar [Algaba et al., J. Chem. Phys. 159, 224707 (2023)] using the solubility method and assuming single occupancy. We extend our previous study and determine the dissociation temperature of the N2 hydrate at different pressures, from 500 to 4500 bar, taking into account the single and double occupancy of the N2 molecules in the hydrate structure. We calculate the solubility of N2 in the aqueous solution as a function of temperature when it is in contact with a N2-rich liquid phase and when in contact with the hydrate phase with single and double occupancy via planar interfaces. Both curves intersect at a certain temperature that determines the dissociation temperature at a given pressure. We observe a negligible effect of occupancy on the dissociation temperature. Our findings are in very good agreement with the experimental data taken from the literature. We have also obtained the driving force for the nucleation of the hydrate as a function of temperature and occupancy at several pressures. As in the case of the dissociation line, the effect of occupancy on the driving force for nucleation is negligible. To the best of our knowledge, this is the first time that the effect of the occupancy on the driving force for nucleation of a hydrate that exhibits sII crystallographic structure is studied from computer simulation.

Accuracy limit of non-polarizable four-point water models: TIP4P/2005 vs OPC. Should water models reproduce the experimental dielectric constant?

The last generation of four center non-polarizable models of water can be divided into two groups: those reproducing the dielectric constant of water, as OPC, and those significantly underestimating its value, as TIP4P/2005. To evaluate the global performance of OPC and TIP4P/2005, we shall follow the test proposed by Vega and Abascal in 2011 evaluating about 40 properties to fairly address this comparison. The liquid-vapor and liquid-solid equilibria are computed, as well as the heat capacities, isothermal compressibilities, surface tensions, densities of different ice polymorphs, the density maximum, equations of state at high pressures, and transport properties. General aspects of the phase diagram are considered by comparing the ratios of different temperatures (namely, the temperature of maximum density, the melting temperature of hexagonal ice, and the critical temperature). The final scores are 7.2 for TIP4P/2005 and 6.3 for OPC. The results of this work strongly suggest that we have reached the limit of what can be achieved with non-polarizable models of water and that the attempt to reproduce the experimental dielectric constant deteriorates the global performance of the water force field. The reason is that the dielectric constant depends on two surfaces (potential energy and dipole moment surfaces), whereas in the absence of an electric field, all properties can be determined simply from just one surface (the potential energy surface). The consequences of the choice of the water model in the modeling of electrolytes in water are also discussed.

Three-phase equilibria of hydrates from computer simulation. II. Finite-size effects in the carbon dioxide hydrate

In this work, the effects of finite size on the determination of the three-phase coexistence temperature (T3) of the carbon dioxide (CO2) hydrate have been studied by molecular dynamic simulations and using the direct coexistence technique. According to this technique, the three phases involved (hydrate-aqueous solution-liquid CO2) are placed together in the same simulation box. By varying the number of molecules of each phase, it is possible to analyze the effect of simulation size and stoichiometry on the T3 determination. In this work, we have determined the T3 value at 8 different pressures (from 100 to 6000 bar) and using 6 different simulation boxes with different numbers of molecules and sizes. In two of these configurations, the ratio of the number of water and CO2 molecules in the aqueous solution and the liquid CO2 phase is the same as in the hydrate (stoichiometric configuration). In both stoichiometric configurations, the formation of a liquid drop of CO2 in the aqueous phase is observed. This drop, which has a cylindrical geometry, increases the amount of CO2 available in the aqueous solution and can in some cases lead to the crystallization of the hydrate at temperatures above T3, overestimating the T3 value obtained from direct coexistence simulations. The simulation results obtained for the CO2 hydrate confirm the sensitivity of T3 depending on the size and composition of the system, explaining the discrepancies observed in the original work by Míguez et al. [J. Chem Phys. 142, 124505 (2015)]. Non-stoichiometric configurations with larger unit cells show a convergence of T3 values, suggesting that finite-size effects for these system sizes, regardless of drop formation, can be safely neglected. The results obtained in this work highlight that the choice of a correct initial configuration is essential to accurately estimate the three-phase coexistence temperature of hydrates by direct coexistence simulations.

Three-phase equilibria of hydrates from computer simulation. I. Finite-size effects in the methane hydrate

Clathrate hydrates are vital in energy research and environmental applications. Understanding their stability is crucial for harnessing their potential. In this work, we employ direct coexistence simulations to study finite-size effects in the determination of the three-phase equilibrium temperature (T3) for methane hydrates. Two popular water models, TIP4P/Ice and TIP4P/2005, are employed, exploring various system sizes by varying the number of molecules in the hydrate, liquid, and gas phases. The results reveal that finite-size effects play a crucial role in determining T3. The study includes nine configurations with varying system sizes, demonstrating that smaller systems, particularly those leading to stoichiometric conditions and bubble formation, may yield inaccurate T3 values. The emergence of methane bubbles within the liquid phase, observed in smaller configurations, significantly influences the behavior of the system and can lead to erroneous temperature estimations. Our findings reveal finite-size effects on the calculation of T3 by direct coexistence simulations and clarify the system size convergence for both models, shedding light on discrepancies found in the literature. The results contribute to a deeper understanding of the phase equilibrium of gas hydrates and offer valuable information for future research in this field.

Three-phase equilibria of hydrates from computer simulation. III. Effect of dispersive interactions in the methane and carbon dioxide hydrates

In this work, the effect of the range of dispersive interactions in determining the three-phase coexistence line of the CO2 and CH4 hydrates has been studied. In particular, the temperature (T3) at which solid hydrate, water, and liquid CO2/gas CH4 coexist has been determined through molecular dynamics simulations using different cutoff values (from 0.9 to 1.6 nm) for dispersive interactions. The T3 of both hydrates has been determined using the direct coexistence simulation technique. Following this method, the three phases in equilibrium are put together in the same simulation box, the pressure is fixed, and simulations are performed at different temperatures T. If the hydrate melts, then T > T3. Conversely, if the hydrate grows, then T < T3. The effect of the cutoff distance on the dissociation temperature has been analyzed at three different pressures for CO2 hydrate: 100, 400, and 1000 bar. Then, we have changed the guest and studied the effect of the cutoff distance on the dissociation temperature of the CH4 hydrate at 400 bar. Moreover, the effect of long-range corrections for dispersive interactions has been analyzed by running simulations with homo- and inhomogeneous corrections and a cutoff value of 0.9 nm. The results obtained in this work highlight that the cutoff distance for the dispersive interactions affects the stability conditions of these hydrates. This effect is enhanced when the pressure is decreased, displacing the T3 about 2-4 K depending on the system and the pressure.

Modeling oceanic sedimentary methane hydrate growth through molecular dynamics simulation

The crystallization process of methane hydrates in a confined geometry resembling seabed porous silica sedimentary conditions has been studied using molecular dynamics simulations. With this objective in mind, a fully atomistic quartz silica slit pore has been designed, and the temperature stability of a methane hydrate crystalline seed in the presence of water and guest molecule methane has been analyzed. NaCl ion pairs have been added in different concentrations, simulating salinity conditions up to values higher than average oceanic conditions. The structure obtained when the hydrate crystallizes inside the pore is discussed, paying special attention to the presence of ionic doping inside the hydrate and the subsequent induced structural distortion. The shift in the hydrate stability conditions due to the increasing water salinity is discussed and compared with the case of unconfined hydrate, concluding that the influence of the confinement geometry and pore hydrophilicity produces a larger deviation in the confined hydrate phase equilibria.

Methane hydrate phase equilibrium considering dissolved methane concentrations and interfacial geometries from molecular simulations

Natural gas hydrates, mainly existing in permafrost and on the seabed, are expected to be a new energy source with great potential. The exploitation technology of natural gas hydrates is one of the main focuses of hydrate-related studies. In this study, a large-size liquid aqueous solution wrapping a methane hydrate system was established and molecular dynamics simulations were used to investigate the phase equilibrium conditions of methane hydrate at different methane concentrations and interfacial geometries. It is found that the methane concentration of a solution significantly affects the phase equilibrium of methane hydrates. Different methane concentrations at the same temperature and pressure can lead to hydrate formation or decomposition. At the same temperature and pressure, in a system reaching equilibrium, the size of spherical hydrate clusters is coupled to the solution concentration, which is proportional to the Laplace pressure at the solid-liquid interface. Lower solution concentrations reduce the phase equilibrium temperature of methane hydrates at the same pressure; as the concentration increases, the phase equilibrium temperature gradually approaches the actual phase equilibrium temperature. In addition, the interfacial geometry of hydrates affects the thermodynamic stability of hydrates. The spherical hydrate particles have the highest stability for the same volume. Through this study, we provide a stronger foundation to understand the principles driving hydrate formation/dissociation relevant to the exploitation of methane hydrates.

Understanding why constant energy or constant temperature may affect nucleation behavior in MD simulations: A study of gas hydrate nucleation

Molecular dynamics simulations have been widely used in exploring the nucleation behavior of many systems, including gas hydrates. Gas hydrates are ice-like solids in which gas molecules are trapped in water cages. During hydrate formation, a considerable amount of heat is released, and previous work has reported that the choice of temperature control scheme may affect the behavior of hydrate formation. The origins of this effect have remained an open question. To address this question, extensive NVE simulations and thermostatted (NPT and NVT) simulations with different temperature coupling strengths have been performed and compared for systems where a water nanodroplet is immersed in a H2S liquid. Detailed analysis of the hydrate structures and their mechanisms of formation has been carried out. Slower nucleation rates in NVE simulations in comparison to NPT simulations have been observed in agreement with previous studies. Probability distributions for various temperature measures along with their spatial distributions have been examined. Interestingly, a comparison of these temperature distributions reveals a small yet noticeable difference in the widths of the distributions for water. The somewhat reduced fluctuations in the temperature for the water species in the NVE simulations appear to be responsible for reducing the hydrate nucleation rate. We further conjecture that the NVE-impeded nucleation rate may be the result of the finite size of the surroundings (here the liquid H2S portion of the system). Additionally, a local spatial temperature gradient arising from the heat released during hydrate formation could not be detected.

Effects of freeze-thaw cycles on methanogenic hydrocarbon degradation: Experiment and modeling

Cold regions are warming much faster than the global average, resulting in more frequent and intense freeze-thaw cycles (FTCs) in soils. In hydrocarbon-contaminated soils, FTCs modify the biogeochemical and physical processes controlling petroleum hydrocarbon (PHC) biodegradation and the associated generation of methane (CH₄) and carbon dioxide (CO₂). Thus, understanding the effects of FTCs on the biodegradation of PHCs is critical for environmental risk assessment and the design of remediation strategies for contaminated soils in cold regions. In this study, we developed a diffusion-reaction model that accounts for the effects of FTCs on toluene biodegradation, including methanogenic biodegradation. The model is verified against data generated in a 215 day-long batch experiment with soil collected from a PHC contaminated site in Ontario, Canada. The fully saturated soil incubations with six different treatments were exposed to successive 4-week FTCs, with temperatures oscillating between -10 °C and +15 °C, under anoxic conditions to stimulate methanogenic biodegradation. We measured the headspace concentrations and ¹³C isotope compositions of CH₄ and CO₂ and analyzed the porewater for pH, acetate, dissolved organic and inorganic carbon, and toluene. The numerical model represents solute diffusion, volatilization, sorption, as well as a reaction network of 13 biogeochemical processes. The model successfully simulates the soil porewater and headspace concentration time series data by representing the temperature dependencies of microbial reaction and gas diffusion rates during FTCs. According to the model results, the observed increases in the headspace concentrations of CH₄ and CO₂ by 87% and 136%, respectively, following toluene addition are explained by toluene fermentation and subsequent methanogenesis reactions. The experiment and the numerical simulation show that methanogenic degradation is the primary toluene attenuation mechanism under the electron acceptor-limited conditions experienced by the soil samples, representing 74% of the attenuation, with sorption contributing to 11%, and evaporation contributing to 15%. Also, the model-predicted contribution of acetate-based methanogenesis to total produced CH₄ agrees with that derived from the ¹³C isotope data. The freezing-induced soil matrix organic carbon release is considered as an important process causing DOC increase following each freezing period according to the calculations of carbon balance and SUVA index. The simulation results of a no FTC scenario indicate that, in the absence of FTCs, CO₂ and CH₄ generation would decrease by 29% and 26%, respectively, and that toluene would be biodegraded 23% faster than in the FTC scenario. Because our modeling approach represents the dominant processes controlling PHC biodegradation and the associated CH₄ and CO₂ fluxes, it can be used to analyze the sensitivity of these processes to FTC frequency and duration driven by temperature fluctuations.

Effect of pressure on the carbon dioxide hydrate-water interfacial free energy along its dissociation line

We investigate the effect of pressure on the carbon dioxide (CO2) hydrate-water interfacial free energy along its dissociation line using advanced computer simulation techniques. In previous works, we have determined the interfacial energy of the hydrate at 400 bars using the TIP4P/Ice and TraPPE molecular models for water and CO2, respectively, in combination with two different extensions of the Mold Integration technique [J. Colloid Interface Sci. 623, 354 (2022) and J. Chem. Phys. 157, 134709 (2022)]. Results obtained from computer simulation, 29(2) and 30(2) mJ/m2, are found to be in excellent agreement with the only two measurements that exist in the literature, 28(6) mJ/m2 determined by Uchida et al. [J. Phys. Chem. B 106, 8202 (2002)] and 30(3) mJ/m2 determined by Anderson et al. [J. Phys. Chem. B 107, 3507 (2002)]. Since the experiments do not allow to obtain the variation of the interfacial energy along the dissociation line of the hydrate, we extend our previous studies to quantify the effect of pressure on the interfacial energy at different pressures. Our results suggest that there exists a correlation between the interfacial free energy values and the pressure, i.e., it decreases with the pressure between 100 and 1000 bars. We expect that the combination of reliable molecular models and advanced simulation techniques could help to improve our knowledge of the thermodynamic parameters that control the interfacial free energy of hydrates from a molecular perspective.

Solubility of carbon dioxide in water: Some useful results for hydrate nucleation

In this paper, the solubility of carbon dioxide (CO2) in water along the isobar of 400 bar is determined by computer simulations using the well-known TIP4P/Ice force field for water and the TraPPE model for CO2. In particular, the solubility of CO2 in water when in contact with the CO2 liquid phase and the solubility of CO2 in water when in contact with the hydrate have been determined. The solubility of CO2 in a liquid-liquid system decreases as the temperature increases. The solubility of CO2 in a hydrate-liquid system increases with temperature. The two curves intersect at a certain temperature that determines the dissociation temperature of the hydrate at 400 bar (T3). We compare the predictions with T3 obtained using the direct coexistence technique in a previous work. The results of both methods agree, and we suggest 290(2) K as the value of T3 for this system using the same cutoff distance for dispersive interactions. We also propose a novel and alternative route to evaluate the change in chemical potential for the formation of hydrates along the isobar. The new approach is based on the use of the solubility curve of CO2 when the aqueous solution is in contact with the hydrate phase. It considers rigorously the non-ideality of the aqueous solution of CO2, providing reliable values for the driving force for nucleation of hydrates in good agreement with other thermodynamic routes used. It is shown that the driving force for hydrate nucleation at 400 bar is larger for the methane hydrate than for the carbon dioxide hydrate when compared at the same supercooling. We have also analyzed and discussed the effect of the cutoff distance of dispersive interactions and the occupancy of CO2 on the driving force for nucleation of the hydrate.

Homogeneous nucleation rate of methane hydrate formation under experimental conditions from seeding simulations

In this work, we shall estimate via computer simulations the homogeneous nucleation rate for the methane hydrate at 400 bars for a supercooling of about 35 K. The TIP4P/ICE model and a Lennard-Jones center were used for water and methane, respectively. To estimate the nucleation rate, the seeding technique was employed. Clusters of the methane hydrate of different sizes were inserted into the aqueous phase of a two-phase gas-liquid equilibrium system at 260 K and 400 bars. Using these systems, we determined the size at which the cluster of the hydrate is critical (i.e., it has 50% probability of either growing or melting). Since nucleation rates estimated from the seeding technique are sensitive to the choice of the order parameter used to determine the size of the cluster of the solid, we considered several possibilities. We performed brute force simulations of an aqueous solution of methane in water in which the concentration of methane was several times higher than the equilibrium concentration (i.e., the solution was supersaturated). From brute force runs, we infer the value of the nucleation rate for this system rigorously. Subsequently, seeding runs were carried out for this system, and it was found that only two of the considered order parameters were able to reproduce the value of the nucleation rate obtained from brute force simulations. By using these two order parameters, we estimated the nucleation rate under experimental conditions (400 bars and 260 K) to be of the order of log₁₀ (J/(m³ s)) = -7(5).

Homogeneous nucleation of crystalline methane hydrate in molecular dynamics transition paths sampled under realistic conditions

Methane hydrates are important from a scientific and industrial perspective, and form by nucleation and growth from a supersaturated aqueous solution of methane. Molecular simulation is able to shed light on the process of homogeneous nucleation of hydrates, using straightforward molecular dynamics or rare event enhanced sampling techniques with atomistic and coarse grained force fields. In our previous work [Arjun, T. A. Berendsen, and P. G. Bolhuis, Proc. Natl. Acad. Sci. U. S. A. 116, 19305 (2019)], we performed transition path sampling (TPS) simulations using all atom force fields under moderate driving forces at high pressure, which enabled unbiased atomistic insight into the formation of methane hydrates. The supersaturation in these simulations was influenced by the Laplace pressure induced by the spherical gas reservoir. Here, we investigate the effect of removing this influence. Focusing on the supercooled, supersaturated regime to keep the system size tractable, our TPS simulations indicate that nuclei form amorphous structures below roughly 260 K and crystalline sI structures above 260 K. For these temperatures, the average transition path lengths are significantly longer than in our previous study, pushing the boundaries of what can be achieved with TPS. The temperature to observe a critical nucleus of certain size was roughly 20 K lower compared to a spherical reservoir due to the lower concentration of methane in the solution, yielding a reduced driving force. We analyze the TPS results using a model based on classical nucleation theory. The corresponding free energy barriers are estimated and found to be consistent with previous predictions, thus adding to the overall picture of the hydrate formation process.