China has become the third-largest shale gas producer in the world, establishing commercial mining in Fuling, Changning–Weiyuan, Zhaotong, and other districts (Dong et al. 2016a, b; Dong et al. 2016a, b; Zou et al. 2015, 2016, 2018). However, the main shale gas producing areas in China are still confined to the Sichuan Basin and its periphery. The production capacity is limited, the differences between wells are significant, input and output have not yet become matched, and large-scale exploitation still faces major difficulties (Chen 2016; Chen et al. 2012; Guo 2016). Paleozoic shale gas reservoirs in south China are over-matured, characterized by deep burial, wide distribution, and a complex reservoir pore structure (Chen et al. 2010; Guo and Zhang 2014; Nie et al. 2009). These characteristics cause difficulties in understanding the pore properties of the reservoirs and the characteristics of shale gas occurrence, presenting a barrier to overcoming the exploration and development bottleneck (Chen et al. 2020).

Shale gas reservoirs exhibit significant heterogeneity in space–time after through multi-stage and multi-type geological transformation (Pandey et al. 2018). Their pore types and structures, too, show great variation on the micro and macro scales (Cao et al. 2015; Jarvie et al. 2007; Loucks et al. 2009; Ross and Bustin 2007, 2008). The pores in shale can be divided into organic pores and inorganic pores. A large amount of methane exists in the pores of organic matter in the form of adsorption. The process of maturation plays a critical role in the development of the organic pores in shale gas reservoirs (Timothy et al. 2020). Thermal evolution experiments with rich organic matter have shown that pore volume and surface area increase with thermal maturity. Micropores are added gradually with the increase of thermal evolution degree, and the proportion of micropores reaches its maximum synchronously with the peak of the hydrocarbon production phase of the thermal evolution process (Chen and Xiao 2014; Mathia et al. 2016). While micropores increase with maturity, the ratio of mesopores and macropores decreases when the temperature and pressure increase. Pore diameters change most obviously in the thermal evolution process between 1.5–7.5 nm and 60–70 nm. Micropores and mesopores are the main contributors to specific surface area (Chen et al. 2018). The combined results of experiments on the thermal evolution of shale and on methane adsorption under various evolutionary sequences have shown that pores are the most vital factor for methane adsorption performance in shale. With an increase in maturity, the TOC content gradually reduces, the microporous content gradually increases, and the increase in the microporous content during the thermal evolution process can significantly improve the capacity of methane adsorption (Singh 2011; Wang and Carr 2012; Singh et al. 2015; Pandey et al. 2018). However, the effect of the numbers of different types of pores appears to vary, and the augmentation effect of mesopores is smaller (Zhong et al. 2015). Thermal simulation experiments on oil shale samples suggest that the main diameter range first decreases and then increases with maturity. They also indicate that the generation of nanoscale cylindrical pores within organic matter is the primary cause of added porosity during mature hydrocarbon generation. Slit porosity at contacts between the organic matter and the skeleton particles, interlayer pores in clay minerals and dissolution pores also began to grow in the mature stage (Xue et al. 2015).

Adsorbed methane is an important part of shale gas, and the adsorption capacity of organic matter pores is stronger. It is of great significance to study the adsorption characteristics of methane in organic pores for improving shale gas productivity. Due to the complex and harsh geological conditions, it is difficult to meet the conditions in the process of experiment. More and more scholars use molecular simulation to study the adsorption characteristics of methane. Monte Carlo method and molecular dynamics method are common in the adsorption simulation system (Wu et al. 2015; Liu et al. 2016; Wang et al. 2016; Zhang et al. 2017). At present, molecular simulation methods have been used to study the adsorption characteristics of methane in organic pores, clay mineral pores and graphite pores (Xiong et al. 2017a, b, c, d), the competitive adsorption effect of CO₂ and CH₄ in kerogen (Sui and Yao. 2016), and the adsorption characteristics of methane on quartz with different wettability (Jiao et al. 2014). Molecular simulation has become an important supplement to experimental methods and an important tool to study adsorption characteristics.

Thermal evolution process leads to the development of multi-type, multi-shape and multi-scale organic pores in the reservoir. However, in the process of thermal evolution, the adsorption characteristics of methane in multi-type and multi-scale organic matter pores have not been fully understood. Based on the geological background and evolution of the Longmaxi Formation in the Changning area of the Sichuan Basin, computer-aided molecular simulations have been performed to investigate reservoir development of multi-type, multi-morphological pores and the molecular methane adsorption mechanism for multi-scale pores under a thermal evolution sequence. In this study, the following three problems are clarified: (1) the adsorption characteristics of methane in different types of pores with the same pore size; (2) the pore size of different types of pores with the strongest adsorption capacity; (3) the adsorption characteristics of methane in different evolution stages.

