# Molecular Simulation of CH4 Nanoscale Behavior and Enhanced Gas Recovery in Organic-Rich Shale

Authors: Yang Li; Lize Lu; Jingyi Zhu; Zhaozhong Yang; Jianhua Qu; Heng Xue; Jingyun Ouyang
Journal: Geofluids
 (2022)

---

## Abstract

Occurrence and transport are two important nanoscale behaviors in the exploitation of shale gas. Nanopores in a realistic shale organic matrix are composed of kerogen molecules, which will have a great impact on surface-gas interactions and gas nanoconfined behavior. Although there are previous studies, the physics of gas transport through shale systems remains ambiguous. In this work, cylindrical nanopore models representing different pore sizes and organic-rich shale were constructed. By applying the molecular dynamics simulation method, the occurrence characteristics and transport characteristics of CH4 in the nanopores of organic-rich shale were studied. At last, the process of the adsorbed CH4 displaced by CO2 and N2 in shale nanopores at the subsurface condition was explored. This work can provide a better understanding of gas nanoscale behavior in shale systems and assist the future design of the CO2 sequestration and enhanced gas recovery technique.


---

## Body

## 1. Introduction

Shale gas is one of the most promising alternative resources to conventional natural gas resources, and the development of shale gas is of strategic importance to energy security. However, the special reservoir conditions and seepage characteristics of shale gas determine the characteristics of rapid natural decline and low natural production capacity of gas reservoirs. Current shale gas reservoirs are mostly developed in a depletion drive, and the recovery rate is generally low. Unlike conventional oil and gas reservoirs, shale gas reservoirs have a complex pore structure and large size range. Typically, shale pores are classified into three categories according to the International Union of Pure and Applied Chemistry (IUPAC): micropores (<2 nm), mesopores (2-50 nm), and macropores (>50 nm) [1, 2]. Zhang et al. [3] studied the microscopic characteristics of shale pore structures, and they noted that 95.6% of shale pores are nanopores, making it difficult to visualize the microscopic flow behavior of CH4.Shale gas mainly consists of free gas in pores and fractures, adsorbed gas in clay surfaces and organic pores, and dissolved gas in the liquid phase [4–6]. The adsorbed gas content is considerable, which can account for 20%-85% of the total volume of shale gas [7]. Previous studies have shown that most of the adsorbed gas is mainly occurred in the organic matter nanopores of shale [8], and the adsorbed gas content is positively correlated with the organic carbon content of shale [9]. Kerogen is a dispersed organic matter in shale, although it only accounts for a small proportion of the shale volume, but due to the large specific surface area and strong affinity for gas molecules, kerogen becomes the main space for shale adsorbed gas [10].In recent years, the method of gas injection to displace CH4 is usually used in engineering to enhance the extraction efficiency [11]. Among them, CO2 injection in shale gas reservoirs not only can sequester CO2 and reduce the greenhouse effect but also can promote CH4 desorption and improve CH4 recovery rate, which has received wide attention in recent years. Therefore, the study of shale gas microflow behavior using kerogen as an adsorbent is an important guideline for enhanced shale gas recovery technology.Thanks to the advantages of molecular dynamics (MD) method at scale, they have been widely used in unconventional oil and gas extraction research, such as shale oil, shale gas, combustible ice, and coalbed methane (CBM) [12, 13]. It is because MD are based on classical mechanics and quantum mechanics method that can effectively reveal the link between structure and properties as well as understand the interactions between physicochemical systems. The pore structure of organic matter in shale is complex and dispersed in inorganic matter, which is difficult to be separated from shale while keeping its pore structure unbroken. The organic matter pores are mainly distributed at the nanoscale, and the pore throat channel size is comparable to the molecular mean free range, which makes it difficult for experimental methods to characterize the dynamic behavior of shale gas in these pore channels. Therefore, MD have become a powerful tool to study the structure and properties of shale at the nanoscale, and the MD method can have a greater advantage in studying the interaction between kerogen and fluids [14].In molecular simulation calculations, two main ways are currently used to study the interaction between shale and gases; one is to select kerogen molecules as a representative to replace the structure of shale organic matter; the other research method is to use graphene sheets or carbon nanotubes to replace the carbon skeleton in shale for fluid adsorption and transport properties. Yiannourakou et al. [15] on the basis of experimental techniques established the molecular fragment of mature kerogen, and the predicted CH4 adsorption isotherm of the obtained microporous material was in good agreement with the experimental data. Collell et al. [16] replaced the molecular structure of kerogen with aromatic hydrocarbon molecular groups and obtained the structure of shale organic matter with a density of 1.24 g/cm3 by MD method. Wang et al. [17] investigated the effect of graphite flakes on the adsorption and transport of CH4 molecules using the Monte Carlo method. The adsorption behavior of pores composed of graphite showed that pores with sizes in the range of 1-50 nm were the main pores for shale gas adsorption, and the adsorption showed a multilayer trend as the pressure increased subject to intermolecular interactions. Kazemi et al. [18, 19] reported that in slit-like organic nanopores of 2 nm, as the pressure gradient increased from 1.8 kPa·nm-1 to 17.65 kPa·nm-1, the shale gas flux increased from 4498 to 37244 kg·m-2·s-1. The contribution of adsorbed molecules could exceed 50%, and with increasing pore size, the adsorption of 10 nm slit organic pores decreased to less than 40%. Cao et al. [20] and Noorian et al. [21] similarly investigated the experimentally observed slip phenomenon through molecular simulations, showing that the slip length was related not only to the Knudsen number but also to the surface roughness. Ungerer et al. [10] proposed six different types of molecular models of low to high maturity keregon that exhibit thermodynamic properties consistent with trends in the literature. These structures are widely used in studies of competing adsorption processes of CO2/CH4 and transport characteristics of mixed fluids [22, 23]. Slit organic pores can also be constructed from these structures for analysis of fluid flow processes [24]. These results show that MD simulation is an effective method to compensate for the shortcomings of experimental methods in characterizing the dynamical processes in the structure of kerogen at nanoscale [10]. The effects of pore size distribution, pore throat size, pore roughness, and pore network connectivity on shale gas transportation can also be effectively investigated by MD [25–27].Through experimental analysis, it is found that the pore channels in shale organic matter are extremely irregular and the pore surface is very heterogeneous, and it has been found that the roughness of the pore surface has a large effect on both the flow rate and state of shale gas. However, so far, the studies have mainly focused on the fluid interaction with the simplified model, which differs significantly from the realistic pore model. In addition, the research on gas injection for CH4 displacement has mainly focused on CBM, and there are fewer studies for the special geological conditions of shale reservoir. There are also more factors affecting the effect of displacement, such as injection gas type, injection rate, and injection pressure and formation pressure, and the influence law of these main parameters on displacement is not very clear. Therefore, the fluid displacement and transport processes in shale organic matter nanopores are not well understood. The physical models constructed in this work embody the complex amorphous structure and strong inhomogeneous surface properties within the shale pores, which can be accurately characterized nanoconfined CH4 flow behavior in realistic shale organic matter. The aim of this paper is to investigate the flow characteristic laws of shale gas at the microscopic scale and to study the gas injection displacement thus helping to select reasonable parameters and improve the extraction efficiency, both of which are essential to guide the development of shale gas reservoirs.


## 2. Models and Methodology

### 2.1. Construction of Nanopore Structure Model

