The simulation system is composed of adsorbate including methane, carbon dioxide and adsorbent, such as montmorillonite, illite, as well as kaolinite. The montmorillonite model adopted in our simulation is sodium-saturated Wyoming-type montmorillonite with the unit cell chemical formula of Na_0.75 [Si_7.75Al_0.25] (Al_3.5Mg_0.5)O₂₀(OH)₄, comprising typical tetrahedral-octahedral-tetrahedral (TOT) layers (Skipper et al., 1991). On the basis of the formula, one Si⁴⁺ is substituted by Al³⁺ every 32 Si⁴⁺ in the tetrahedral sheet, while one Al³⁺ is substituted by Mg²⁺ every 8 Al³⁺ in the octahedral sheet. And the interlayer Na⁺ can be used to equilibrium negative charge caused by these isomorphic substitution. Illite is also a typical 2:1 clay mineral, and the structure parameters of illite unit used in this work are established by Drits et al. (2010), with the chemical formula of KAl₄(Si₇Al)O₂₀(OH)₄. Isomorphic substitutions are achieved by substituting Si⁴⁺ by Al³⁺ every 8 Si⁴⁺, and the interlayer K⁺ is used to balance the negative layer charge. Kaolinite is composed of alumina octahedral and silica tetrahedron, showing typical characteristics of 1:1 clay mineral, with the composition of Si₄Al₄O₁₀(OH)₈. The structure parameters of the kaolinite unit are determined by Bish and Von Dreele (1989). The simulation box contains two layers of the clay sheets with different basal spacings (1nm, 2nm and 4 nm) in between, forming slit-like pore structure.

The ClayFF force field, based on precise representation of the metal-oxygen interactions between hydrated crystalline compounds and aqueous solutions (Cygan et al., 2004), is selected for the clay minerals. The potential model taken for the CO₂ molecule is from the study of Cygan et al. (2012). This potential model can not only better reproduce the physicochemical properties of CO₂, but also accurately describe the interaction between CO₂ and silicate minerals. The charge of carbon and oxygen atom in CO₂ molecule is +0.6512e and −0.3256e, respectively. Furthermore, the potential model adopted for CH₄ molecule is obtained from the TraPPE force field, which is created to describe thermodynamic characteristics of alkane (Martin and Siepmann, 1998). The hydrogen and carbon atom in CH₄ molecule are treated as a united atom without charge.

The adsorption characteristics of pure CO₂ and CH₄ as well as their mixtures in slit-like pores of clay minerals are determined using the GCMC method. Generally, the GCMC simulation is conducted in the μ V T ensemble, in which the system temperature T, system volume V and chemical potential μ and are fixed as constant. Note that the chemical potential μ is related to the fugacity, and the value of fugacity is achieved by Peng-Robinson equation of state (Peng and Robinson, 1976). In this work, the simulation temperature is set to 313 K and pressures up to 23 MPa. The van der Waals force interaction in the simulation is calculated by atom based method, and the Coulombic interaction is obtained by Ewald and Group method. During each simulation, a total of 2×10⁷ cycles are conducted, wherein the first 1×10⁷ cycles are performed to guarantee equilibrium and the remaining cycles are for statistical adsorption amount. According to previous research of Rani et al. (2019) and Aringhieri (2004), the development of micropores and small mesopores in the internal crystal layer of clay minerals is one of the major contributors to the surface area of shale and plays critical role in the occurrence of gases. Thus, the basal spacings of 1nm, 2nm and 4 nm are adopted to discuss the effect of pore size on the adsorption behaviors of gases. In this work, the cations and the clay sheets in the models are treated as rigid to simplify the simulation. Furthermore, periodic boundary conditions are applied to mimic crystalline periodicity.

