Adsorption and Diffusion Behaviors of CO₂ and CH₄ Mixtures in Different Types of Kerogens and Their Roles in Enhanced Energy Recovery

Abstract

CO₂ geological sequestration in subsurface shale formations is a promising strategy to store CO₂ and to increase shale gas production. The understanding of gas adsorption and diffusion mechanisms in microporous media is critical for CO₂ storage-enhanced gas recovery (CS-EGR). The type of kerogens is one of the important factors that influence the adsorption and diffusion behaviors of gases. In this work, the Grand Canonical Monte Carlo and Molecular Dynamics simulations were utilized to develop kerogen models and further investigate gas and water adsorption and diffusion behavior on the type IA, IIA, and IIIA kerogen models. The results indicated that the adsorption and diffusion capacities of CO₂ are larger than those of CH₄. The adsorption and diffusion capacity decreased with increasing water content. However, the CO₂/CH₄ adsorption selectivity increased with the increase in water content. Type IIIA demonstrated the best potential for adsorption and diffusion. This study provides insights into the role of the adsorption and diffusion behavior of CO₂ and CH₄ mixtures on kerogens of different types under different water contents at a microscopic scale, and can facilitate further understanding of the processes involved in CO₂ storage coupled with enhanced energy recovery.

1. Introduction

Emissions of greenhouse gases, especially CO₂, have been identified as the main contributor to global warming and climate change. Carbon capture and storage (CCS) is considered one of the most effective carbon reduction options, and can be achieved in several ways, primarily through geological sequestration, deep-sea storage, and mineral carbonation [1,2]. CO₂ geological sequestration is considered the most feasible method of storage, thanks to its unique economics, site availability, and environmental safety advantages. As a developing country, China has planned to reduce 40–45% of its CO₂ emissions per unit GDP by 2020 based on the 2005 level [3,4]. Most of the potential and ongoing sites for CO₂ geological sequestration in China include saline aquifers (Ordos, Daqing), oilfields (Daqing, Songyuan, Dongying, Puyang, Yanchang, Tianjin, Ordos), gas reservoirs (Sichuan, Chongqing), and coal fields (Jincheng) [5,6]. Carbon dioxide sequestration in oil and gas reservoirs is accepted as one of the most effective storage options. First, the infrastructure for oil and gas exploration can be directly used in the engineering of CO₂ storage after only minor modifications, which reduces the cost of construction [7,8]. Second, experience with CO₂-EOR can be applied to carbon dioxide storage projects, which further reduces the cost of technology research and development for CO₂ sequestration [9]. Third, CO₂ sequestration in oil and gas reservoirs enhances oil and gas production, reducing overall engineering costs [10,11]. Shale formations are potential sites for CO₂ geological sequestration [12], as the geologically porous structure of shale is suitable for gas storage. Meanwhile, the sequestrated CO₂ might enhance methane energy recovery due to the interaction between gases and increased pressure in shale formations.

According to the van Krevelen diagrams [13], the kerogen type is classified by the elemental ratios of hydrogen/carbon (H/C), oxygen/carbon (O/C), and sulfur/carbon (S/C). Type I kerogen is from the lacustrine anoxic environment, type II kerogen is from marine shale and continental planktons, type III kerogen is from plants in tertiary and quaternary coals, and type IV kerogen is from older sediments redeposited after erosion. The differences in kerogen are due to sedimentary organic matter maturation during the burial process [14]. Based on the International Union of Pure and Applied Chemistry, the different pores that exist in the reservoir can be classified into three basic categories: micropores (<2 nm), mesopores (2 to 50 nm), and macropores (>50 nm) [15,16]. Guo et al. used CO₂ adsorption, N₂ adsorption, nano-CT scanning, and micro-CT scanning to image the Longmaxi shale. Micropores (<2 nm), mesopores (2–50 nm), macropores (50 nm–2 µm), and microfracture (2–98.5 μm) accounted for 33.59%, 36.28%, 14.04%, and 16.09% of pore volume, respectively [17]. The mesopores were mainly hosted in the minerals, and the micropores were mainly hosted in the organic matter [18]. Moreover, organic matter is typically kerogen [19]. Kerogen is one of the major constituents (up to 30 vol%) of organic-rich rocks [20].

Shale gas can mainly be stored in three ways: free gas in fractures and macropores, adsorbed gas on kerogen and clay-particle surfaces, and dissolved gas in kerogen and bitumen. The contribution of dissolved is was about 22%. Almost all of the gas is compressed gas and adsorbed gas [21,22]. The CS-EGR technique is based on the competitive adsorption of gas and CO₂. In shale, the clay mineral is widely accepted as hydrophilic, while the organic matter is traditionally accepted as hydrophobic [23]. However, recent experimental and simulation work found that kerogen has both hydrophilic and hydrophobic features due to the carbon skeleton being hydrophobic, while the carboxylic and hydroxylic functional groups are hydrophilic [24,25]. In addition, the hydraulic fracturing process increases the water content of the reservoir by injecting a large quantity of water-based fluids. H₂O molecules can block the throat of the shale and occupy the gas molecule adsorption sites, thereby rapidly reducing the adsorption capacity [26]. Water content is a first-order controller of gas adsorption. The impact on gas adsorption is more significant than temperature or maturity [24,27]. To achieve large–scale field applications of CO₂ storage and improve the efficiency of the EGR process, the mechanism of CO₂ and CH₄ adsorption and diffusion should first be understood. However, this is difficult because of the chemically heterogeneous structure of natural kerogen and its complex gas adsorption and diffusion behavior, especially the complicated processes occurring in interfaces between the gas and kerogen with water content.