Although simulations can be performed based on simplified models equivalent to shale organic nanopores, simulations of gas microscopic behavior using realistic shale organic nanopore models are more necessary because the chemical composition, pore structure, and surface roughness of kerogen matrix have essential effect on the gas microscopic behavior in shale organic nanopores. Based on the investigation of the pore genesis and development characteristics of shale, it is found that the organic pores inside real shale are generally circlelike based on massive scanning electron microscopy microscopic images of shale samples [28–30].One of the nanopore structure model construction method is kerogen molecule randomly placed into a simulation cell with periodic boundaries, and then, a shale matrix model with microporous characteristics was established by the MD method under an isothermal-isobaric ensemble with a constant number of molecules. It should be noted that the shale organic nanopore structure model established by this method has complex basic properties, such as pore size distribution, morphological characteristics, topological characteristics, and surface roughness. All of these factors together lead to the corresponding gas microscopic behavior results. If the shale organic nanopore model is used directly to study the gas microscopic behavior under nanoconfined conditions, the specific effects of each factor need to be determined and quantified. To solve this problem, this paper generalizes the shale organic nanopores into capillary bundles, and the capillary bundle pore model assumes that the microscopic pores of porous media are composed of capillary bundles with different radius parallel to each other. It is a good solution to draw a single capillary and use it as single shale organic nanopore for the study of gas microscopic behavior.The most straightforward method to generate nanopores in organic solids is to remove the atoms in the target region [16]. This method results in accurate pore size and smooth pore surface after removing the solid atoms from the nanopore, but it also exposes a major problem in that the molecular topology and the charge distribution of the nanopore surface are disrupted. Collell et al. [16] proposed a method for inserting repulsive dummy atoms during the model construction, forming spherical nanopores with a radius of about 1 nm. This method compensates for the shortcomings of the atom removal method, however, it is not possible to form nanopores with other geometries such as slit pore or cylindrical pore, and it is also difficult to form pores with large radius. Based on Collell’s work, Zhou et al. [31] proposed an MD-based cutter atom pore generation method. This method can construct nanopores of arbitrary shape and size, and the constructed nanopores have a reasonable physical meaning of porous matrix. Therefore, this paper chooses to use the cutter atom method to construct the nanopore structure model of shale organic matter.As aforementioned, the kerogern molecule is carbon-based macromolecule representing shale organic matter, besides that, some representative kerogen molecular models had proposed by a recent review work [32]. Based on the different element ratios of kerogen with different maturation and burial histories, Bousige et al. [33] divided kerogen into four main types: type-I kerogen, type-II kerogen, type-III kerogen, and type-IV kerogen. As shale reservoirs are basically not generated by type-IV kerogen [34], therefore, in this work, we utilized three typical types (I, II, and III) of kerogen molecules to represent authentic shale organic matter with different thermal maturity [10]. The chemical formula for type-I kerogen is C251H385O13N7S3, type-II kerogen is C252H294O24N6S3, and type-III kerogen is C233H204O27N4; the molecular structures of kerogen molecules are present in Figure 1.Figure 1
The kerogen molecular model of different types. (a) The type-I kerogen C251H385O13N7S3; (b) the type-II kerogen C252H294O24N6S3; (c) the type-III kerogen C233H204O27N4. Color atom correspondence: gray for C, white for H, red for O, blue for N, and yellow for S.
(a)(b)(c)First, modeling the single-nanopore structure of shale organic matrix needs to determine the number of molecules required for different models. After determining the number of kerogen molecules, the kerogen molecules were randomly placed into the simulation cell with fixed cutter atoms, which were composed of carbon atoms forming a clipping region and located at the center of the simulation cell, as the initial state of the shale organic matrix skeleton. Then, a series of MD simulations were performed for the whole system to make the configuration relax. During the relaxation process, the cutter atoms remained fixed and the kerogen molecules were continuously compressed and cannot enter the clipping region, thus forming a nanopore surrounded by kerogen molecules.Finally, the cutter atoms were removed and the solid shale organic nanopores with the desired pore size and kerogen type were constructed by the above process. The pore size can be controlled by adjusting the volume of the clipping area formed by the cutter atoms. The process of constructing a nanopore with the target kerogen type and pore size is illustrated in Figure2.Figure 2
Schematic diagram of shale organic nanopore construction process.


### 2.2. Simulation Details

In this work, the MD of the microoccurrence behavior of CH4 was performed using the large-scale atomic/molecular massively parallel simulator (LAMMPS) from Sandia National Laboratories, USA [35]. The force field was chosen the all-atomic pcff+ force field developed by the Materials Design team[15], which uses the LJ 9-6 potential to describe dispersion-repulsion interactions and uses point charges to describe Coulomb interactions, classical bond stretching, bending, and torsion terms. Periodic boundary conditions were applied in the x, y, and z directions. The relaxation process was carried out by a series of heating and cooling MD simulations; the specific steps are listed in Table 1.Table 1
MD simulation steps for constructing shale organic nanopores.
EnsemblePressure (MPa)Temperature (K)Simulation time (ns)NVT-10003.0NPT6010006.0NPT607006.0NPT604006.0NVT-3383.0During the simulations, the cut-off distance set by the van der Waals force was 12 Å for short-range interactions, and long-range electrostatic interactions were calculated using the particle-particle-particle-mesh (PPPM) method with an accuracy set to1×10−4. The Wäldmann-Hagler mixing rule was implemented to describe the different atomic interactions [36]. The time step of the MD simulations was set to 1 fs for the microcanonical ensemble and 1×10−1 fs for the isothermal-isobaric ensemble, where constant temperature and constant pressure were achieved by coupling to the atomic velocity and unit dimension, respectively. The constant temperature was applied to the translational degrees of freedom of the atoms, and the Nosé-Hoover algorithm was used to control the temperature of the system [37]. The constant pressure was isotropically coupled to the entire volume of the simulated cell without changing its geometry [38]. To control the temperature and pressure, the damping coefficients were set to 100 and 1000 fs, respectively, and the equations of motion were calculated using the velocity Verlet method. After the relaxation process, the cutter atoms were removed and the configuration was followed by a short time MD equilibrium simulation under the canonical ensemble, which was used to reduce the possible stress concentration on the surface of the pore wall.After obtaining the structure model of the shale organic nanopore with the desired pore size and kerogen type, the gas microoccurrence behavior within the nanoconfined space will be simulated by applying the equilibrium molecular dynamics (EMD) simulation method. Using the cylindrical nanopore model as a flow channel, the dynamic characteristics of CH4 transport in shale nanopores were revealed by applying nonequilibrium molecular dynamics (NEMD) simulation method. At last, the MD simulation method was used to explore the process of the adsorbed CH4 displaced by CO2 and N2 in shale organic nanopores at the subsurface condition.


## 2.1. Construction of Nanopore Structure Model