In general, the adsorption amount achieved directly by experimental measurement is the excess adsorption amount, while it is worthy noting that the output of GCMC adsorption denotes the total loading number of molecules known as total adsorption amount. According to relevant literature (Rouquerol et al., 1999), the excess adsorption amount is defined as the product of adsorbed volume with the difference of adsorbed density and bulk density. Therefore, during the simulation, the excess adsorption amount per surface area of the adsorbent could be achieved by the Eq 1 (Pang et al., 2019). n e x = n t N A − ρ b × V f M a d / S a (1) where n e x represents the excess adsorption amount, N A means the Avogadro constant with the value of 6.02 × 10²³mol^-1, n t denotes the total loading number of the adsorbate, and the M a d refers to the mole mass of adsorbate. The bulk density ρ b at different simulation conditions is from the National Institute of Standards and Technology (NIST) Chemistry WebBook, while the free volume V f and surface area S a of the slitlike nanopore in the clay mineral models are achieved by the “Connolly surface” method using corresponding gas molecules. In this work, the excess adsorption isotherms obtained from the simulation are described by the Dubinin-Radushkevich-based model, which has been widely applied to study the CH₄ and CO₂ adsorption behaviors in coal and shales as shown in Eq 2(Dubinin, 1960; Sakurovs et al., 2007; Ozdemir, 2016): n e x = n 0 1 − ρ b ρ a exp − D ln ρ a ρ b m + k ρ b 1 − ρ b ρ a (2) where n 0 denotes the adsorption capacity of gases, ρ b and ρ a means the bulk density and adsorbed density, respectively. D represents a constant value associated with the affinity of adsorbate, k is the interaction coefficient between adsorbate and adsorbent, m is a constant, generally taking an integer in the range of 1–6. When m = 2, it is the modified supercritical DRk model established by Sakurovs et al. (2007).

The reasonability of these clay mineral models as well as the accuracy of simulation method need to be verified by comparing the simulation amount with experimental results. In general, the excess adsorption isotherms obtained by GCMC method are normalized by surface area of the basal spacing to unify the same comparison standard with the experimental measurement (Chen et al., 2017). Figure 1 shows excess adsorption isotherms of pure CH₄ and CO₂ on illite calculated in this work and relevant experimental data documented by previous literatures (Heller and Zoback, 2014; Jeon et al., 2014). It is observed that these excess adsorption isotherms obtained from experimental test and simulation are in the same order of magnitude and exhibit good consistency. Meanwhile, it should be noted that the slight deviation between the simulated and experimental isotherms is a normal phenomenon, because the clay crystals used in the simulation are ideal models, which may deviate from the minerals in the actual experiment. Besides, the pore size used for each adsorption simulation is a fixed value, while the pores in natural clay minerals are featured by multi-scale characteristics. In addition, the values of surface area used in the experimental tests were obtained from low-pressure N₂ adsorption method, which cannot reflect the characteristics of micropores. Overall, the simulation results are acceptable, and these clay mineral models and GCMC method can be applied to further investigate the interactions between gases and clay minerals.

Figure 2 depicts the adsorption isotherms as well as the fitting curves of CH₄ using DR-based model in clay minerals with varying pore sizes (1nm, 2nm and 4 nm) at 313 K. It is observed that the shape of the adsorption curves can be classified as Type-I according to the classification of Aranovich and Donohue (1998), reflecting the characteristic of microporous adsorbent. The total adsorption isotherms of CH₄ in different clay minerals display a similar variation tendency, that is, the total adsorption amount of CH₄ gradually increases with the increase of pressure and pore size. This is because limited storage space of smaller pores is not conducive to the occurrence of gas, and a larger pore diameter can provide more space for the loading of molecules. The excess isotherms demonstrate a maximum with increasing pressure, showing typical characteristic of high-pressure adsorption (Zhou et al., 2018). However, it should be noted that the increase of total adsorption amount is not entirely from the contribution of adsorbed molecules. According to the fitting results of DR-based model as shown in Table 1, the increase of pore size will lead to the continuous decrease of CH₄ adsorption capacity and the adsorbed phase density in various clay minerals, which may be due to the decrease of the coupling surface-gas interaction from the two sides of the walls as basal spacing increases from 1 nm to 4 nm.