2.1 Background

Previous studies have analyzed the geological conditions, such as geological structure, organic geochemistry and pore structure parameters, of Longmaxi Formation shale in Sichuan Basin. (Chen et al. 2017) (Fig. 1). According to IUPAC, shale pore size is classified into micropores (< 2 nm), mesoporous pores (2–50 nm), and macropores (> 50 nm). Histories of thermal evolution and hydrocarbon generation indicate a division into the immature stage, with R_o (vitrinite reflectance) values of < 0.5%, a mature stage mainly producing condensate and moisture (0.5% < R_o < 1.2%), and an over-mature stage mainly producing dry gas (R_o > 0.5%) (Singh et al. 2013, 2016). The conditions used on the simulation were based these classification schemes as they relate to the Longmaxi Formation in the study area.

Geological background and condition settings for the simulation ((b) is modified from (Chen et al. 2017))

2.2 Simulation methods

2.2.1 Models

Due to prolong sedimentation, diagenesis, structure and other multi-type, multi-stage evolutionary effects of geological history, shale reservoirs develop extremely complex and heterogeneous pore systems (Fan and Liang 2019; Tang et al. 2019; Yang et al. 2016). It can be observed in the argon-ion polished scanning electron microscope image of shale rich in organic matter in the Longmaxi Formation, south Sichuan, that cylindrical pores, slit-shaped pores and slit-shaped pores with marginal grooves are ubiquitous in the reservoir. On the basis of the pore morphology found in the shale reservoir in this area, the pores are modeled as three simplified types: a homogeneous slit-pore model, a slit-pore model with grooves along its walls, and a cylindrical pore model. Cylindrical pores can be regarded as slit pores with curved sides. In order to make the physical properties of cylindrical pores and slit pores uniform, the physical properties (density) of homogeneous slit pores was taken as a reference, and the number of molecules on the upper and lower walls of a uniform slit pore were used as the basis for construction of a cylindrical pore (Fig. 2).

Pore models and parameters

2.2.2 Simulation details

The sorption module in Materials Studio (Accelrys) software was utilized to research methane adsorption behavior with the macro-regular Monte Carlo method (GCMC), and a Compass force field was adopted as the force field type (Mosher et al. 2013; Sun 1998; Liu et al. 2018, 2019; Hu et al. 2010; Zhang et al. 2019). The interactions in the simulation include van der Waals force interaction and Coulomb force interaction. For the specific expression of van der Waals force interaction and Coulomb force interaction, refer to Zhang et al. (2019) and Liu et al. (2019).

The temperature points of the adsorption simulation were set depending on the process of maturation and evolution of the Longmaxi Formation (Fig. 1a). The pore size sequences were divided into less than 2 nm (micropores), 2–5 nm (mesopores) and greater than 5 nm (macropores) in light of the IUPAC classification (Fig. 1c). There were six pore size points and seven temperature points (Fig. 1b, c). Pressure ranged from 0.1 to 30 MPa, and constant temperature and constant pressure were used to calculate on a point-to-point basis during the simulation process, which was divided into 12 pressure points for calculation. The thermodynamic characteristics and distribution of methane in pores were studied by molecular dynamics simulation using the stable configuration obtained after adsorption. The calculations for the molecular dynamics were made with the Focite module, adopting the same force field and intermolecular force settings as used previously. The canonical ensemble (NVT) was adopted, and an Andersen thermal bath was used to control the temperature. A 0.5-ns molecular dynamics calculation with a step length of 1 fs was carried out to determine the thermodynamic properties.

2.2.3 Excess adsorption amount