Although simulations can be performed based on simplified models equivalent to shale organic nanopores, simulations of gas microscopic behavior using realistic shale organic nanopore models are more necessary because the chemical composition, pore structure, and surface roughness of kerogen matrix have essential effect on the gas microscopic behavior in shale organic nanopores. Based on the investigation of the pore genesis and development characteristics of shale, it is found that the organic pores inside real shale are generally circlelike based on massive scanning electron microscopy microscopic images of shale samples [28–30].One of the nanopore structure model construction method is kerogen molecule randomly placed into a simulation cell with periodic boundaries, and then, a shale matrix model with microporous characteristics was established by the MD method under an isothermal-isobaric ensemble with a constant number of molecules. It should be noted that the shale organic nanopore structure model established by this method has complex basic properties, such as pore size distribution, morphological characteristics, topological characteristics, and surface roughness. All of these factors together lead to the corresponding gas microscopic behavior results. If the shale organic nanopore model is used directly to study the gas microscopic behavior under nanoconfined conditions, the specific effects of each factor need to be determined and quantified. To solve this problem, this paper generalizes the shale organic nanopores into capillary bundles, and the capillary bundle pore model assumes that the microscopic pores of porous media are composed of capillary bundles with different radius parallel to each other. It is a good solution to draw a single capillary and use it as single shale organic nanopore for the study of gas microscopic behavior.The most straightforward method to generate nanopores in organic solids is to remove the atoms in the target region [16]. This method results in accurate pore size and smooth pore surface after removing the solid atoms from the nanopore, but it also exposes a major problem in that the molecular topology and the charge distribution of the nanopore surface are disrupted. Collell et al. [16] proposed a method for inserting repulsive dummy atoms during the model construction, forming spherical nanopores with a radius of about 1 nm. This method compensates for the shortcomings of the atom removal method, however, it is not possible to form nanopores with other geometries such as slit pore or cylindrical pore, and it is also difficult to form pores with large radius. Based on Collell’s work, Zhou et al. [31] proposed an MD-based cutter atom pore generation method. This method can construct nanopores of arbitrary shape and size, and the constructed nanopores have a reasonable physical meaning of porous matrix. Therefore, this paper chooses to use the cutter atom method to construct the nanopore structure model of shale organic matter.As aforementioned, the kerogern molecule is carbon-based macromolecule representing shale organic matter, besides that, some representative kerogen molecular models had proposed by a recent review work [32]. Based on the different element ratios of kerogen with different maturation and burial histories, Bousige et al. [33] divided kerogen into four main types: type-I kerogen, type-II kerogen, type-III kerogen, and type-IV kerogen. As shale reservoirs are basically not generated by type-IV kerogen [34], therefore, in this work, we utilized three typical types (I, II, and III) of kerogen molecules to represent authentic shale organic matter with different thermal maturity [10]. The chemical formula for type-I kerogen is C251H385O13N7S3, type-II kerogen is C252H294O24N6S3, and type-III kerogen is C233H204O27N4; the molecular structures of kerogen molecules are present in Figure 1.Figure 1
The kerogen molecular model of different types. (a) The type-I kerogen C251H385O13N7S3; (b) the type-II kerogen C252H294O24N6S3; (c) the type-III kerogen C233H204O27N4. Color atom correspondence: gray for C, white for H, red for O, blue for N, and yellow for S.
(a)(b)(c)First, modeling the single-nanopore structure of shale organic matrix needs to determine the number of molecules required for different models. After determining the number of kerogen molecules, the kerogen molecules were randomly placed into the simulation cell with fixed cutter atoms, which were composed of carbon atoms forming a clipping region and located at the center of the simulation cell, as the initial state of the shale organic matrix skeleton. Then, a series of MD simulations were performed for the whole system to make the configuration relax. During the relaxation process, the cutter atoms remained fixed and the kerogen molecules were continuously compressed and cannot enter the clipping region, thus forming a nanopore surrounded by kerogen molecules.Finally, the cutter atoms were removed and the solid shale organic nanopores with the desired pore size and kerogen type were constructed by the above process. The pore size can be controlled by adjusting the volume of the clipping area formed by the cutter atoms. The process of constructing a nanopore with the target kerogen type and pore size is illustrated in Figure2.Figure 2
Schematic diagram of shale organic nanopore construction process.


## 2.2. Simulation Details

In this work, the MD of the microoccurrence behavior of CH4 was performed using the large-scale atomic/molecular massively parallel simulator (LAMMPS) from Sandia National Laboratories, USA [35]. The force field was chosen the all-atomic pcff+ force field developed by the Materials Design team[15], which uses the LJ 9-6 potential to describe dispersion-repulsion interactions and uses point charges to describe Coulomb interactions, classical bond stretching, bending, and torsion terms. Periodic boundary conditions were applied in the x, y, and z directions. The relaxation process was carried out by a series of heating and cooling MD simulations; the specific steps are listed in Table 1.Table 1
MD simulation steps for constructing shale organic nanopores.
EnsemblePressure (MPa)Temperature (K)Simulation time (ns)NVT-10003.0NPT6010006.0NPT607006.0NPT604006.0NVT-3383.0During the simulations, the cut-off distance set by the van der Waals force was 12 Å for short-range interactions, and long-range electrostatic interactions were calculated using the particle-particle-particle-mesh (PPPM) method with an accuracy set to1×10−4. The Wäldmann-Hagler mixing rule was implemented to describe the different atomic interactions [36]. The time step of the MD simulations was set to 1 fs for the microcanonical ensemble and 1×10−1 fs for the isothermal-isobaric ensemble, where constant temperature and constant pressure were achieved by coupling to the atomic velocity and unit dimension, respectively. The constant temperature was applied to the translational degrees of freedom of the atoms, and the Nosé-Hoover algorithm was used to control the temperature of the system [37]. The constant pressure was isotropically coupled to the entire volume of the simulated cell without changing its geometry [38]. To control the temperature and pressure, the damping coefficients were set to 100 and 1000 fs, respectively, and the equations of motion were calculated using the velocity Verlet method. After the relaxation process, the cutter atoms were removed and the configuration was followed by a short time MD equilibrium simulation under the canonical ensemble, which was used to reduce the possible stress concentration on the surface of the pore wall.After obtaining the structure model of the shale organic nanopore with the desired pore size and kerogen type, the gas microoccurrence behavior within the nanoconfined space will be simulated by applying the equilibrium molecular dynamics (EMD) simulation method. Using the cylindrical nanopore model as a flow channel, the dynamic characteristics of CH4 transport in shale nanopores were revealed by applying nonequilibrium molecular dynamics (NEMD) simulation method. At last, the MD simulation method was used to explore the process of the adsorbed CH4 displaced by CO2 and N2 in shale organic nanopores at the subsurface condition.


## 3. Results and Discussion

### 3.1. CH4 Distribution