The adsorption and diffusion of gases have been studied experimentally in recent years. Experimental results demonstrate that CO₂ is favorably adsorbed over CH₄ when the two gas species coexist in shale [28]. The TOC content was the most important control factor for methane adsorption capacity, and hydrophobic and hydrophilic sites might influence the water distribution [29,30]. Studies have indicated that the differences in adsorption capacity among the different kerogen models could be attributed to the different chemical structures of the kerogens [31,32]. The methane adsorption capacity of kerogen as a function of burial depth was established based on the Langmuir model, indicating that the adsorption capacity was affected by temperature and pressure [33] Both the mineral and kerogen components of shale were studied by comparing shale and the corresponding isolated kerogens to evaluate the relative contributions of these components. Microscopic studies have revealed that the porous systems in clay-rich, organic-rich, and microfossil-rich parts of shale are very different, as well as the importance of the inter-granular organic-mineral interface [34]. Moreover, effective pore diffusion is sensitive to the atomic size of gas molecules [22]. Furthermore, surface diffusion gradually becomes a primary contributor to gas production with decreasing pore pressure and pore diameter as well as increasing surface diffusivity and adsorbed capacity of kerogen [35].

It is challenging to directly understand the adsorption and diffusion properties using experimental methods, especially in high-temperature and high-pressure conditions. Molecular simulation is a powerful tool for studying the physical and chemical properties at the molecular level. Molecular simulations, such as the density functional theory (DFT), grand canonical Monte Carlo (GCMC) method, and dynamics simulation (MD) have come into wide use for studying CO₂ and CH₄ adsorption and diffusion. The DFT is a simulation method for calculating a phenomenon from an atomic level [36,37]. The Grand Canonical Monte Carlo (GCMC) method is a simulation method for solving a phenomenon from a microscopic level, with the configuration generated in the simulation changes randomly until it approaches the statistical mechanical ensemble [38,39]. The molecular dynamics (MD) method tries to solve structural problems at the atomic length scale and obtain the dynamic properties of gas molecules, such as gas diffusion in kerogen [40,41]. Experimental and molecular simulations have investigated the adsorption behaviors of methane on organic and inorganic matter. The shale adsorption capacity decreased in the following order: organic pores > clay mineral pores > quartz pores [42]. The adsorption and diffusion behavior of pure CH₄ and pure CO₂ as a single component in dry kerogen type IIA were investigated as well. The adsorption capacity of pure CO₂ was higher than that of pure CH₄. The CH₄/CO₂ self-diffusion coefficient increased with the increase in temperature at the same pressure, and decreased with increasing pressure at the same temperature [43]. In addition to pure gas research, the adsorption of mixtures of CO₂ and CH₄ has been investigated. The results indicated that the type of kerogen was one of the important influencing factors, and that kerogen IIIA was the optimized organic type for CS-EGR due to its large and stable CO₂ storage capacity [44]. In natural shale formations, in addition to the main component, methane, the adsorption and diffusion behaviors of other gas components are important as well. The diffusivity of CH₄ and C₂H₆ in the kerogen matrix have been calculated at 298.15 and 398.15 K at various pressures; methane diffused ∼two times faster than C₂H₆, and the diffusivity increased with increasing temperature [45]. However, the detailed adsorption and diffusion mechanisms of CO₂ and CH₄ mixture on different kerogen models remain undetermined, and knowledge about the effect of kerogen type on CS-EGR application is limited.

In this work, the models for kerogens of type IA, type IIA, and type IIIA were constructed, and the adsorption and diffusion behaviors of CO₂ and CH₄ mixture on kerogens under different water contents at a microscopic scale were investigated using GCMC and MD methods. In addition, the pore structure of different kerogen models, adsorption capacity, isosteric heat of adsorption, CO₂/CH₄ selectivity, RDFs between gases and elements on kerogen models, and diffusion of gases of kerogen are discussed in this paper.

2. Calculation Method and Theory

2.1. Unit Model of Kerogen

Kerogen type is considered one of the main factors affecting the adsorption and diffusion behavior of gases. Type I and type II have low oxygen and high hydrogen content, are oil-prone generators, and can produce gas at high maturation [46,47]. Type III has a high oxygen and low hydrogen content, and is a gas-prone generator [46,47]. In comparison, type IV kerogen attracts less research and commercial attention because of its overmatured organic matter, and has few possibilities for oil and gas generation [48]. The unit models of kerogen types IA, IIA, and IIIA utilized in this work were taken from previous research [49,50]. Their structures and compositions showed good agreement with the experimental results of a combination of solid-state ¹³C NMR, Sulfur X-ray Absorption Near-Edge Structure (S-XANES), and X-ray Photoelectron Spectroscopy (XPS) techniques [49]. The chemical formulas of the kerogen models are roughly assumed to be C₂₅₁H₃₈₅O₁₃N₇S₃, C₂₅₂H₂₉₄O₂₄N₆S₃, and C₂₃₃H₂₀₄O₂₇N₄, respectively.