On the other hand, it is found that the maximum CH₄ adsorption capacity of kaolinite is larger than that of illite and montmorillonite under the same pore size. This is an interesting phenomenon. As a non-polar molecule, the interactions between clay minerals and CH₄ are dominated by van der Waals force (Wang and Huang, 2019). Different from montmorillonite and illite crytals, the surface of kaolinite model used in this work contains a large number of hydroxyl groups, which may strengthen the interaction mechanism of CH₄ with kaolinite. Similar observations have been noticed by previous literature (Liu and Hou, 2020). In order to gain insights into these different adsorption behaviors of CH₄ in clay minerals, the isosteric adsorption heat of CH₄ in montmorillonite, illite and kaolinite under 20 MPa is quantitatively evaluated by Eq 3(Fokion and (3Alan) L, 1991): Q s t = R T + μ i n t r a N − E (3) where the Q_st denotes the isosteric adsorption heat between adsorbent and adsorbate, T is the temperature in the system, R is the universal gas constant, μ i n t r a represents intramolecular chemical potential, N and E denote the ensemble averaged molecular number of sorbate and total energy, respectively. Figure 3 illustrates that the value of isosteric adsorption heat between kaolinite and CH₄ is larger than that between montmorillonite, illite and CH₄ under the same basal spacing, demonstrating that the interaction between kaolinite and CH₄ is stronger. Furthermore, the adsorption heat gradually decreases with the increase of pore size, resulting from the decrease in overlapping effect of the two walls, which is consistent with the decreasing trend of the maximum adsorption capacity of CH₄.

A comparison of CO₂ isotherms in three clay minerals with different basal spacings is illustrated in Figure 4 as a function of pressure. Interestingly, similar Type-I shape presented in CH₄ isotherm is also observed in CO₂ adsorption curves. The variation in CO₂ adsorption amount of the different clay crystals exhibits a similar pattern. It is clearly seen that the loading number of CO₂ molecules obviously increase with increasing pressure. For the nanopores with 1 nm or 2 nm basal spacing, the adsorption amount of CO₂ is more easily to approach saturation as the adsorption pressure increases from 1 MPa to 9 MPa. In addition, it appears that the pore size of the clay crystals plays a more complex role in CO₂ adsorption behavior. A careful analysis of Figure 4 further demonstrates that in low pressure range the total adsorption amount of CO₂ in micropore (1 nm) is larger than that in mesopore (2nm and 4 nm), which may result from the overlapping effect of the stronger interactions between CO₂ and the cation exchange near the two clay planes. This phenomenon is not unusual and similar observation has been reported by Yang et al. (2015). Moreover, it is also seen from Figure 4 that the CO₂ excess adsorption amount first reaches a maximum at around 10–15 MPa, and then decreases with the increase of pressure, even becomes negative at higher pressures. The nonmonotonic behavior of CO₂ excess isotherm with maximum has been widely observed on various adsorbents by both simulation and experiment (Pini et al., 2010; Liu and Wilcox, 2012; De Silva and Ranjith, 2014; Merey and Sinayuc, 2016).

Table 2 gives the maximum CO₂ adsorption capacity of the three clay minerals with different basal spacings. It appears that the interaction mechanisms between CO₂ and clay minerals are more complicated compared with CH₄. As can be seen from Figure 4; Table 2 that the sequence in maximum CO₂ adsorption capacity obtained by the DR-based model is in the order of montmorillonite > illite > kaolinite. Interestingly, kaolinite shows the smallest adsorption capacity to CO₂ while the highest adsorption capacity to CH₄. In addition, in contrast to CH₄, the maximum CO₂ adsorption capacity of the clay minerals does not decrease continuously with the increase of pore size. As seen from Table 2, as pore size increases from 1 nm to 2nm, the maximum CO₂ adsorption capacity of montmorillonite, illite and kaolinite increases from 0.276 cm³/m², 0.211 cm³/m² and 0.179 cm³/m² to 0.374 cm³/m², 0.304 cm³/m² and 0.248 cm³/m², respectively. However, when the basal spacing further increases to 4 nm, the maximum CO₂ adsorption capacity of the clay crystals all show an obvious decreasing trend, indicating that the optimal storage space of CO₂ in clay minerals may be around 2 nm.