Absolute adsorption is defined as the quantity of gas present only in an adsorbed state (Mosher et al. 2013). When the temperature is lower than the critical temperature, the absolute adsorption amount represents the real adsorption amount. In shale formations, the methane adsorption becomes excessive when the temperature exceeds the critical temperature of methane (The critical temperature of methane is 191 K and the critical pressure is 4.6 MPa). In other words, absolute adsorption amount includes both gas molecules adsorbed on the pore wall in the form of adsorption phase and the gas molecules existing in the form of gas phase in the adsorption zone (near the pore wall) when the temperature exceeds the critical temperature. For excessive adsorption, the absolute adsorption amount is not the real adsorption amount. Therefore, under supercritical conditions, the excess adsorption capacity should be used to characterize the adsorption amount rather than the absolute adsorption capacity. Mosher and Gibbs have provided precise definitions of excessive adsorption (Mosher et al. 2013; Gibbs 1928). Excessive adsorption is defined as the additional amount of gas adsorbed per unit pore volume compared with the amount of gas in the same volume of a given pore without pore walls (Mosher et al. 2013). The excess adsorption amount only includes gas molecules adsorbed on the wall in the form of adsorption phase. The excess adsorption amount is equal to the absolute adsorption amount when the experimental temperature is lower than the critical temperature of the gas, whereas the excess adsorption amount is not equal to the absolute adsorption amount when the experimental temperature is higher than the critical temperature of the gas. Also, the difference between excess adsorption and absolute adsorption amount increases with the increasing experimental pressure when the temperature is higher than the critical temperature of the gas (Zhang et al. 2019). In this paper, the excess adsorption amount is used to express the adsorption amount. It can be found that it is reasonable to use excess adsorption amount to represent the real adsorption amount, which is in good agreement with the experimental results (Fig. 3). The calculation formula of excess adsorption amount is as follows:

Excess adsorption isotherm from simulation and experiments (the experimental data were reported by Bi et al. (2014, 2016))

where, N is the total amount of gas, which is fixed for the fixed model, mmol/g; V_p is the adsorbed phase volume which was occupied by some gas phase molecules, cm³; V_c is the gas phase volume, cm³; V_g is the free volume, approximately equal to the free volume of the pore, which is fixed for the fixed model and can be calculated by Connolly Surface in the simulation software, cm³; ρ_g is the density of the free gas, which is calculated by the P–R equation of state and is only related to temperature, g/cm³. It is worth mentioning that ρ_g will increases with the increase of pressure and may exceed N. Therefore, it is possible that the excess adsorption capacity is negative (Figs. 8b and 9e). The fundamental reason for this phenomenon is that the volume of adsorbed phase occupied by gas molecules varies with temperature and pressure, and we use fixed free volume to represent the volume of gas phase and the volume of adsorbed phase occupied by gas molecules. The reason for this is that we cannot determine the volume of adsorbed phase occupied by gas molecules, which is a difficult problem in adsorption science.

2.2.4 Probe molecule

Free volume is an essential parameter for determining the excess adsorption capacity in GCMC adsorption simulation. In GCMC simulation, confirmation of the free volume is usually achieved by introducing diverse molecules one by one into the pore interior as probes. The choice of probe molecules has a great influence on the accuracy of the simulation results. Since helium gas was utilized as the probe molecule in methane isothermal adsorption experiment used to obtain the free volume of experimental samples (Ji et al. 2012; Zhang et al. 2012), the free volume of pores was obtained by selecting helium as the probe molecule in this simulation to maintain consistency between simulation and experiment.

3.1 Analysis of simulation accuracy

The results for isothermal methane adsorption at T₁-type (slit) pores in the simulation and experiments were compared (Fig. 3). The simulated excess adsorption capacity matched well with the experimental results, confirming the reliability of the pore model used in the simulation.

The molecular radius of the helium probe used in this study is 0.13 nm, which is smaller than that of methane (0.19 nm). A helium probe can detect a finer pore space than can methane. Therefore, the accessible volume measured by helium probe is greater than the volume accessible for methane, and the calculated free volume is correspondingly larger than that of methane (Fig. 4).

Pores seen by different probes

Helium gas and methane were therefore respectively used as probe molecules to acquire diverse data on free pore volume, and corresponding adsorption curves for the two probes were compared (Fig. 5). The amount of adsorption in T_2-1 (grooved slit) pores calculated with a helium probe is less than with a methane probe. The amount of adsorption calculated for T₁-type pores using the two kinds of probes is almost the same; this is because T₁-type pores have a homogeneous structure, and there are no tiny spaces within them that methane molecules are not able to enter, unlike in the interior of T_2-1-type pores.

Excess adsorption amount with different probe molecules

3.2 Methane adsorption characteristics of different pore types

Figure 6 illustrates the comparative methane adsorption behaviors of the five morphological types of pores with the same pore size conditions. The molecular distribution of methane has two states. One is where methane molecules are gathered and accumulated near the wall of the pore area by the potential energy function of pore wall adsorption. Methane in this area includes adsorption-phase methane molecules and some free-phase methane molecules. The free-phase methane in the region is the cause of the decline in the adsorption curve. The molecular distribution pattern of methane is similar to that of the pore wall surface. Pore-groove methane has a concave distribution in the groove space, similar to the shape of the groove space (Fig. 6a). Within cylindrical pores, methane near the pore surface has a circular distribution along the pore wall (Fig. 6c). The other molecular distribution state is where methane molecules are less affected by the pore wall potential energy. This occurs where inter-molecular methane interactions occur far away from the pore wall. These molecules are dispersed, in free phase and have random distribution pattern.