2.2. Modeling of Kerogen

The type IA, IIA, and IIIA kerogen models were developed using Materials Studio software. The COMPASS force field was utilized to build the kerogen models [51]. COMPASS is one of the all-atom force fields with high accuracy in calculating thermodynamic parameters of condensed models of organic and inorganic matters [52,53]. Moreover, CH₄, CO₂, and H₂O geometry parameters optimized by COMPASS are consistent with previous experimental research [54].

The geometry optimization function was used to relax the structures for the single unit of type IA, IIA, and IIIA. During the process, the Smart algorithm, which integrated the steepest descent method, ABNR method, and quasi-Newton method, was used to optimize the structures [55,56,57]. For both Coulombic and van der Waals energy summation, a fine cutoff distance of 15.5 Å was chosen for atom-based interactions. The Annealing function was used to obtain the stable structures with the NVT ensemble. The single unit of type IA, IIA, and IIIA was annealed with ten cycles in a temperature range from 300 to 800 K for 400 ps.

To obtain the bulk kerogen models for adsorption and diffusion simulation, ten stable unit kerogen models were formed into an amorphous cell with a density of 0.1 g/cm³ in the Amorphous Cell module as initial bulk kerogen models. First, to relax the kerogen models, the initial bulk kerogen models were annealed from 300 to 800 K with an NVT ensemble for 200 ps at 20 MPa. Then, to obtain the kerogen models similar to real kerogen, the relaxed kerogen models were annealed with NPT ensemble at 900, 700, 500, and 300 K for 300 ps at 20 MPa, respectively. Finally, to obtain the stable bulk kerogen models, the kerogen models were annealed at 300 K for 400 ps at 0.1 MPa. The Andersen temperature thermostat and Berendsen pressure barostat were chosen to generate the positions and velocities during the annealing performance.

For the moist kerogen models, a certain number of H₂O molecules were adsorbed on dry kerogen models as different water contents (0.6, 1.2, 1.8, and 2.4 wt.%) using the Fixed loading function in the Sorption module.

2.3. Adsorption Simulation Details

To investigate the adsorption behaviors of CO₂ and CH₄, the GCMC method was used in the Sorption module. The COMPASS was the force field. Furthermore, the Ewald and Atom-based methods were used for the Coulombic and van der Waals interactions, respectively [43,58]. In this work, the fugacity, which combined the Dalton law and the Peng–Robinson equation instead of pressure, was used during the GCMC simulation [59]. For each specific pressure, a total of 1 × 10⁷ calculation steps were performed to obtain data related to the adsorption process for further analysis. The first 5 × 10⁶ steps were used to relax the kerogen models to the equilibrium state, and the other 5 × 10⁶ steps were used to obtain the simulation data.

To compare our simulation results with previous experimental data, the absolute adsorption capacity was converted to excess adsorption capacity by the formula as follows:

$$n^e=n^a-v\rho$$

(1)

where $n^a$ is the absolute adsorption, $n^e$ is the excess adsorption, $v$ is the free volume, and $\rho$ is the density of particles calculated by the Peng−Robinson equation at a specific temperature and pressure [59].

Isosteric heat of adsorption is an important thermodynamic parameter for evaluating the interaction between CO₂ or CH₄ molecules and the kerogen model. It can be determined by the Clausius–Clapeyron equation [60]:

$$q_{st}=R{T}^2{\left(\frac{\partial \ln p}{\partial T}\right)}_n$$

(2)

where $q_{st}$ is the isosteric heat of absolute adsorption of $n$, kJ/mol; $R$ is the universal gas constant, kJ/(mol·K); $T$ is the temperature, K; and $p$ is the pressure, kPa.

Adsorption selectivity is an important parameter used to evaluate the ability of an adsorbent to preferentially adsorb substances in a binary mixture system. An adsorption selectivity value$S_{i/j}$ greater than one means that component i is preferentially adsorbed by the other components j in a binary mixture system. The adsorption selectivity of CO₂ over CH₄, $S_{C{O}_2/C{H}_4}$, is defined as follows:

$$S_{C{O}_2/C{H}_4}=\frac{x_{C{O}_2}/x_{C{H}_4}}{y_{C{O}_2}/y_{C{H}_4}}$$

(3)

where $x_i$ is the mole fraction of gas $i$ in the adsorbed phase and $y_i$ is the mole fraction of gas $i$ in the kerogen bulk phase.

To evaluate the accessibility of kerogen models we use the term “degree of connectivity”, defined as follows [61]:

$$Degree of connectivity=\frac{accessible pore volume}{total pore volume}\times 100\%$$

(4)

The adsorption properties, such as porosity, degree of connectivity, cumulative pore volume (CPV), and pore size distribution (PSD), were calculated by inserting probes using the Atom Volumes Surface Tool function. Probes with different radii from 0.5–3 Å were randomly inserted into the bulk kerogen models. The space where the probe space did not overlap the bulk kerogen model was regarded as free pore volume.

Radial Distribution Function (RDF) is a physical quantity used to investigate the distribution of CO₂, CH₄, and H₂O in bulk kerogen models. The expression for RDF is as follows:

$$g_{ab}\left(r\right)=\frac{dN}{4\pi {\rho }_b{r}^{2}dr}$$

(5)

where $dN$ is the number of b particles at a distance r to r + dr from a particles and ρ_b is the bulk density of b particles.

2.4. Diffusion Simulation Details

The self-diffusion coefficient is commonly measured by experimental techniques such as Quasi-Elastic Neutron Scattering (QENS) and Pulsed-Field Gradient-Nuclear Magnetic Resonance (PFG-NMR) [62]. In addition, the self-diffusion is determined by Mean Square Displacement (MSD) and the Einstein method. It can be calculated as follows:

$$D_s=\frac{1}{6N}\underset{t\to \infty }{\lim }\frac{d}{dt}{\displaystyle \sum }_{i=1}^N{\left[r_i\left(t\right)-r_i\left(0\right)\right]}^2$$

(6)

where $D_s$ is the self-diffusion coefficient, N is the number of diffusion molecules, t is the calculation time, $r_i\left(t\right)$ is the displacement vector at time t, and $r_i\left(0\right)$ is the displacement vector at the initial time. The angular brackets denote that the quantity is an ensemble average property.

Based on the adsorption results, Molecular Dynamics simulation was performed to obtain CO₂, CH₄, and H₂O diffusion coefficients in kerogen models. The COMPASS force field was utilized in the Dynamics module. The Atom-based and Ewald methods were used for the van der Waals and electrostatic interactions. To relax the system, the kerogen models were annealed with an NVT ensemble for 2000 ps at 338 K, and a Nose thermostat was used to generate the positions. Then, the kerogen models were annealed with an NPT ensemble for 2000 ps at 338 K to ensure that the data reached the linear region for further analysis. The Nose temperature thermostat and Berendsen pressure barostat were used to generate the positions and velocities.

This was work performed with long simulation durations and enough sampling steps to ensure sufficient data for further analysis. In addition, we constructed three different configurations for each dry or moist kerogen model to eliminate possible systematic errors. The adsorption and diffusion results were from the average of three simulation calculations.

3. Results and Discussion

3.1. The Rationality of Kerogen Models

The changes in density during model construction are shown in Figure 1a. The initial bulk kerogen models were annealed from 300 to 800 K with an NVT ensemble for 200 ps at 20 MPa. Then, the relaxed kerogen models were annealed with NPT ensemble at 900, 700, 500, and 300 K for 300 ps at 20 MPa, respectively. Finally, the kerogen models were annealed at 300 K for 400 ps at 0.1 MPa. The simulated density of kerogen type IA, type IIA, and type IIIA was 1.01 ± 0.01 g/cm³, 1.13 ± 0.01 g/cm³, and 1.20 ± 0.01 g/cm³, respectively, under specific conditions. The simulated results were reasonably close to the experimental densities of kerogen type IA, type IIA, and type IIIA at 1.04 g/cm³ [63], 1.20 g/cm³ [64], and 1.25 g/cm³ [65], respectively. The differences between simulated results and experimental results from the literature were possibly due to different initial configurations.

Three schematic illustrations of porosity are presented in Figure 1b. We utilized helium as a probe to determine the porosity, which is the same as in experimental studies [31]. The porosity of dry kerogen type IA, type IIA, and type IIIA were 9.19%, 11.74%, and 16.7%, respectively. The simulated results were consistent with the previous results, such as the porosity of type IA 9.8% [66], type IIA 12.5% [67], and type IIIA (4.45–22.50%) [68]. The difference between our results and previous results might be attributed to the test temperature and the initial constituents.

To validate the excess adsorption isotherms of CH₄, we selected previous research which had similar experimental conditions to our simulation conditions. The excess adsorption isotherms of CH₄ calculated from Equation (1) and compared with documented experimental results are shown in Figure 1c. The documented experimental data included Green River Shale [69], Woodford Shale [70], and Cameo coal [71]. Furthermore, the Vitrinite Reflectance (VRr) for these kerogen models ranges from 0.5–0.58%, indicating early mature kerogen types, and thus are representative of typical type IA [72], type IIA [73], and type IIIA [31], respectively. For comparison purposes, the experimental data were normalized by TOC to study the effect of organic matter. The CH₄ excess adsorption capacity of kerogen type IA at 14 MPa 338 K (0.22 ± 0.009 mmol/g) was close to the experimental results for Green River Shale at 14 MPa 338 K (0.489 mmol/g) [31]. The CH₄ excess adsorption capacity of kerogen type IIA at 13 MPa 338 K (0.34 ± 0.016 mmol/g) was close to the experimental results for Woodford Shale at 14 MPa 338 K (0.667 mmol/g) [31]. The CH₄ excess adsorption capacity of kerogen type IIIA at 13 MPa 338 K (0.95 ± 0.044 mmol/g) was close to the experimental results for Cameo coal at 13 MPa 338 K (1.12 mmol/g) [31]. The magnitude and tendency of the excess isotherms of CH₄ were consistent with experimental results. However, the excess isotherms in the simulation were not perfectly fitted with the data from the literature, which might be attributed to the different structures between the molecular models and experimental samples; inorganic matter in the experimental samples might be another impact factor.