The effect of crystal types on CO₂ adsorption behaviors and the variation of the maximum CO₂ adsorption capacity with varying pore sizes are related to the characteristics of CO₂ and different crystal structures. CO₂ is a polar molecule with strong quadrupole moment (Jin and Firoozabadi, 2013; Deng et al., 2023). During the process of CO₂ adsorption, in addition to van der Waals interaction energy, the electrostatic energy between clay minerals and CO₂ also plays a key role. Figure 5 displays the contribution of different forces in interaction energy between the 2-nm clay minerals and CO₂ molecules under 20 MPa. As illustrated in Figure 5, both electrostatic energy and van der Waals energy contribute to the interaction energy. For montmorillonite and illite with cation exchange, the electrostatic energy contributes almost half of the total energy. Therefore, the limited pore volume of smaller pore (such as 1 nm) may limit the effect of strong electrostatic adsorption energy induced by CO₂ and charged clay crystals. As the basal spacing increases to 2 nm, the CO₂ adsorption capacity is further enhanced due to enough space, while when the pore size continuously increases to 4 nm, the larger width will weaken the overlapping effect of the two mineral surfaces, leading to the reduction of the adsorption capacity.

Furthermore, the stronger electrostatic interaction induced by CO₂ and the cation exchange in illite and montmorillonite can significantly promote the adsorption behavior of CO₂, especially at lower pressure. It is also seen from Figure 4 that compared with kaolinite, the total adsorption amount of CO₂ on montmorillonite and illite is obviously enhanced in 1∼3 MPa pressure range due to the interlayer cations. These interlayer cations can also greatly effect the morphology of CO₂ in the nanopore. Figure 6 presents the adsorption configuration of CO₂ molecules in 2-nm crystal pores at 20 MPa. It is found that in the montmorillonite model, the CO₂ molecules near the mineral surface are concentrated around Na⁺ cations. The distribution pattern of CO₂ molecules in illite is similar to that of CO₂ molecules in montmorillonite, in which the CO₂ molecules close to the walls of the crystal are oriented toward K⁺ cations. However, different from illite and montmorillonite, the CO₂ molecules in kaolinite close to the model surface are almost parallel to the crystal sheets, showing an orientation preference. In addition, hydrogen bonds are observed between some CO₂ molecules and the hydroxyl groups in kaolinite crystal model.

Based on the results in Tables 1, 2, the ratio of maximum adsorption capacity of pure CO₂ to CH₄ in different clay minerals is further compared and analyzed as illustrated in Figure 7. Obviously, the maximum adsorption capacity of CO₂ in different pore sizes of various clay minerals is always larger than that of CH₄. From the perspective of CS-EGR project, this conclusion is quite favourable. Specifically, the CO₂ adsorption capacity of clay minerals is almost 2–8 times than that of CH₄, and among these clay minerals the montmorillonite shows the best performance of preferential adsorption capacity, followed by illite and kaolinite. This conclusion is in line with previous experimental measurements (Kang et al., 2011; Heller and Zoback, 2014). According to the research by Heller and Zoback (2014) and Kang et al. (2011), for pure gases, the maximum adsorption capacity of CO₂ is almost 2–10 times greater than that of CH₄ in shale samples. Moreover, an in-depth analysis from Figure 7 reveals that the preferential CO₂ adsorption capacity of the clay crystals is further enhanced due to the increase of basal spacing, while the extent of improvement under different pore sizes exists discrepancies. For example, it is interesting to note that the adsorption capacity ratio of CO₂ to CH₄ shows an obvious improvement from 4.0 to 7.0 as the basal spacing rises from 1 to 2 nm, while the ratio slowly increases to 7.6 when the pore size further reaches to 4nm, implying that it is favorable for CO₂ to replace CH₄ molecules in smaller mesopore.