In this section, through Grand Canonical Monte Carlo (GCMC) and EMD simulation of CH4 molecules in shale organic nanopores, the distribution of CH4 in shale organic nanopores of different pore sizes and kerogen types under different pressure conditions were obtained. At the same time, the model was validated, and the CH4 distribution characteristics in the nanoconfined pores were analyzed. For shale reservoir, which is a multicomponent, heterogeneous porous medium, high temperature, and high pressure system [39–41], CH4 adsorption in in situ reservoir is particularly complex, which has become the bottleneck of shale gas occurrence mechanism research. To reveal the occurrence and transport characteristics of CH4, the influence of temperature and pressure at the subsurface condition will be evaluated in our simulation work.First, a certain number of CH4 molecules were put into the established nanopore space, and then, the system was simulated 3~10 ns under the canonical ensemble to make the system reach the equilibrium state, to obtain the stable pore pressure. Regarding the number of CH4 molecules put in, based on the ideal gas equation of state, the number of molecules in the shale organic nanopore space under the target pore pressure is preliminarily calculated. Then, it was adjusted based on the results of the EMD simulation, and finally, the final number of molecules by the average pressure of the CH4 molecules in the system was determined. The equilibrium pressure parameters obtained by the nanopore models of different pore sizes and kerogen types under different numbers of CH4 molecules are shown in Table 2.Table 2
The basic parameters of CH4 EMD simulation under nanoconfined conditions.
Kerogen typePore diameter (nm)CH4 molecular numberPressure (MPa)Type-I1759.812519.819050.2211010.221520.527051.345209.580019.7130052.1Type-II18010.211020.319049.8212010.022019.929051.945309.981019.6120050.5Type-III1859.212021.318051.8212010.621020.928050.945109.580019.9130050.2As can be seen from Table2, a total of 27 simulation schemes with different kerogen types, pressures, and pore sizes were set up to study the distribution state of CH4 in the nanopores. After EMD simulations, the CH4 density distribution of all schemes is displayed in the form of 2D cloud maps in Figures 3–5. In order to facilitate the comparison of the CH4 distribution states under different conditions, the maximum value of the color scale of all density distribution maps was set to 0.5 g/cm3.Figure 3
CH4 density distribution in 1 nm shale organic nanopore under different pressure conditions.Figure 4
CH4 density distribution in 2 nm shale organic nanopore under different pressure conditions.Figure 5
CH4 density distribution in 4 nm shale organic nanopore under different pressure conditions.It is obvious from the CH4 density distribution maps that the distribution of CH4 within the shale organic nanopores is heterogeneous and does not exactly follow the central symmetric, which is quite different from the previous simulation result using a simplified model based on graphene. Due to the strong chemical heterogeneity of the kerogen matrix, the nanopores have different chemical structures at different positions on the pore surface, which further indicates that the interaction strength of CH4 on the surface of the nanopores is different. Therefore, CH4 will accumulate on surfaces with high-energy adsorption sites. In contrast, on some weak-energy adsorption sites, the adsorption capacity is weak and the accumulation of CH4 will not occur. The complex surface properties are based on the inherent physical properties of kerogen molecules to construct shale organic nanopores, which is the necessary reason to explain the three-dimensional heterogeneous distribution of CH4 in shale organic nanopores. This is the critical point distinguishing the findings of this paper and other existing studies, and also, the proof of the claim that the physical model constructed in this paper shares the most similarity with the realistic shale organic matrix.In addition to observing the heterogeneous distribution of CH4 within the shale organic nanopores, it can also be clearly seen that there is high-density region near the nanopore wall and the bulk region inside the nanopores. CH4 not only adheres to the surface of shale organic nanopore, but also there is a small portion of CH4 which enters the kerogen matrix and stays there.As illustrated in the density distribution of CH4 inside shale organic nanopores, there is no clear boundary near the surface of the nanopores established by kerogen molecules compared to that established by graphene-based nanopores. It can also be analyzed from Figures 3–5 that the overall density of CH4 inside the nanopore increases with increasing pore pressure, both the near-surface region and the central bulk region. Notably, due to the inherent surface composition of the shale organic nanopore, distinct peak density sites can be found in specific nanopores, and the relative position of the peak density remains constant. In addition, only one main density peak can be observed in the shale organic nanopores, which implies that the adsorption behavior of CH4 within the shale organic nanopores is monolayer adsorption. It can be used as evidence of the application of monolayer adsorption theory to shale. Moreover, the thickness of the CH4 high-density region in the nanopores is almost the same regardless of the kerogen type, pore pressure, and pore size. This further suggests that the proportion of high-density region in small nanopores will be higher than that in large nanopores, which has an important influence on the average density in the nanopores.In order to further validate the model, the shale organic nanopores with the same pore size, different kerogen types, and pressure conditions were selected following, and the bulk-region density simulation results will be compared with the National Institute of Standards and Technology (NIST) database to prove the accuracy of the simulation method and model in this study.As can be observed in Figure6(a), the density of CH4 molecules in the bulk region (4 nm pore diameter) collected by EMD simulation can achieve excellent agreement with the NIST data. Based on the above content, it is shown that the molecular simulation conducted in this study successfully described the distribution characteristics of CH4 in the shale organic nanopores. It further shows that the force field parameter assignment on each particle in the simulation system was correct. Therefore, it can reflect the accuracy of the simulation method and model of this research to a certain extent.Figure 6
Comparison of CH4 density simulation results. (a) CH4 bulk density simulation results vs. NIST data, (b) density ratio vs. pressure with the same kerogen type, and (c) density ratio vs. pressure with the same pore size.
(a)(b)(c)From the above analysis and other recently simulation work [42, 43], it can be observed that the peak density in the high-density region is related to the pore size, while the bulk-region density is not affected by the pore size but mainly by the temperature and pressure to which the molecules within the nanopore were subjected. Therefore, the adsorption layer affects the average density within the nanopore. Moreover, the density distribution is a key factor affecting the fluid flow capacity from a microscopic point of view and the gas production performance of a gas well from a macroscopic point of view. Therefore, it is worthwhile to further analyze the density characteristics under the effect of adsorption within the nanoporeThe average density (ρave) inside the nanopore divides the bulk-region density (ρbulk) within the nanopore defined as the density ratio (Rdens), which is described as follows:
(1)Rdens=ρaveρbulk.The relationship between the density ratio and pressure of the type-II kerogen models with different pore sizes was analyzed, and the results are shown in Figure6(b).It can be seen from Figure6(b) that the density ratio varies roughly from 1.16 to 4.69 and under the same pressure condition increases as the pore size decreases. This can be explained by the fact that the thickness of the high-density adsorption region does not vary with the pore size, so that the proportion of high-density region in small-sized nanopores is higher than that in large-sized nanopores, thus leading to an increase in the average density with decreasing pore size. It can also be seen in Figure 6(b) that the density ratio of the same pore size decreases with increasing pressure and the magnitude of the change gradually decreases. Thus, the contribution of the high-density adsorption region dominates under low-pressure conditions, while this contribution gradually becomes weaker under high-pressure conditions.Considering that kerogen type also has an effect on the adsorption capacity, the variation relationship between density ratio and pressure was analyzed for the three kerogen type models with the same pore size (2 nm), and the results are shown in Figure6(c). It can be observed that the density ratio of the same kerogen type decreases with the increase of pressure. Besides, at a certain pressure, the density ratios show an increasing trend with the increase of maturity, in the order of high- to low-maturity kerogen. This phenomenon is related to the adsorption capacity of different kerogen type, and the high-maturity kerogen nanopores have a strong adsorption capacity, which means that at a certain pressure, compared with other kerogen type, the density peak formed by the adsorption layer will be higher.


### 3.2. CH4 Transport Characteristic