According to the density, porosity, and excess adsorption isotherms, the different types of kerogen models in our work were reasonable for use in further investigation of the adsorption and diffusion behavior of CH₄ and CO₂.

3.2. Adsorption Behavior on Dry Kerogen Models

3.2.1. Pore Structure of Kerogen Models

The adsorption capacity showed a close relationship with the pore structure of kerogens. The porosity probed by helium is shown in Figure 2a. The dry type IIIA kerogen models were able to provide the most significant space for CO₂/CH₄ adsorption and diffusion, as the porosity of kerogen is type IA (9.19%) < type IIA (11.74%) < type IIIA (16.7%).

The accessible pore width was defined as the kinetic diameter above 4 Å [44]. The “degree of connectivity” was defined as the accessible pore width. The degree of connectivity, as plotted in Figure 2a, was 4.66%, 9.77%, and 21.16% for type IA, IIA, and IIIA kerogen, respectively. Type IIIA showed the highest degree of connectivity. This indicates that type IIIA kerogen is able to provide more connected channels for gas adsorption and diffusion.

Moreover, the cumulative pore volume (CPV) and pore size distribution (PSD) were calculated, as shown in Figure 2b. The CPV of dry kerogen type IA, IIA, and IIIA increased from 0.65 vol.% to 99.71 vol.%, 0.66 vol.% to 98.76 vol.%, and 0.54 vol.% to 91.63 vol.%, respectively, with the increase of the pore width from 0.2 to 6 Å. Although the cumulative pore volume showed a similar trend, dry type IA showed a higher volume than type IIA and type IIIA. The pore width from 0.2 to 6 Å covered more than 98% free volume for gas adsorption for type IA and type IIA. However, when the pore width was more than 6 Å, there was an 8.37% free volume of type IIIA. Figure 2b shows PSD (dV/dr) under different pore widths of dry kerogen of types IA, IIA, and IIIA. The PSD of the dry kerogen models shows a single-peaked type. The peak of the PSD for type IA is around 1.4 Å, indicating that the diameter of most pores in type IA kerogen was around 1.4 Å. The peak of the PSD for type IIA and IIIA is about 1.6 Å, indicating that most of the pores in type IIA and IIIA kerogen were around 1.6 Å [44]. Although the peaks are almost the same, type IIA has a higher peak than type IIIA, indicating that the pores of type IIA were more concentrated around 1.6 Å than in type IIIA.

3.2.2. RDFs on Dry Kerogen Models

The RDFs between CH₄/CO₂ ($y_{C{O}_2}=0.5$) and elements in kerogen of dry type IA, IIA, and IIIA models at 25 MPa and 338 K were calculated and plotted in Figure 3. Figure 3a,c,e shows that the close contact peaks between CH₄ and the C, H, and O elements of dry kerogen type IA, IIA, and IIIA, respectively, are not noticeable, implying weak interaction between CH₄ and carbon, hydrogen, and oxygen functional groups. The close contact peak between CH₄ and S elements of kerogen type IA and IIA are sharper than between CH₄ and other elements, indicating that the sulfur functional groups provided high energy positions for the adsorption of CH₄ in dry kerogen type IA and IIA. The close contact peak between CH₄ and N elements of kerogen IIIA is sharper between CH₄ and other elements, indicating that the nitrogen functional groups provided high energy positions for the adsorption of CH₄ in dry kerogen type IIIA.

Figure 3b,d,f shows that the close contact peaks between CO₂ and C and H elements of dry kerogen type IA, IIA, and IIIA are not noticeable, implying weak interaction between CO₂ and carbon and hydrogen functional groups. The close contact peak between CO₂ and O, N, and S of kerogen type IA and IIA is sharper between CO₂ and other elements, indicating that the oxygen, nitrogen, and sulfur functional groups provided high energy positions for adsorption of CO₂ in dry type IA and IIA kerogen models. The close contact peak between CO₂ and the O and N elements of kerogen IIIA is sharper between CO₂ and other elements, indicating that the oxygen and nitrogen functional groups provided high energy positions for the adsorption of CO₂ in dry kerogen type IIIA.

3.2.3. Adsorption Capacity on Different Dry Kerogen Models

The adsorption capacity of the CO₂ and CH₄ mixture ($y_{C{O}_2}=0.5$) on dry kerogen type IA, IIA, and IIIA at 338 K is shown as a 3D surface in Figure 4a and as 2D with error bands in Figure 4b–d. In general, the adsorption capacity increased with increasing pressure for all three types of kerogen models. The adsorption capacity of CO₂ was much larger than that of CH₄. Moreover, the adsorption capacity of CO₂ and CH₄ followed the sequence type IA < type IIA < type IIIA. The standard deviation error band showed acceptable uncertainty (less than 5%) for the CH₄ and CO₂ adsorption capacity. The uncertainty was more significant for CO₂ than for CH₄, and the uncertainty increased with increasing pressure.