The GCMC simulation is also applied to further determine competitive adsorption behaviors of CO₂/CH₄ binary mixtures. Considering the differences in chemical and physical properties of the different clay minerals, the montmorillonite and kaolinite crystals are selected to compare the competitive adsorption mechanisms of the binary gas mixtures with and without the effect of cation exchange. Figure 8 shows the adsorption isotherms of gas mixtures (50.0 mol%:50.0 mol%) in 2-nm montmorillonite and kaolinite nanopore with increasing pressure. It is found that for montmorillonite the CO₂ isotherm shows a steep rise in low pressure range before approaching saturation, while a relatively slow rising trend is observed in the CO₂ isotherm for kaolinite at low pressure. Comparatively, CH₄ isotherms of both montmorillonite and kaolinite exhibit a gradual increase tendency with increasing loading. These different behaviors of isotherms are in accord with the performance of pure CH₄ and CO₂ in corresponding clay mineral. Furthermore, the adsorption amount of CO₂ in both kaolinite and montmorillonite is obviously higher than that of CH₄ in the whole adsorption pressure range, which is related to the interaction of charged clay structures and CO₂ molecules, especially the stronger electrostatic interaction induced by CO₂ and Na⁺ cations in montmorillonite. Such a large contrast in CO₂ and CH₄ adsorption amount is further favourable to implement the CS-EGR project. Moreover, in order to quantitatively study the selective adsorption capacity of binary gas mixtures in clay minerals, the adsorption selectivity S C O 2 / C H 4 is calculated based on the Eq. 4 as follows (Liu and Hou, 2020): S C O 2 / C H 4 = x c o 2 / x C H 4 y C O 2 / y C H 4 (4) where x_i is the molar fractions of component i in the clay pores, y_i is the molar fractions of component i in the bulk gas reservoir. Figure 9 depicts the adsorption selectivity of CO₂/CH₄ mixtures in 2-nm montmorillonite and kaolinite pores as a function of pressure. Notably, the values of adsorption selectivity in both montmorillonite and kaonilite are always much larger than one over the whole pressure range, confirming that the clay minerals have a stronger affinity for CO₂. Moreover, the montmorillonite exhibits larger selectivity than kaolinite under the same adsorption pressure, implying that montmorillonite can store more CO₂ with more CH₄ being recovered compared with kaolinate at the same condition. Similar results have been found by Liu and Hou (2020). It is observed from Figure 9 that the selectivity of CH₄/CO₂ in montmorillonite reduces sharply as the pressure increases from 0 MPa to 5 MPa, and then decreases slowly with further increase of pressure. It appears that for shale gas reservoirs rich in montmorillonite, lower pressure conditions are more favorable for CO₂ to replace CH₄. In other words, the depleted montmorillonite-rich formations are the best candidate for the implementation of CS-EGR. However, the selectivity of CH₄/CO₂ in kaolinite presents a slight fluctuation trend with increasing pressure, indicating that for kaolinite-rich reserviors the influence of pressure on the competitive adsorption behavior of CO₂/CH₄ mixtures is relatively insignificant.

Figures 10, 11 show the configuration snapshots of CO₂/CH₄ mixtures in the 2-nm montmorillonite and kaolinite nanopores, respectively. Interestingly, it is observed from Figure 10 that the CO₂ molecules preferentially occupy the adsorption sites near Na⁺ at lower coverage, forcing CH₄ molecules away from the wall surface. As a consequence, almost few CH₄ molecules can be found near the cation exchange, and the existing CH₄ molecules are randomly dispersed in the pore space. With increasing coverage, CO₂ molecules start to adsorb in the area away from the walls and coexist with small amount of CH₄ molecules. Finally, the CH₄ molecules are aggregated in the center of pore because of the extrusion effect of CO₂ molecules as the adsorption pressure increases to 20 MPa. On the other hand, although CO₂ molecules hold most of the adsorption sites in kaolinite surface, a small amount of CH₄ molecules can still adsorb near the crystal walls under different adsorption pressures as shown in Figure 11, indicating that kaolinite can provide a small amount of adsorption sites for CH₄ during the competitive adsorption process. With increasing loading, both CH₄ and CO₂ molecules gradually increase in proportion and distribute irregularly in the pore space. These changes in Figure 10 and Figure 11 can clearly reveal the essence for the differentially preferential adsorption amount of CO₂ over CH₄ in montmorillonite and kaolinite.