Methane density and snapshot of different pore types

During adsorption, methane density is always higher at the position where the potential energy is the lowest (Liu et al. 2018; Charoensuppanimit et al. 2015; Hasanzadeh et al. 2010). There are peak densities at symmetrical positions at both ends and inside the slit-shaped pore methane density distribution curve because the potential energy of the pore wall is strongest near the pore wall region. The density peak reflects the methane adsorption layer near the pore wall with the most methane adsorption (Fig. 6a). In addition, we can also see that T_2-1, T_2-2 and T_2-3-type pores also have several significant subsidiary peak positions inside the pore wall. T_2-1-type pores exhibit only one subsidiary peak, as methane molecules cannot enter the second-layer groove space inside these pores because a pore space needs to be greater than 0.492 nm for methane molecules to be adsorbed. T_2-2 and T_2-3-type pores each have two subsidiary peak positions, corresponding to their two-layer groove space; the larger the amount of methane adsorption in the groove space, the higher the methane density peak. The methane density distribution curve for cylindrical pores only has a single density peak at symmetrical positions at both ends, reflecting the adsorbed methane layer along the pore wall. This reaches a significantly smaller density value than do the curves for slit-shaped pores (Fig. 6b).

Figure 7a shows isotherms of the methane adsorption capacity in the five types of pores with three different morphologies. Under different pressure conditions, the methane adsorption capacity of the slit pore with groove space is larger than that of the slit pore and cylindrical pore without groove space. It indicates that the presence of groove space has a significant additional effect on methane adsorption. When the growth rate of the four slit-shaped pore models is compared with the cylindrical pores under equal pressure conditions, it can be seen that the existence of groove space greatly influenced the pore methane adsorption capacity. The incremental growth rate compared with that in cylindrical pores reached more than 120% when the groove space opening size was 1.722 nm, showing a major improvement in methane adsorption capacity (Fig. 7b). Therefore, it can be considered that the micro grooves on the pore wall can greatly improve the methane adsorption capacity of the pores. The adsorption capacity of methane increases with the increase of groove space.

Methane adsorption characteristics of different pore types

For temperatures greater than 398 K (corresponding to R_o > 1.2%), the adsorption curve at pressures greater than 10 MPa was used to carry out further fitting analysis (Fig. 7c, d). From this, it is found that when groove space is formed and enlarges inside a pore, the slope of the curve fitting equation and the absolute value of the quadratic coefficient follow an auxetic law, showing that the presence of groove space enhances the pressure-sensitivity of the methane adsorption capacity of the pore. Multi-stage evolution over the course of the geological history of a reservoir results in a complex, irregular pore morphology. Such pores are more sensitive to pressure, and more irregular pores are more affected by pressure.

3.3 Methane adsorption characteristics at different pore scales

Figure 8a and b show isotherms of excess methane adsorption in different pore size bands. For slit-shaped pores less than 2 nm in diameter, when the pore size is less than 2 nm, when the pore size increases, the methane adsorption capacity increases; for those larger than 2 nm, methane adsorption capacity decreases as the pore size increases. Cylindrical pores, conversely, exhibit no obvious regularity.

Variation in the amount of adsorption with pore size

However, an insufficient division of pore diameter can obscure the variation law of the amount of pore adsorption in some pore diameter intervals. Therefore, more accurate relationship between adsorption capacity and pore size can be obtained by further subdividing the pore size range (Fig. 8c, d). As we can see, when using the more subdivided pore size ranges, the variation law of adsorption capacity observed is different from that concluded above. The variation in the adsorption capacity of slit-shaped pores with pore size shows a ternary-form distribution: within the aperture range of 0.7–0.9 nm, when the pore size increased, the adsorption changed little and increased slightly; it then shot up when the pore size is 1 nm before suddenly slowing down its increasing trend after passing 1.3 nm; in the pore size range above 1.5 nm, the adsorption capacity followed a law of stable decrease with the increase of pore size. The variation law of cylindrical pores with pore diameter is different from that of slit-shaped pores, though it also shows a ternary-form distribution: within the diameter range of 1–1.9 nm, the amount of adsorption shows a certain volatility with an increase in pore diameter; after 1.9 nm, the adsorption decreases rapidly and linearly with an increase in pore size. However, it can be implied from the combined results for slit-shaped and cylindrical pores that the pore size range most favorable for methane adsorption capacity is around 2 nm.

3.4 Pore adsorption characteristics under the thermal evolution sequence