To study the flow pattern of CH4 within the nanopore, we used external force nonequilibrium molecular dynamics (EF-NEMD) simulations to deviate the occurrence system from equilibrium[44]. The specific procedure is that, based on the equilibrium state configuration obtained from EMD simulations, a force of constant magnitude is applied to the CH4 molecules inside the pore, which leads to NEMD simulations and makes the system enter the steady state. The force value was set so that the acceleration of each molecule was kept at 10-4 to 10-3 nm/ps2 to ensure a linear response. Due to the additional force on the gas, a certain pressure drop ΔP is generated in the flow direction, which can be calculated from
(2)ΔP=NG×FEAcross,whereNG,FE,andAcross denote the number of molecules applying external force, the external force, and the nanopore cross-sectional area.It should be noted that the fluid flow rate during the NEMD simulation refers to the center-of-mass motion of each fluid atom, and considering the degrees of freedom in the direction of fluid flow leads to an excessive calculated fluid temperature. Therefore, it is more reasonable to consider only the degrees of freedom perpendicular to the direction of the driving force when calculating the fluid temperature in this thermostatic way. Once the flow regime within the shale nanopores reaches a steady state, the fluid flux (Jz) can be obtained by measuring the average density ρ¯ and average flow velocity v¯ of the regime in the steady state, where JZ=ρ¯⋅v¯.In order to present the NEMD simulation results in detail and to further analyze the transport characteristics of CH4 molecules within the nanopores, the velocity distribution of CH4 molecules in shale organic nanopores under different conditions was analyzed below using control variates method. Figure 7 show the results of the velocity distribution of CH4 transport within type-I, type-II, and type-III kerogen nanopores with different pore pressures (10, 20, and 50 MPa), pore size (1, 2, and 4 nm), pressure drop (0.25, 0.75, and 1.25 MPa), and ambient temperature (308, 338, and 368 K). In addition, in order to control the single variable in different tests, the physical parameters of the benchmark one are 20 MPa (pore pressure), 2 nm (pore size), 0.75 MPa (pressure drop), and 338 K (ambient temperature).Figure 7
Velocity distribution of CH4 in shale organic nanopores under different conditions. Type I: (a−d), Type II: (e−h), Type III: (i−l). Pore pressure (a, e, and i), pressure drop (b, f, and j), pore size (c, g, and k), and ambient temperature (d, h, and l).
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)It can be seen that the variation of velocity distribution with different tests is essentially the same within the shale nanopore model for different kerogen types. For the impact of pore pressure (Figures7(a), 7(e), and 7(i)), when the pressure was 50 MPa, the shape of the velocity distribution curves is a parabolic shape with nonslip boundaries. As the pore pressure decreases, the velocity profiles within the different shale nanopores show a common feature, the generation of a slip boundary. When the pressure was 20 MPa, the shape of the velocity distribution curve in the nanopore is close to the classical parabolic shape and the velocity peak at the pore wall disappears, but the velocity inflection points in the adsorption layer region and the central region of the pore can still be seen. Since different shale organic nanopores at different pore pressures have similar velocity distribution characteristics, the velocities at different pore pressures were normalized using type-II kerogen shale organic nanopore as an example (Figure 8(a)), and the normalized velocity values were obtained based on the average velocities at different pressures within the pore.Figure 8
Normalized velocity and CH4 density distribution in the nanopores of type-II kerogen under different pore pressures.
(a)(b)As can be seen from Figure8(a), the velocity distribution within the shale nanopores can be divided into three regions. The first region (Np-wall) is adjacent to the shale nanopore wall, where gas atoms are strongly repelled by the solid wall and almost no atoms or velocities appear; the second region (Np-adsorbed) is around the adsorption layer where complex phenomena exist, the slip velocity decreases with increasing pore pressure and disappears at high pressures, and the velocity values are reduced at the potential well where density peaks exist compared to the other regions; the third region (Np-bulk) is the central bulk phase region of the nanopore, where the transport behavior mainly follows the classical parabolic pattern, but with different curvatures due to the rarefaction effect in the nanopores[45].It can be concluded that the CH4 molecules in the adsorption layer exhibit completely different transport characteristics compared to the bulk phase region in the center of the nanopore. This is mainly due to the different density distributions of the adsorbed layer and the bulk phase region within the shale nanopores leading to their different velocity distributions. As an example to illustrate the results of the density distribution curve of type-II kerogen (Figure 8(b)), when the pressure was low (10 MPa), the low peak at the nanopore wall in the density profile indicates that the adsorption and nanospatial confinement of the nanopore wall were weak, and fewer CH4 molecules were adsorbed on the nanopore surface, so it corresponds to Figure 8(a), which the transport velocity is higher near the nanopore wall. As the pore pressure increases, the density of CH4 molecules increases within the nanopore, and the difference in density between the Np-adsorbed region and the Np-bulk region decreases. When the pressure reached 50 MPa, due to the strong confinement effect in the shale nanopores under high pressure conditions, the velocity profile showing a nonslip boundary feature. That can be elaborated as the pressure increases; the main kinetic behavior of CH4 molecules in shale nanopores changes from collisions between gas molecules and nanopore wall to gas intermolecular collisions, which leads to significant differences in gas transport characteristics at different pressures.Figures7(b), 7(f), and 7(j) show the results of the velocity distribution under the influence of different pressure drop. It can be seen that the velocity distribution within the shale nanopore is positively correlated with the pressure drop. The approximate parabolic shape of the simulated CH4 velocity profile obtained in the center of the nanopore again indicates the presence of viscous flow during the transport of supercritical state CH4 within the shale nanopore. Furthermore, by observing closely at the different velocity profile values, especially in the Np-bulk region, it can be found a near-linear relationship between pressure drop and the velocity values. When the velocity values about 0.17 Å/ps at 1.25 MPa pressure drop are 1.7 times larger than ~0.1 Å/ps at 0.75 MPa pressure drop and are 5 times larger than ~0.03 Å/ps at 0.25 MPa pressure drop.Moreover, the results of the velocity distribution of CH4 molecules under the influence of different pore sizes are analyzed in Figures 7(c), 7(g), and 7(k). It can be clearly seen that the velocity profile of CH4 molecules within the shale nanopores is highly dependent on the pore size, with higher values of flow velocity for larger pore sizes. And further comparative analysis will reveal that the velocity profile shows a nonlinear relationship with pore size, as exemplified by the type-II kerogen (Figure 7(g)), where the maximum velocity in the central region of the nanopore (~0.27 Å/ps) is about 2.7 times higher for the 4 nm pore size condition than for the 2 nm pore size case (~0.1 Å/ps) and 6.75 times (~0.04 Å/ps) for the 1 nm pore size case. This illustrates the large effect of pore size on the gas flow capacity within shale nanopores, where several-fold differences in pore size can lead to flow rate differences of almost more than an order of magnitude. This is consistent with continuous flow theory (Poiseuille’s law), where the gas flow rate within the channel is proportional to the square of the channel size and the underlying mechanisms were revealed by previous research results [46, 47].Finally, the simulation results of CH4 molecular velocity distribution within the shale nanopores at different ambient temperatures were compared. From the velocity distribution of CH4 molecules within shale nanopores in Figures 7(d), 7(h), and 7(l), for the effect of ambient temperature, it can be found that the difference between velocity distributions under different ambient temperature conditions is small, indicating that the increase in temperature mainly accelerates the thermal movement of molecules and has little enhancement for the velocity in the transport direction. Therefore, the transport capacity of the gas can be roughly considered to be the same under different shale reservoir temperature conditions.In summary, the above analysis of the transport simulation results illustrates that the transport pattern of CH4 in shale nanopores is complex and different from the behavior in homogeneous carbon nanotubes or graphene nanochannels. The nonhomogeneous nanopore walls result in velocity profiles that are not symmetrically distributed, and the chemical structures of different kerogen types alter the interaction between CH4 and shale nanopore surface, which affects the molecular transport velocity. Conventional understanding ignores the shale gas transport pattern with adsorption, which may underestimate the gas permeability in shale reservoirs. In a future research work, we will compare the MD results with the classical theoretical model values and then develop a new model for the apparent permeability of gas in shale organic nanopore, trying to applicate our simulation work to real shale matrices with multiscale complex pore networks.


### 3.3. Displacement of the Adsorbed CH4