The isosteric heats of CH₄ and CO₂ of dry type IA, type IIA, and type IIIA at 338 K were calculated, and are plotted in Figure 5a. The adsorption of CH₄ and CO₂ were both physical adsorptions as the isosteric heats were lower than 42 kJ/mol. This conclusion agrees with previous experimental results [44,72]. The isosteric heat of adsorption of CO₂ was larger than that of CH₄, indicating that CO₂ had a higher affinity to kerogen, which is consistent with previous results [44]. The order of the isosteric heats of adsorption was type IA < type IIA < type IIIA, which shows the same tendency as previous experimental results [31] and simulation results [44]. The isosteric heat might be influenced by two competing processes, namely, gas–kerogen and gas–gas interaction [44], and may be influenced by the O/C ratio as well [44]. For the gas–kerogen process, the isosteric heat decreased with increasing pressure. The gas molecules might prefer to adsorb to high-energy positions at a low pressure range, while they would tend to weak energy positions at a high pressure range. For the gas–gas process, the isosteric heat increased with increasing pressure, indicating more gas molecular collisions at a higher pressure range. The isosteric heat decreased from 0–5 MPa, which might be explained as the gas—kerogen interaction dominating the process. From 5 to 30 MPa, the isosteric heat increased slightly, implying that the gas–gas interaction dominated the process. Type IIIA had the highest isosteric heat, as it provided the largest space for gas–gas interaction. Moreover, the isosteric heat increased with the increase of the O/C ratio. The O/C ratio of kerogen was type IA (5.18%) < type IIA (9.52%) < type IIIA (11.59%).

The CO₂/CH₄ adsorption selectivity of dry type IA, IIA, and IIIA and the first derivative of selectivity are shown in Figure 5b. The CO₂/CH₄ adsorption selectivity is larger than 1, indicating that CO₂ had a larger adsorption capacity than CH₄ on kerogen. The adsorption selectivity was from 1.82 to 4.92, close to the simulation results on shale (1.87 to 6.97) [44] and coal (2.3 to 8.9) [58]. The general trend of selectivity of the type IA, IIA, and IIIA kerogen models increased with increasing pressure when pressure was less than 5 MPa, and decreased quickly to reach a constant. The first derivative of CO₂/CH₄ selectivity shows this tendency clearly. When the first derivative was more than zero the selectivity increased, as the pressure was less than 5 MPa. When the first derivative was less than zero, the selectivity decreased with increasing pressure. Moreover, the first derivative of selectivity reached a minimum of around 6 MPa, indicating that the gradient of selectivity became large at around 6 MPa. There was a nearly linear correlation between the polarizability of the gas molecule and the resulting enthalpy of adsorption of the adsorbate [44]. Although the polarizability of CO₂ (26.3 × 10⁻²⁵ cm³) and CH₄ (26.0 × 10⁻²⁵ cm³) was similar, CO₂ had a quadrupole moment of 13.4 ± 0.4 × 10⁻⁴⁰·C m² [74], while the quadrupole moment of CH₄ is a non-polar molecule [75,76]. The polar surface functional groups on kerogen increased the CO₂/CH₄ selectivity by more vital gas–kerogen interactions at the beginning. As the pressure increased, the competition between CH₄ and CO₂ molecules began to dominate the process, and the CO₂/CH₄ selectivity decreased and tended to be stable [77]. The CO₂/CH₄ selectivity sequence was type IIIA > type IIA > type IA. The difference in pore space for CO₂/CH₄ adsorption might be the main reason for this. Type IIIA had the most considerable degree of connectivity, and a lower degree of connectivity decreased CO₂/CH₄ selectivity.

3.3. Diffusion on Dry Kerogen Models

To characterize the migration of CO₂ and CH₄ molecules, the self-diffusion coefficients of CO₂ and CH₄ on different type kerogen models in the dry state at 338 K of 5 MPa, 15 MPa, and 25 MPa are shown in Figure 6.

The self-diffusion coefficient of CO₂ was always more significant than that of CH₄ in dry type IA, IIA, and IIIA kerogen models. The self-diffusion coefficient order of CH₄ and CO₂ molecules was D_s (CO₂) > D_s (CH₄), which was in the opposite order to their kinetic diameters of ${\sigma }_{C{H}_4}$ (3.8 Å) > ${\sigma }_{C{O}_2}$ (3.3 Å). The kinetic diameter is a sensitive measure when moving in highly restrictive environments [78]. The CO₂ molecules had a smaller diameter, and therefore diffused more easily in micropores, dominating the diffusion process.

Moreover, the diffusion of CO₂ and CH₄ showed a reduction trend with increasing pressure. The higher pressure indicates a higher concentrations of gas molecules, which led to strong interactions among gas molecules and slowed down gas diffusion.

The diffusion of CO₂ and CH₄ of different kerogen followed the order type IIIA > type IIA > type IA. Because the degree of connectivity followed the order type IIIA > type IIA > type IA, self-diffusion increased with increasing amounts of CO₂ and CH₄ molecules. This tendency is consistent with previous research [79].

3.4. Adsorption Behavior on Moist Kerogen Models

3.4.1. Pore Structure of Moist Kerogen Models