The adsorption simulation results for the various thermal evolution temperature point settings are shown in Fig. 9. The results show that the forms of the methane adsorption curves for the three types of pores are roughly the same, all showed a tendency of increasing first and then decreasing with an increase in pressure. For temperature points corresponding to the thermal evolution sequence, the amounts of adsorption of the three types of pores all decreased with an increase in temperature. This occurs because methane adsorption is an exothermic process, and an increase in temperature will aggravate the molecular movement of methane, thus hindering methane adsorption.

Amount of methane adsorption under conditions corresponding with a thermal evolution sequence

The methane adsorption capacity of the grooved slit-shaped pores is larger than that of homogeneous slit-shaped pores, and the adsorption capacity increases gradually with an increase in groove space. Earlier studies have shown that methane cannot enter space with an opening size less than or equal to 0.492 nm and that pore space openings contributing to the adsorption capacity need to be larger than 0.492 nm (Zhang et al. 2019). For the slit-shaped pores with groove space, there is a difference in curve morphology between temperatures less than 398 K (corresponding to R_o < 1.2%) and greater than 398 K (corresponding to R_o > 1.2%). When the temperature is lower than 398 K, the maximum adsorption capacity is reached only at 5 MPa, and the process of reaching the maximum adsorption capacity is faster than that when the temperature is higher than 398 K. After reaching the maximum adsorption capacity, the adsorption capacity began to decrease, and the decrease rate decreased with the increase of temperature. When the temperature was greater than 398 K, the adsorption curve became relatively gentle, and the maximum adsorption capacity occurred at a pressure of 12 MPa or 16 MPa. The adsorption capacity then began to decline very gradually (Fig. 9a–d). This indicates that the pore is strongly adsorbant in the low-temperature stage of thermal evolution, but the amount of gas that can be adsorbed in this stage is greatly affected by pressure and the capacity for gas storage is weak, that is, a change in lithostatic stress will have a great influence on its gas-bearing capacity. In the high-temperature stage of thermal evolution, pore adsorption performance decreases, but the adsorbed gas content is less affected by pressure in this stage, and the gas storage performance is stronger and more stable (Zhang et al. 2019).

The amount of methane adsorption in cylindrical pores is smaller than in slit-shaped pores, indicating that cylindrical pores have weaker adsorption performance under the same temperature and pressure state. Unlike with slit-shaped pores, When the pressure was 4 MPa, the adsorption amount reached a maximum in cylindrical pores, then began to decline, and the rate of decline decreased with an increase in temperature (Fig. 9e). At the same temperature, adsorption in cylindrical pores decreased faster than in slit-shaped pores, confirming that cylindrical pores have weaker adsorption ability than slit-shaped pores. Similarly to slit-shaped pores, the curve morphology differs between temperatures less than and more than 398 K (corresponding to R_o < 1.2% and R_o > 1.2%, respectively). When the temperature is greater than 398 K, the pressure of maximum adsorption capacity increases. There is an intersection point in the cylindrical pore methane adsorption curve. The amount of adsorption is negatively correlated with temperature at the pressure less than at the intersection, whereas it is positively correlated with temperature at the pressure greater than that at the intersection. It is also found that the intersection point shifts to the right with increasing temperature when the temperature is greater than 398 K (corresponding to R_o > 1.2%), again showing that pore methane adsorption performance has higher stress sensitivity under low-temperature conditions, consistent with the conclusion in Sect. 3.2.

This study shows that the pore morphology and the change of temperature and pressure in the process of thermal evolution have an influence on the adsorption characteristics of organic pores. The groove space inside the organic pores can greatly improve the methane adsorption capacity. The pore with pore size of about 2 nm has the strongest adsorption capacity for methane. The pore adsorption capacity first increases and then decreases with the increase of pressure, and increases with the increase of temperature. In the early stage of thermal evolution, pore adsorption capacity is strong and pressure sensitivity is weak, while in the late stage of thermal evolution, pore adsorption capacity is on the contrary.

This will help researchers understand the adsorption characteristics of adsorbed gas in the process of thermal evolution, and lay a foundation for understanding the dynamic changes of adsorbed gas in the process of thermal evolution and the dynamic transformation of adsorbed gas and free gas. Based on this study, we believe that in the study of molecular diffusion and percolation mechanism of shale gas based on various adsorption desorption models, the morphological parameters of pore structure should be considered to make the model more accurate. At present, the understanding of pore structure evolution in the process of thermal evolution is not clear. In the process of thermal evolution, the change of pore morphology, temperature and pressure will affect methane adsorption. Therefore, the further research should pay attention to the coupling of these three factors and other factors on the dynamic control effects of adsorption, so as to get the research results that fit the real geological conditions as much as possible.