To improve the recovery of adsorbed CH4, gas injection displacement is an efficient method. According to the research investigation, the commonly chosen injection gas in enhanced gas recovery engineering is CO2 and N2 [48–50]. Compared to CH4, CO2 has a higher adsorption capacity and shale matrix has strong surface heterogeneity, in which the oxygen-containing functional groups on the surface of shale nanopore enhance the displacement of CO2 to CH4. Not only that, CO2 can be sequestered in the shale formation by means of gas injection. In addition to CO2, N2 is also a common choice for gas injection. According to Dalton’s partial pressure law, with the injection of N2, the partial pressure of CH4 will decrease. Besides, N2 has many advantages, such as being abundant, readily available, and nontoxic and noncorrosive.This section focuses on the displacement of adsorbed CH4 by the injection gas, which is still implemented using a MD method and with the help of the LAMMPS open source program. The model schematic diagram is shown in Figure 9; the injection gas is either CO2 or N2. In the initial state, both boxes were not connected to the shale organic nanopores. The injection gas was equilibrated in the left box while the CH4 was adsorbed in the nanopore. In this section, the selected working condition is that the injection gas pressure is 12 MPa while the adsorbed CH4 pressure is 10 MPa. In addition, the temperature of the whole system was controlled using a Nosé-Hoover thermostat at 338 K. Subsequently, create the entrance between the boxes and nanopore, and dig a hole (remove atoms) in the wall between the nanopore and the boxes on both sides, which is the size of the cross section of the nanopore. Then, the injected gas entered the nanopore for the displacement under the pressure difference. In subsequent simulations, the content of each component of the output gas in the right-hand box was counted to analyze the displacement mechanism and process.Figure 9
Molecular model of CH4 displacement by injection gas. The gray wall at the left-hand side represents piston to maintain the injection pressure in box. The injection gas is either CO2 or N2 and in the schematic diagram is CO2.The process of CH4 displacement by CO2 and N2 was studied comparatively. Figure 10 gives a graph of the variation in the concentration of each component of the output gas after the displacement of CH4 by CO2 and N2. When CO2 was injected, the mole fraction change curve shows a trend of long breakthrough time and sharp front, as shown in Figure 10(a). Since CO2 tends to adsorb on the shale nanopore surface, the amount of CO2 in the bulk state is less and its passage through the nanopore was slower. When the injected gas comes in contact with the adsorbed CH4, a front was formed. When CO2 breaks through the nanopore space, the concentration of output CH4 rapidly decreases to zero. This also indicates a sharper CO2 front within the nanopore space. On the contrary, when N2 was injected, the mole fraction change curve shows a trend of short breakthrough time and a flat front, as shown in Figure 10(b). Since N2 has a weaker adsorption capacity, it can pass through the nanopore faster and accelerate CH4 desorption earlier. Even though N2 has already broken through the nanopore, CH4 continues to be produced. When enough N2 was injected, the concentration of output CH4 gradually tends to zero. These findings suggest that both CO2 and N2 can efficiently displace the adsorbed CH4. The above simulation results are consistent with the experimental results of related literature, which once again validates the accuracy of the simulation method.Figure 10
Displacement gas composition recorded from the MD as a function of displacement time. (a) CO2 displacing CH4; (b) N2 displacing CH4.
(a)(b)


## 3.1. CH4 Distribution

In this section, through Grand Canonical Monte Carlo (GCMC) and EMD simulation of CH4 molecules in shale organic nanopores, the distribution of CH4 in shale organic nanopores of different pore sizes and kerogen types under different pressure conditions were obtained. At the same time, the model was validated, and the CH4 distribution characteristics in the nanoconfined pores were analyzed. For shale reservoir, which is a multicomponent, heterogeneous porous medium, high temperature, and high pressure system [39–41], CH4 adsorption in in situ reservoir is particularly complex, which has become the bottleneck of shale gas occurrence mechanism research. To reveal the occurrence and transport characteristics of CH4, the influence of temperature and pressure at the subsurface condition will be evaluated in our simulation work.First, a certain number of CH4 molecules were put into the established nanopore space, and then, the system was simulated 3~10 ns under the canonical ensemble to make the system reach the equilibrium state, to obtain the stable pore pressure. Regarding the number of CH4 molecules put in, based on the ideal gas equation of state, the number of molecules in the shale organic nanopore space under the target pore pressure is preliminarily calculated. Then, it was adjusted based on the results of the EMD simulation, and finally, the final number of molecules by the average pressure of the CH4 molecules in the system was determined. The equilibrium pressure parameters obtained by the nanopore models of different pore sizes and kerogen types under different numbers of CH4 molecules are shown in Table 2.Table 2
The basic parameters of CH4 EMD simulation under nanoconfined conditions.
Kerogen typePore diameter (nm)CH4 molecular numberPressure (MPa)Type-I1759.812519.819050.2211010.221520.527051.345209.580019.7130052.1Type-II18010.211020.319049.8212010.022019.929051.945309.981019.6120050.5Type-III1859.212021.318051.8212010.621020.928050.945109.580019.9130050.2As can be seen from Table2, a total of 27 simulation schemes with different kerogen types, pressures, and pore sizes were set up to study the distribution state of CH4 in the nanopores. After EMD simulations, the CH4 density distribution of all schemes is displayed in the form of 2D cloud maps in Figures 3–5. In order to facilitate the comparison of the CH4 distribution states under different conditions, the maximum value of the color scale of all density distribution maps was set to 0.5 g/cm3.Figure 3
CH4 density distribution in 1 nm shale organic nanopore under different pressure conditions.Figure 4
CH4 density distribution in 2 nm shale organic nanopore under different pressure conditions.Figure 5
CH4 density distribution in 4 nm shale organic nanopore under different pressure conditions.It is obvious from the CH4 density distribution maps that the distribution of CH4 within the shale organic nanopores is heterogeneous and does not exactly follow the central symmetric, which is quite different from the previous simulation result using a simplified model based on graphene. Due to the strong chemical heterogeneity of the kerogen matrix, the nanopores have different chemical structures at different positions on the pore surface, which further indicates that the interaction strength of CH4 on the surface of the nanopores is different. Therefore, CH4 will accumulate on surfaces with high-energy adsorption sites. In contrast, on some weak-energy adsorption sites, the adsorption capacity is weak and the accumulation of CH4 will not occur. The complex surface properties are based on the inherent physical properties of kerogen molecules to construct shale organic nanopores, which is the necessary reason to explain the three-dimensional heterogeneous distribution of CH4 in shale organic nanopores. This is the critical point distinguishing the findings of this paper and other existing studies, and also, the proof of the claim that the physical model constructed in this paper shares the most similarity with the realistic shale organic matrix.In addition to observing the heterogeneous distribution of CH4 within the shale organic nanopores, it can also be clearly seen that there is high-density region near the nanopore wall and the bulk region inside the nanopores. CH4 not only adheres to the surface of shale organic nanopore, but also there is a small portion of CH4 which enters the kerogen matrix and stays there.As illustrated in the density distribution of CH4 inside shale organic nanopores, there is no clear boundary near the surface of the nanopores established by kerogen molecules compared to that established by graphene-based nanopores. It can also be analyzed from Figures 3–5 that the overall density of CH4 inside the nanopore increases with increasing pore pressure, both the near-surface region and the central bulk region. Notably, due to the inherent surface composition of the shale organic nanopore, distinct peak density sites can be found in specific nanopores, and the relative position of the peak density remains constant. In addition, only one main density peak can be observed in the shale organic nanopores, which implies that the adsorption behavior of CH4 within the shale organic nanopores is monolayer adsorption. It can be used as evidence of the application of monolayer adsorption theory to shale. Moreover, the thickness of the CH4 high-density region in the nanopores is almost the same regardless of the kerogen type, pore pressure, and pore size. This further suggests that the proportion of high-density region in small nanopores will be higher than that in large nanopores, which has an important influence on the average density in the nanopores.In order to further validate the model, the shale organic nanopores with the same pore size, different kerogen types, and pressure conditions were selected following, and the bulk-region density simulation results will be compared with the National Institute of Standards and Technology (NIST) database to prove the accuracy of the simulation method and model in this study.As can be observed in Figure6(a), the density of CH4 molecules in the bulk region (4 nm pore diameter) collected by EMD simulation can achieve excellent agreement with the NIST data. Based on the above content, it is shown that the molecular simulation conducted in this study successfully described the distribution characteristics of CH4 in the shale organic nanopores. It further shows that the force field parameter assignment on each particle in the simulation system was correct. Therefore, it can reflect the accuracy of the simulation method and model of this research to a certain extent.Figure 6
Comparison of CH4 density simulation results. (a) CH4 bulk density simulation results vs. NIST data, (b) density ratio vs. pressure with the same kerogen type, and (c) density ratio vs. pressure with the same pore size.
(a)(b)(c)From the above analysis and other recently simulation work [42, 43], it can be observed that the peak density in the high-density region is related to the pore size, while the bulk-region density is not affected by the pore size but mainly by the temperature and pressure to which the molecules within the nanopore were subjected. Therefore, the adsorption layer affects the average density within the nanopore. Moreover, the density distribution is a key factor affecting the fluid flow capacity from a microscopic point of view and the gas production performance of a gas well from a macroscopic point of view. Therefore, it is worthwhile to further analyze the density characteristics under the effect of adsorption within the nanoporeThe average density (ρave) inside the nanopore divides the bulk-region density (ρbulk) within the nanopore defined as the density ratio (Rdens), which is described as follows:
(1)Rdens=ρaveρbulk.The relationship between the density ratio and pressure of the type-II kerogen models with different pore sizes was analyzed, and the results are shown in Figure6(b).It can be seen from Figure6(b) that the density ratio varies roughly from 1.16 to 4.69 and under the same pressure condition increases as the pore size decreases. This can be explained by the fact that the thickness of the high-density adsorption region does not vary with the pore size, so that the proportion of high-density region in small-sized nanopores is higher than that in large-sized nanopores, thus leading to an increase in the average density with decreasing pore size. It can also be seen in Figure 6(b) that the density ratio of the same pore size decreases with increasing pressure and the magnitude of the change gradually decreases. Thus, the contribution of the high-density adsorption region dominates under low-pressure conditions, while this contribution gradually becomes weaker under high-pressure conditions.Considering that kerogen type also has an effect on the adsorption capacity, the variation relationship between density ratio and pressure was analyzed for the three kerogen type models with the same pore size (2 nm), and the results are shown in Figure6(c). It can be observed that the density ratio of the same kerogen type decreases with the increase of pressure. Besides, at a certain pressure, the density ratios show an increasing trend with the increase of maturity, in the order of high- to low-maturity kerogen. This phenomenon is related to the adsorption capacity of different kerogen type, and the high-maturity kerogen nanopores have a strong adsorption capacity, which means that at a certain pressure, compared with other kerogen type, the density peak formed by the adsorption layer will be higher.