We investigated the effect of H₂O molecules on the pore structure of the different types of kerogen models by calculating the porosity and degree of connectivity. The results are plotted in Figure 7a. The order of the porosity and degree of connectivity was type IA < type IIA < type IIIA, and the tendency was not affected by H₂O molecules. With the increase in water content, the porosity and degree of connectivity decreased for each kerogen type. It is possible that the H₂O molecules might occupy the free volume for gas adsorption and diffusion, while the clustered H₂O molecules might block the channels for adsorption and diffusion. The dry type IIIA model had the maximum porosity, and decreased from 16.70% to 13.32% with increasing water contents from 0 to 2.4 wt.%. The dry type IIIA model had the maximum degree of connectivity as well, which decreased from 8.00% to 4.40% with increasing water content from 0 to 2.4 wt.%. By comparison, the type IA model had the minimum porosity, and it decreased from 9.19% to 6.98% with increasing water content from 0 to 2.4 wt.%. The degree of connectivity decreased from 1.70% to 0.81% with increasing water content from 0 to 2.4 wt.%. While type IA might have free volume for adsorption, the connectivity would not be effective due to the difficulty of forming connected channels for adsorption and diffusion.

The cumulative pore volume (CPV) and pore size distribution (PSD) of different types of kerogens with different water contents were simulated, and the results are shown in Figure 7b–d. With increasing water content, the CPV in the different types of kerogens increased. The CPV values were from 0 to ~99% when the pore widths were from 0 to 6 Å for type IA and type IIA. This implies that a pore width of 6 Å would cover all the free volume for adsorption and diffusion. For dry type IIIA, there was about 8% free volume of pore width above 6 Å, although this decreased to 1.5% when the water content increased to 2.4 wt.%. Therefore, clustered H₂O molecules might occupy the free volume with a pore width larger than 6 Å.

The PSD of different kerogen models with H₂O molecules was of a single-peaked type. With the increase of water content, the peak value of the surface area increased, as the H₂O molecules make the pore structure slightly more uniform. The peak value was about 1.4 Å for type IA and about 1.6 Å for type IIA and IIIA. The structures of the kerogen models were not changed significantly by H₂O molecules, as H₂O molecules affect the free volume of dry kerogens. H₂O molecules would not occupy all the free volume for gas adsorption.

3.4.2. RDFs on Moist Kerogen Models

The RDFs between CH₄/CO₂/H₂O and elements in moist kerogen models of type IA, IIA, and IIIA (water content 1.8 wt.%) at 25 MPa 338 K ($y_{C{O}_2}=0.5$) were calculated, and are shown in Figure 8. Sharper peaks between CH₄ and S elements of kerogen type IA and IIA were obtained in Figure 8a,d, which indicated that the sulfur functional groups provided higher energy positions for adsorption of CH₄ in moist kerogen type IA and IIA than those from carbon, hydrogen, and oxygen functional groups. In kerogen IIIA, CH₄ preferred to interact with nitrogen functional groups.

Figure 8b,e,h shows that the close contact peaks between CO₂ and C, H elements of moist kerogen type IA, IIA, and IIIA are not noticeable, which implies weak interaction between CO₂ and carbon and hydrogen functional groups. The close contact peaks between CO₂ and the N and S elements of moist kerogen type IA are sharper than between CO₂ and other elements, implying that the nitrogen and sulfur functional groups provided high energy positions for adsorption of CO₂ in moist kerogen type IA. The close contact peaks between CO₂ and the O, N, and S groups of moist kerogen type IIA are sharper than those between CO₂ and other elements, implying that the oxygen, nitrogen, and sulfur functional groups provided high energy positions for adsorption CO₂ in moist kerogen type IIA. The close contact peak between CO₂ and the O and N elements of kerogen IIIA are sharper than those between CO₂ and other elements, indicating that the oxygen and nitrogen functional groups provided high energy positions for the adsorption of CO₂ in moist kerogen type IIIA.

Figure 8c,f,i shows that the close contact peaks between H₂O and C, H elements of moist kerogen type IA, IIA, and IIIA are not noticeable, implying weak interaction between H₂O and carbon and hydrogen functional groups for kerogen type IA, IIA, and IIIA. The close contact peak between H₂O and the N and S elements of kerogen type IA are sharper than those between H₂O and other elements, indicating that the nitrogen and sulfur functional groups provided high energy positions for adsorption of H₂O in kerogen type IA. The close contact peak between H₂O and the O, N, and S of moist kerogen type IIA are sharper than those between H₂O and other elements, indicating that the oxygen, nitrogen, and sulfur functional groups provided high energy positions for adsorption of H₂O in kerogen type IIA. For type IIIA kerogen, the close contact peaks between H₂O molecules and the O and N elements on the kerogen are higher than those between H₂O and other elements, indicating that H₂O molecules might occupy the positions around the oxygen and nitrogen functional groups, thereby hindering the adsorption of CO₂ or CH₄ molecules.

Based on these results, oxygen, nitrogen, and sulfur functional groups were most likely to provide high-energy competitive adsorption sites between CO₂ and CH₄. However, H₂O occupied sites for gas adsorption, especially the site for CO₂ adsorption. This indicates that H₂O has a more significant effect on CO₂ than CH₄ from an elemental perspective.

3.4.3. Adsorption Capacity on Moist Kerogen Models

