key: cord-0045741-0jfhilos authors: Zhou, Ji; Jiang, Wenbin; Lin, Mian; Ji, Lili; Cao, Gaohui title: Impact of Water on Methane Adsorption in Nanopores: A Hybrid GCMC-MD Simulation Study date: 2020-06-15 journal: Computational Science - ICCS 2020 DOI: 10.1007/978-3-030-50417-5_14 sha: 95d99eb076499a6d5b171728505ac67466c00086 doc_id: 45741 cord_uid: 0jfhilos Adsorbed methane is an important component of shale gas. Shale generally contains a certain amount of primary water, and isothermal adsorption experiments on wet samples show that water inhibits methane adsorption. Researches on methane adsorption mainly focus on the conditions of low pressure and water content. In this study, a hybrid GCMC-MD simulation method is proposed to study methane adsorption characteristics under high pressure and water content in pores of different sizes. This method can obtain the bulk pressure of the system while ensuring the simultaneous movement of methane and water molecules, and has high efficiency and reliability. It is found that the existence of water does not change the morphology of excess isotherm, and the relative decrease of adsorption capacity due to the existence of water is not sensitive to temperature. In ≤3 nm pores, water molecules form water clusters and partially occupy wall adsorption sites, and the adsorption amount decreases linearly with increasing water saturation. In the 5 nm wide pore with 40% water saturation, water films formed and methane adsorption is strongly suppressed. It is expected these findings could provide guidance for the evaluation of the amount of adsorbed methane with primary water. Methane is the main component of natural gas in low-porosity and low-permeability shale gas reservoir. Isothermal adsorption experiments reveal that the capacity of methane adsorption in shale is considerable, and the amount of adsorbed gas could be equivalent to that of free gas, accounts for 20%-85% of the total gas content in the USA and China [1] [2] [3] . In addition, adsorbed gas will influence the behavior of gas flow and diffusion [3] . Thus, it is very important for the evaluation of shale gas content and recovery to obtain a deeper understanding of methane adsorption characteristics. Currently, methane isothermal adsorption experiments are mainly conducted on dry shale by many departments of industry and researchers. It is found that: 1) the excess adsorption amount firstly increases and then decreases with the increase of pressure, and the isotherm could not be fitted by the standard Langmuir model; 2) the maximum adsorption amount decreases with the increase of temperature; 3) the maximum adsorption amount is proportional to the specific surface. Through high-resolution scanning electron microscopy (SEM) imaging and low-temperature gas adsorption analysis, abundant of nanopores with high specific surface area are found in shale [4] . Based on this finding, researchers simulated the adsorption of methane molecules in nanometer channels with side walls of different compositions, such as graphite (substitute of organic matter), different types of kerogen, quartz, illite and montmorillonite [5, 6] . A high-methane-density layer was observed near the wall and the excess loading can be determined, and the excess isotherms show a similar trend as that obtained from experiments. Such consistency delineates that the van der Waals force and physical adsorption play the leading role in methane adsorption in shale. Molecular simulation is proved to be able to reveal the microscopic mechanism of adsorption. Gas-bearing shale always contains a certain amount of primary water. In recent years, some researchers conducted isothermal adsorption experiments on shales with moisture, and found that the moisture could reduce the adsorption capacity of methane by more than 50% [7] [8] [9] . In recent years, some molecular simulations of competitive adsorption of methane and water have also been carried out. Billemont et al. [10] used grand canonical Monte Carlo (GCMC) to simulate the methane sorption with three different water contents in the graphite pore at 300 K, and observed that the preloaded water notably decreased the adsorption capacity of methane. Sui et al. [11] conducted GCMC simulations on methane adsorption in montmorillonite plates with and without Na + cations containing 0 mmol/cm 3 , 10 mmol/cm 3 and 20 mmol/cm 3 water at 298 K and pressures of 0-20 MPa, and the adsorbed amount of methane in montmorillonite pores with or without Na + both significantly decrease with the increase of water molecules. A similar result is also found by Zhao et al. [12] , who used GCMC to simulate the methane isothermal adsorption with 0-3% water content at temperatures of 298 K, 323 K, 348 K and pressure ranges from 0-20 MPa. In summary, molecular simulations of competitive adsorption of methane and water were primarily limited to relatively low pressure and water content, and did not cover the range of pressure and water saturation of shale reservoir deeper than 2000 m. GCMC is widely adopted in the simulation of adsorption isotherm as the pressure is given as an input parameter. However, in high-density systems, the probability of acceptance of random molecules operations is relatively low. When considering mixed sorbates, different molecules need to be operated separately to make them move, and the computational consumption will double. MD is capable to simulate the movement of different molecules simultaneously via solving the Newtonian equations directly with high parallelism and efficiency. The limited of MD lie in that the pressure near the wall is hard to calculate, and the pressure and density of bulk phase (external system connected with the simulation cell with the same chemical potential) are difficult to determine accurately. In this work, we proposed a hybrid GCMC-MD method to simulate competitive adsorption of methane and water in nanopores with water saturation up to 40% under pressure up to 50 MP. The rests of the paper are arranged as follows: the detail of molecular simulation method is introduced in the next section. Section 3 shows the validation of simulation method and models. The results and some discussions are addressed in Sect. 4. Finally, the conclusions of this study are drawn. K. Lin et al. [13] found that the adsorption characteristics of methane in the pores with parallel graphite walls were similar to those in shale. Considering the focus of this work is the effect of different pore sizes and water contents on the methane adsorption, it is more convenient to perform a single factor analysis in pores with parallel graphite walls to disclose internal mechanism. A rectangular simulation cell is built with periodicities in the x, y and z directions. Four graphene sheets separated by a distance of 0.335 nm form the pore wall to meet the minimum mirror criterion. The bond length of C-C is 0.142 nm and the angle is 120°. The dimension along the x-y surface is 9.116 nm  7.681 nm. The pore size H is determined by the length of vacuum in z direction. In this study, a united atom model is used to represent the methane molecule [5] . The water molecule is described as SPC/E model [14] . The bond length between the hydrogen atom and the oxygen atom is 0.1 nm, and the bond angle is 109.47°. The schematic representation of the slit-like graphite pore is presented in Fig. 1 . Two interactions included van der Waals force and Coulomb force are taken into account. The Lennard-Jones (L-J) potential model is used to describe van der Waals force: while the Coulomb force is represented by the following equation where r ij is the distance between atom i and j; q i and q j are the charges of atom i and j; e ij and r ij are the L-J well depth and the L-J zero-energy separation distance, respectively; r cut is the cutoff distance, in this study it is set to be 4r max , and r max is the maximum value of L-J zero-energy separation distances of all types of atoms. The long-range effectively-infinite Coulombic pairwise interactions within the cutoff distance are computed directly; the interactions outside this distance are computed in reciprocal space using the particle-particle particle-mesh solver (pppm) solver with an accuracy of 10 −4 . Table 1 lists the value of these parameters of the particles involved in the simulation. Lorentz-Berthelot mixing rule is utilized to compute the parameters of unlike particle's interaction. In this work, we propose a simulation method (GCMC-MD) in which GCMC and MD are performed alternately to simulate the adsorption of methane with a certain water saturation. GCMC operation is applied only on methane molecules to ensure the chemical potential of the system is consistent with the target value, and the bulk phase density is known. MD computation is applied on both methane and water molecules to simulate their movements based on potential parameters and Newtonian equations. This method combines the advantages of both GCMC and MD. Comparing with pure GCMC simulation, this model reduces GCMC operations by halves as the number of water molecules is predetermined and unchanged, therefore, it is more efficient. On the other side, this method overcomes the difficulty of determining the bulk density in MD simulation through GCMC operations on methane molecules. It is implemented using LAMMPS [15] . The procedure of GCMC-MD is described as follow: 1) Performing one step of GCMC operation after every 50 fs, 100 random exchanges (insertions or deletions) of methane molecules are attempted; 2) MD simulation is performed in every time step. the NVT ensemble (atom number N, volume V and temperature T are constant) is applied with a time step of 2 fs. The Nose/Hoover thermostat with a relaxation time of 200 fs is adopted to maintain the system temperature; 3) Repeating 1) and 2) until the total simulation time step reaches the preset value. The total computation time is set to be 1.0 ns in this work. For a given simulation condition with a prefix pressure, temperature, water saturation and pore width, the first 2/3 time steps are for equilibration and the data in the last 1/3 time steps are recorded for analysis. The GCMC-MD simulation flow chart is shown in Fig. 2 . In general, GCMC operations include exchanges (insertions or deletions) of atoms or molecules with an imaginary reservoir and Monte Carlo (MC) moves within the simulation cell. In this method, only exchanges of atoms or molecules are performed in the GCMC step, and the following MD steps start from a sample that is generated by a successful attempt. Simulations performed with MD algorithm in the NVT ensemble are equivalent to successful MC operations. Overall, the simulated samples belong to the grand canonical ensemble. The total number of methane molecules N tol at equilibrium state can be obtained by averaging the number of methane molecules in the time steps of data recording. The mass of methane molecules in bulk phase is computed by q free V free where q free g=ml ð Þis the bulk phase density of methane and V free ml ð Þ is the volume of space that methane can enter. Then the excess adsorption amount of methane per unit surface area q ex (ml/m 2 ) is computed following the equation. where V std is the volume of 1 mol methane under the standard conditions, ml; N A ¼ 6:022  10 23 mol À1 is the Avogadro constant; M CH 4 is the molar mass of methane, g/mol; S is the inner surface area of the simulation cell, m 2 . Considering the extremely weak adsorption capacity of helium gas, it is used to measure the dead volume in the adsorption experiment. Similarly, helium is used to replace methane during the simulation in order to obtain the number of helium molecules in the system. And then the V free can be calculated based on the following equation. We performed the bulk phase density simulations on methane, water and helium at temperature of 313 K with GCMC-MD method, respectively. The simulated densities are close to the values from the database of NIST (National Institute of Standards and Technology) (Fig. 3) , indicating that the force field parameters of methane, water and helium used in the model are reasonable. The efficiency of GCMC-MD is also investigated. A comparison of methane bulk phase density attained from GCMC simulation and GCMC-MD simulation to investigate the difference of efficiency. The density obtained from the GCMC-MD simulation is nearly identical to that of GCMC simulation under the same simulation conditions (the interaction parameters, the time step and the size of simulation system, etc.). While the running time of GCMC is about 30 times as large as that of GCMC-MD (Table 2 ) on the same computer with 4 core 4 thread (i5-4690k) CPU. [10] conducted methane isothermal adsorption experiments on the dry and wet F400 activated carbon with 90% carbon content whose dominant pore size is 0.8 nm. We perform GCMC-MD simulation of methane adsorption in the graphite nanopore (H = 0.8 nm) with the same water content. In order to make the simulation procedure consistent with the experimental steps, water molecules are injected into the pore at first, and MD simulation are conducted for 0.1 ns to make water molecules reach equilibrium state. Then the methane molecules are injected into the pore to participate competitive adsorption through GCMC-MD simulation. The experimental and simulated adsorption isotherms under dry and wet condition are in good agreement with each other as depicted in Fig. 4 , which indicates that the parameter setting of force field between different molecules is reasonable and the simulation process is correct and reliable. In a pore with a size of H under given temperature (T), pressure (P) and water saturation (S w ) conditions, the excess adsorption amount q ex S w ; P; T; H ð Þcan be obtained through GCMC-MD simulation. Single factor analysis is performed to investigate the influence of each factor. We simulate the methane adsorption in a 3 nm pore under 313 K, 25 MPa with S w = 0-40%. The results exhibited in Fig. 5 include snapshots of the equilibrium distribution of methane (in green) and water (in red) molecules in the slit-like graphite pore, and the number density profile of methane along the vertical direction. The number density of methane molecules has a peak near the wall surface, and peaks decrease with increasing water saturation. The excess adsorption amount decline linearly with the increase of water saturation (see Fig. 6 ). This trend is similar to the experimental result obtained by Hu et al. [8] . It can be seen from the molecular distribution in the system that the molecules occupy the wall surface in the form of water clusters, occupying the adsorption sites and causing the adsorption amount to decrease. With the increase of water saturation, the occupied adsorption sites gradually increase, so the adsorption capacity continues to decrease. Some water molecules exist around the center of the channel, resulting in zigzag fluctuations in the center of the density profile. Methane adsorption simulations in 3 nm pore under 313 K, 25 MPa, S w = 0-40% with pressure ranging from 2 MPa to 50 MPa are further conducted. The adsorption isotherms with different S w show the same trend (Fig. 7) i.e. the excess adsorption amount grow with the increase of pressure and then abruptly decreasing. It means that water saturation does not affect the supercritical adsorption characteristics of methane. Therefore, we select the maximum excess adsorption amount q ex;max S w ; T; H ð Þas the characteristic value which represents the adsorption capacity of the nanopore of a certain pore size at a certain temperature and water saturation. In the 3 nm pore, we further simulate the adsorption of methane with S w = 0-40%, P = 2-50 MPa at 313 K, 353 K, 393 K, and a series of excess adsorption isotherms are obtained and the maximum adsorption capacity q ex,max (0-40%, 313-393 K, 3 nm) can be determined. Let d S w ; T; H ð Þ¼q ex;max S w ; T; H ð Þ =q ex;max 0%; T; H ð Þbe the relative decrease in methane adsorption capacity with different water content. The variations of d S w ; T; 3 nm ð Þwith S w and T are shown in Fig. 8 . At different temperatures, the relative decrease d of adsorption capacity has the same trend with increasing water content. Under the same S w , the relative decrease d at different temperatures are relatively close. Because the water molecules are always liquid within the studied temperature and pressure ranges, the state of occupying the adsorption site is less affected by temperature, so the relative decrease in different temperatures is close. Methane adsorption at 25 MPa, 313 K, with 40% water saturation in pores of 1 nm, 3 nm, and 5 nm are simulated, respectively, to investigate the difference in pores with different size. The distribution of methane and water molecules in nanopores with different pore sizes and the density profile of methane along the vertical direction are shown in Fig. 9 . It can be seen that in the 1 nm and 3 nm pores, water partially occupies the wall surface in the form of water clusters, and methane still has a higher density area at the wall surface; in the 5 nm pore, water molecules are completely spread out on the wall, forming a water film completely covering the wall surface, blocking the adsorption of methane, the methane density near the wall surface is lower than the center. The reasons for these differences are analyzed from the perspective of intermolecular interaction forces. In the 1 nm and 3 nm pores, the water clusters are in contact with the upper and lower wall at the same time, which results in finite wall area being occupied. In the 5 nm pores, the upper and lower walls are far away, which makes the interactions between water molecules adsorbed on the upper and lower walls are too weak to form clusters. In order to study the characteristics of methane adsorption under the conditions of high pressure, high water content and a certain water saturation in pores with different sizes, we propose a hybrid GCMC-MD simulation method, and perform molecular simulations of methane adsorption under different temperature, pressure and water saturation conditions in pores with different pore sizes. The conclusions reached are as follows: (1) The hybrid simulation method can obtain the bulk pressure of the system while ensuring that the methane and water molecules move under the interaction force, which has high efficiency and can generate reliable results; (2) In 3 nm pores, water molecules form water clusters and partially occupy wall adsorption sites. As the water content increases, the adsorption amount decreases linearly. In the graphite pore with a width of 5 nm under 40% water saturation, water molecules completely cover the surface and form a water film, and the amount of methane adsorption is greatly reduced. (3) The water-containing conditions do not change the tendency of the adsorption isotherm to rise first and then decrease with increasing pressure. The relative decrease of methane adsorption capacity due to the existence of water is not sensitive to temperature. It is expected these findings could provide guidance for the evaluation of the amount of adsorbed methane with primary water. Fractured shale-gas systems Barnett shale gas-in-place volume including sorbed and free gas volume Geological features of Mesozoic lacustrine shale gas in south of Ordos Basin Mississippian Barnett Shale: Lithofacies and depositional setting of a deep-water shale-gas succession in the Fort Worth Basin Shale gas-in-place calculations part I: new pore-scale considerations Molecular dynamics investigation of conversion methods for excess adsorption amount of shale gas Shale gas potential of the lower Jurassic Gordondale member, northeastern British Columbia Influence of reservoir primary water on shale gas occurrence and flow capacity Experimental study on methane adsorption characteristics of hydrous shale An experimental and molecular simulation study of the adsorption of carbon dioxide and methane in nanoporous carbons in the presence of water Molecular simulation of shale gas adsorption and diffusion in clay nanopores Molecular simulation of methane adsorption on type II kerogen with the impact of water content Using graphene to simplify the adsorption of methane on shale in MD simulations A consistent empirical potential for water-protein interactions Fast parallel algorithms for short-range molecular dynamics