## 3.2. CH4 Transport Characteristic

To study the flow pattern of CH4 within the nanopore, we used external force nonequilibrium molecular dynamics (EF-NEMD) simulations to deviate the occurrence system from equilibrium[44]. The specific procedure is that, based on the equilibrium state configuration obtained from EMD simulations, a force of constant magnitude is applied to the CH4 molecules inside the pore, which leads to NEMD simulations and makes the system enter the steady state. The force value was set so that the acceleration of each molecule was kept at 10-4 to 10-3 nm/ps2 to ensure a linear response. Due to the additional force on the gas, a certain pressure drop ΔP is generated in the flow direction, which can be calculated from
(2)ΔP=NG×FEAcross,whereNG,FE,andAcross denote the number of molecules applying external force, the external force, and the nanopore cross-sectional area.It should be noted that the fluid flow rate during the NEMD simulation refers to the center-of-mass motion of each fluid atom, and considering the degrees of freedom in the direction of fluid flow leads to an excessive calculated fluid temperature. Therefore, it is more reasonable to consider only the degrees of freedom perpendicular to the direction of the driving force when calculating the fluid temperature in this thermostatic way. Once the flow regime within the shale nanopores reaches a steady state, the fluid flux (Jz) can be obtained by measuring the average density ρ¯ and average flow velocity v¯ of the regime in the steady state, where JZ=ρ¯⋅v¯.In order to present the NEMD simulation results in detail and to further analyze the transport characteristics of CH4 molecules within the nanopores, the velocity distribution of CH4 molecules in shale organic nanopores under different conditions was analyzed below using control variates method. Figure 7 show the results of the velocity distribution of CH4 transport within type-I, type-II, and type-III kerogen nanopores with different pore pressures (10, 20, and 50 MPa), pore size (1, 2, and 4 nm), pressure drop (0.25, 0.75, and 1.25 MPa), and ambient temperature (308, 338, and 368 K). In addition, in order to control the single variable in different tests, the physical parameters of the benchmark one are 20 MPa (pore pressure), 2 nm (pore size), 0.75 MPa (pressure drop), and 338 K (ambient temperature).Figure 7
Velocity distribution of CH4 in shale organic nanopores under different conditions. Type I: (a−d), Type II: (e−h), Type III: (i−l). Pore pressure (a, e, and i), pressure drop (b, f, and j), pore size (c, g, and k), and ambient temperature (d, h, and l).
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)It can be seen that the variation of velocity distribution with different tests is essentially the same within the shale nanopore model for different kerogen types. For the impact of pore pressure (Figures7(a), 7(e), and 7(i)), when the pressure was 50 MPa, the shape of the velocity distribution curves is a parabolic shape with nonslip boundaries. As the pore pressure decreases, the velocity profiles within the different shale nanopores show a common feature, the generation of a slip boundary. When the pressure was 20 MPa, the shape of the velocity distribution curve in the nanopore is close to the classical parabolic shape and the velocity peak at the pore wall disappears, but the velocity inflection points in the adsorption layer region and the central region of the pore can still be seen. Since different shale organic nanopores at different pore pressures have similar velocity distribution characteristics, the velocities at different pore pressures were normalized using type-II kerogen shale organic nanopore as an example (Figure 8(a)), and the normalized velocity values were obtained based on the average velocities at different pressures within the pore.Figure 8
Normalized velocity and CH4 density distribution in the nanopores of type-II kerogen under different pore pressures.
(a)(b)As can be seen from Figure8(a), the velocity distribution within the shale nanopores can be divided into three regions. The first region (Np-wall) is adjacent to the shale nanopore wall, where gas atoms are strongly repelled by the solid wall and almost no atoms or velocities appear; the second region (Np-adsorbed) is around the adsorption layer where complex phenomena exist, the slip velocity decreases with increasing pore pressure and disappears at high pressures, and the velocity values are reduced at the potential well where density peaks exist compared to the other regions; the third region (Np-bulk) is the central bulk phase region of the nanopore, where the transport behavior mainly follows the classical parabolic pattern, but with different curvatures due to the rarefaction effect in the nanopores[45].It can be concluded that the CH4 molecules in the adsorption layer exhibit completely different transport characteristics compared to the bulk phase region in the center of the nanopore. This is mainly due to the different density distributions of the adsorbed layer and the bulk phase region within the shale nanopores leading to their different velocity distributions. As an example to illustrate the results of the density distribution curve of type-II kerogen (Figure 8(b)), when the pressure was low (10 MPa), the low peak at the nanopore wall in the density profile indicates that the adsorption and nanospatial confinement of the nanopore wall were weak, and fewer CH4 molecules were adsorbed on the nanopore surface, so it corresponds to Figure 8(a), which the transport velocity is higher near the nanopore wall. As the pore pressure increases, the density of CH4 molecules increases within the nanopore, and the difference in density between the Np-adsorbed region and the Np-bulk region decreases. When the pressure reached 50 MPa, due to the strong confinement effect in the shale nanopores under high pressure conditions, the velocity profile showing a nonslip boundary feature. That can be elaborated as the pressure increases; the main kinetic behavior of CH4 molecules in shale nanopores changes from collisions between gas molecules and nanopore wall to gas intermolecular collisions, which leads to significant differences in gas transport characteristics at different pressures.Figures7(b), 7(f), and 7(j) show the results of the velocity distribution under the influence of different pressure drop. It can be seen that the velocity distribution within the shale nanopore is positively correlated with the pressure drop. The approximate parabolic shape of the simulated CH4 velocity profile obtained in the center of the nanopore again indicates the presence of viscous flow during the transport of supercritical state CH4 within the shale nanopore. Furthermore, by observing closely at the different velocity profile values, especially in the Np-bulk region, it can be found a near-linear relationship between pressure drop and the velocity values. When the velocity values about 0.17 Å/ps at 1.25 MPa pressure drop are 1.7 times larger than ~0.1 Å/ps at 0.75 MPa pressure drop and are 5 times larger than ~0.03 Å/ps at 0.25 MPa pressure drop.Moreover, the results of the velocity distribution of CH4 molecules under the influence of different pore sizes are analyzed in Figures 7(c), 7(g), and 7(k). It can be clearly seen that the velocity profile of CH4 molecules within the shale nanopores is highly dependent on the pore size, with higher values of flow velocity for larger pore sizes. And further comparative analysis will reveal that the velocity profile shows a nonlinear relationship with pore size, as exemplified by the type-II kerogen (Figure 7(g)), where the maximum velocity in the central region of the nanopore (~0.27 Å/ps) is about 2.7 times higher for the 4 nm pore size condition than for the 2 nm pore size case (~0.1 Å/ps) and 6.75 times (~0.04 Å/ps) for the 1 nm pore size case. This illustrates the large effect of pore size on the gas flow capacity within shale nanopores, where several-fold differences in pore size can lead to flow rate differences of almost more than an order of magnitude. This is consistent with continuous flow theory (Poiseuille’s law), where the gas flow rate within the channel is proportional to the square of the channel size and the underlying mechanisms were revealed by previous research results [46, 47].Finally, the simulation results of CH4 molecular velocity distribution within the shale nanopores at different ambient temperatures were compared. From the velocity distribution of CH4 molecules within shale nanopores in Figures 7(d), 7(h), and 7(l), for the effect of ambient temperature, it can be found that the difference between velocity distributions under different ambient temperature conditions is small, indicating that the increase in temperature mainly accelerates the thermal movement of molecules and has little enhancement for the velocity in the transport direction. Therefore, the transport capacity of the gas can be roughly considered to be the same under different shale reservoir temperature conditions.In summary, the above analysis of the transport simulation results illustrates that the transport pattern of CH4 in shale nanopores is complex and different from the behavior in homogeneous carbon nanotubes or graphene nanochannels. The nonhomogeneous nanopore walls result in velocity profiles that are not symmetrically distributed, and the chemical structures of different kerogen types alter the interaction between CH4 and shale nanopore surface, which affects the molecular transport velocity. Conventional understanding ignores the shale gas transport pattern with adsorption, which may underestimate the gas permeability in shale reservoirs. In a future research work, we will compare the MD results with the classical theoretical model values and then develop a new model for the apparent permeability of gas in shale organic nanopore, trying to applicate our simulation work to real shale matrices with multiscale complex pore networks.