The adsorption capacities of the CO₂ and CH₄ mixture ($y_{C{O}_2 }=0.5$) at 25 MPa and 338 K are shown in Figure 9a. For both CO₂ and CH₄, the adsorption capacity decreases with increased water content. However, water has a more significant effect on the adsorption capacity of CO₂ than that of CH₄. Moreover, the reductions are more significant in dry states to 0.6 wt.% water contents. At low water content (<0.6 wt.% water content), the H₂O preferred to adsorb on the oxygen, sulfur, or nitrogen functional groups, which had the same tendency as CO₂ molecules. When the water content increased, H₂O preferred to aggregate into clusters in the pore space, leaving more space for CO₂ and CH₄ molecules. The sequence of the reductions affected by water contents was type IA < type IIA < type IIIA. The reduction effect was stronger on type IIIA, as type IIIA had more space for H₂O, CO₂, and CH₄ molecules to adsorb.

The isosteric heats of adsorption of CO₂ and CH₄ mixture ($y_{C{O}_2}=0.5$) on moist kerogen models at 25 MPa and 338 K are plotted in Figure 9b. Generally, the isosteric results indicate that CO₂ has a higher affinity to kerogens. The order of the isosteric heats of adsorption was type IA < type IIA < type IIIA, the same as when dry. The isosteric heat of CH₄ and CO₂ decreased when the water content increased to 1.2 wt.% and increased when the water content increased from 1.2 wt.% to 2.4 wt.%. This decrease might be attributed to the H₂O molecules tending to adsorb on the oxygen, sulfur, and nitrogen functional groups. The increase can be interpreted as the H₂O becoming aggregated into clusters and more oxygen, sulfur, and nitrogen becoming available again.

The CO₂/CH₄ selectivity with water is plotted in Figure 9b. The CO₂/CH₄ selectivity increases slightly with increasing water content from dry to 0.6 wt.%, then increases more quickly from 0.6 wt.% to 2.4 wt.%. At a low water content range, the H₂O and CO₂ both preferred to adsorb on the oxygen, sulfur, and nitrogen functional groups (type IA and type IIA) and oxygen and nitrogen functional groups (type IIIA). In contrast, CH₄ only preferred to adsorb on the sulfur functional groups (type IA, type IIA) and the nitrogen functional groups (type IIIA). Thus, the reduction in CO₂ was higher than that of CH₄. At a higher water content range, H₂O preferred to migrate and aggregate, thereby releasing more sites for CO₂ and CH₄ adsorption. As a result, CO₂ was reduced more quickly than CH₄, which is in agreement with experimental results [80]. The CO₂/CH₄ adsorption selectivity increased with increasing water content. Because the adsorption selectivity of the dry state was about 2–3, if selectivity is kept constant, CO₂ should decrease 2–3 times more than CH₄. However, the effect of H₂O was not strong enough to maintain this selectivity. The CO₂/CH₄ adsorption selectivity order was type IIIA > type IIA > type IA, indicating that type IIIA is the better option for the CS–EGR application than type IA and type IIA.

3.5. Diffusion Behavior on Moist Kerogen Models

To investigate the effect of H₂O on CO₂ and CH₄ diffusion behavior, the diffusion coefficient of CO₂, CH₄, and H₂O on kerogen models with different types and water contents at 25 MPa and 338 K are shown in Figure 10. Generally, the order of the self–coefficient of CO₂, CH₄, and H₂O followed type IA < type IIA < type IIIA, the same order as the degree of connectivity. With increasing water content, the diffusion coefficient of both CH₄ and CO₂ decreased, as H₂O molecules could occupy the pore space and hinder the diffusion channels of CO₂ and CH₄ [81,82]. When the water content was below 1.2 wt.%, the CO₂ and CH₄ dominated the pore space, leaving little room for the diffusion of H₂O. When the water content was above 1.8 wt.%, H₂O molecules dominated the pore space, and the diffusion coefficient of H₂O increased. The diffusion of H₂O molecules was relatively more significant than CO₂ diffusion, as H₂O molecules had ample space for movement. These results are consistent with those of the previous study ($D_{H_{2}O}$ (1.6 × 10⁻¹⁰ m² s⁻¹) > $D_{C{O}_2}$ (1.5 × 10⁻¹⁰ m² s⁻¹) > $D_{C{H}_4}$ (1.37 × 10⁻¹⁰ m² s⁻¹) at 298.15 K) [83].

4. Conclusions

In summary, we studied the competitive adsorption and diffusion of CH₄, CO₂, and H₂O on the different types of kerogens based on the GCMC and MD methods. The CO₂ adsorption capacity on the same kerogen model was more significant than CH₄. With the increase in the C/H ratio of the kerogen, both the CO₂ and CH₄ adsorption capacity increased. The CO₂/CH₄ selectivity increased at low pressure and then decreased to become stable with increasing pressure. The diffusion rate of CO₂ was more significant than that of CH₄, as the CO₂ molecules had a smaller diameter. The porosity, adsorption capacity, and diffusion rate decreased with increasing water content, as the H₂O molecules aggregated as clusters and blocked the pore channel. Although H₂O molecules had an adverse effect on adsorption and diffusion behaviors, the adsorption and diffusion capacity nonetheless followed the order IA < IIA < IIIA. While the water content reduced the adsorption capacity, the CO₂/CH₄ adsorption selectivity increased with increasing water content, which is beneficial for CS-EGR applications.