## 3.3. Displacement of the Adsorbed CH4

To improve the recovery of adsorbed CH4, gas injection displacement is an efficient method. According to the research investigation, the commonly chosen injection gas in enhanced gas recovery engineering is CO2 and N2 [48–50]. Compared to CH4, CO2 has a higher adsorption capacity and shale matrix has strong surface heterogeneity, in which the oxygen-containing functional groups on the surface of shale nanopore enhance the displacement of CO2 to CH4. Not only that, CO2 can be sequestered in the shale formation by means of gas injection. In addition to CO2, N2 is also a common choice for gas injection. According to Dalton’s partial pressure law, with the injection of N2, the partial pressure of CH4 will decrease. Besides, N2 has many advantages, such as being abundant, readily available, and nontoxic and noncorrosive.This section focuses on the displacement of adsorbed CH4 by the injection gas, which is still implemented using a MD method and with the help of the LAMMPS open source program. The model schematic diagram is shown in Figure 9; the injection gas is either CO2 or N2. In the initial state, both boxes were not connected to the shale organic nanopores. The injection gas was equilibrated in the left box while the CH4 was adsorbed in the nanopore. In this section, the selected working condition is that the injection gas pressure is 12 MPa while the adsorbed CH4 pressure is 10 MPa. In addition, the temperature of the whole system was controlled using a Nosé-Hoover thermostat at 338 K. Subsequently, create the entrance between the boxes and nanopore, and dig a hole (remove atoms) in the wall between the nanopore and the boxes on both sides, which is the size of the cross section of the nanopore. Then, the injected gas entered the nanopore for the displacement under the pressure difference. In subsequent simulations, the content of each component of the output gas in the right-hand box was counted to analyze the displacement mechanism and process.Figure 9
Molecular model of CH4 displacement by injection gas. The gray wall at the left-hand side represents piston to maintain the injection pressure in box. The injection gas is either CO2 or N2 and in the schematic diagram is CO2.The process of CH4 displacement by CO2 and N2 was studied comparatively. Figure 10 gives a graph of the variation in the concentration of each component of the output gas after the displacement of CH4 by CO2 and N2. When CO2 was injected, the mole fraction change curve shows a trend of long breakthrough time and sharp front, as shown in Figure 10(a). Since CO2 tends to adsorb on the shale nanopore surface, the amount of CO2 in the bulk state is less and its passage through the nanopore was slower. When the injected gas comes in contact with the adsorbed CH4, a front was formed. When CO2 breaks through the nanopore space, the concentration of output CH4 rapidly decreases to zero. This also indicates a sharper CO2 front within the nanopore space. On the contrary, when N2 was injected, the mole fraction change curve shows a trend of short breakthrough time and a flat front, as shown in Figure 10(b). Since N2 has a weaker adsorption capacity, it can pass through the nanopore faster and accelerate CH4 desorption earlier. Even though N2 has already broken through the nanopore, CH4 continues to be produced. When enough N2 was injected, the concentration of output CH4 gradually tends to zero. These findings suggest that both CO2 and N2 can efficiently displace the adsorbed CH4. The above simulation results are consistent with the experimental results of related literature, which once again validates the accuracy of the simulation method.Figure 10
Displacement gas composition recorded from the MD as a function of displacement time. (a) CO2 displacing CH4; (b) N2 displacing CH4.
(a)(b)


## 4. Conclusion

Three shale organic nanopore models of different maturity were established based on real kerogen molecules. The GCMC, EMD, and NEMD simulation methods were used to study the microbehavior of CH4 under nanoscale conditions. First, the CH4 occurrence characteristics in shale organic nanopores were analyzed through obtaining the density distribution under different pressure conditions. Moreover, the transport characteristic was elaborated through the velocity distribution of CH4 molecules in shale organic nanopores under different conditions using control variates method. Lastly, the displacement of adsorbed CH4 by the injection gas was investigated by simulating the process of CH4 displacement by CO2 and N2 comparatively. The following specific conclusions were as follows:
(1)
CH4 is mainly adsorbed as a monolayer in shale organic nanopores; in the occurrence of gas within the nanopores, there are high-density adsorption layer region and bulk phase region; density ratios are negatively correlated with both pressure and pore size; the reduction of the pore size will help to improve the gas storage capacity in the nanopores(2)
The velocity distribution of CH4 within the shale organic nanopore can be divided into three regions; the flow velocity is positively correlated with pressure drop, pore size, and ambient temperature and negatively correlated with pore pressure; the rapid flow of CH4 is more favorable at lower pore pressures, at higher pressure drop or in larger pore size(3)
Both CO2 and N2 can effectively displace the adsorbed CH4; it is found that when CO2 was injected, the breakthrough time is long and the front is sharp and when N2 was injected, the breakthrough time is short and the front is flat



---

*Source: 2420869-2022-03-07.xml*