Shale gas is trapped in shale formations with three statuses: adsorbed in the micropores and fine meso-pores, dissolved in brine and bitumen, and compressed in pores and fractures (Long et al. 2018; Jin and Firoozabadi 2016a; Qajar et al. 2015; Wang et al. 2020; Li et al. 2018a). These three kinds of methane together represent the gas-in-place (GIP), which is a significant parameter for evaluation of a shale reservoir and development of an extraction plan. The widely distributed nanoscale pores in shale formation and the high density of adsorbed methane can contribute up to 85% of total GIP (Xiong et al. 2021; Chen and Wang 2022; Curtis 2002). Compared with dissolved methane and compressed free methane, the understanding about adsorption capacity of shale is essential to both GIP evaluation and well productivity. On the other hand, accurate assessment of the adsorption behavior of shale gas will enable a better understanding of not only the transport of shale gas through nanopores (Li et al. 2016; Yuan et al. 2014; Xu et al. 2020; He et al. 2019), but also the thermodynamic characteristics of shale gas reservoirs (Chen et al. 2019; Zhou et al. 2019; Dang et al. 2020; Gao et al. 2020).

In the Sichuan Basin of southwest China, shale gas wells are commonly drilled to depths of around 4000 m, and the deepest shale gas well was drilled deeper than 7000 m, resulting in shale reservoir temperatures of up to 200 °C and reservoir pressures of up to 70 MPa, respectively. Obviously, the pressure and temperature of most shale formations exceed the critical pressure and temperature of methane, with corresponding values of 4.59 MPa and − 82.6 °C, and exceeds the measurement range of most isothermal adsorption instruments. The adsorption of methane in shale is classified as supercritical adsorption, which is characterized by unique physical properties distinct from its subcritical counterpart. Supercritical fluids lack distinct liquid and gas phases, as well as a saturated vapor pressure. As a result, the adsorption behavior of supercritical and subcritical gases differs significantly, posing a challenge in the study of supercritical methane adsorption in shale reservoirs. Furthermore, the adsorption potential of shale is subject to several parameters, including total organic carbon (TOC), thermal maturity, clay content, pore size distribution, temperature, and moisture content. The interplay of these factors makes the adsorption of shale gas a challenging endeavor.

Many attempts have been made to investigate the adsorption characteristics of various temperatures of methane on shale under in situ conditions. Table 1 lists experimental studies in recent years of methane adsorption on shale over a wide range of pressures and temperatures on one or more defined samples.

As summarized in Table 1, the gravimetric method (magnetic suspension balance) was mainly used for the experiments. The magnetic suspension balance in the gravimetric method enables the simultaneous measurement of the bulk phase density and excess adsorption by measuring the masses at three different positions under vacuum and pressure equilibrium conditions (Hu et al. 2022a), and eliminates error accumulation. This technique provides a significant advantage over other methods and was widely applied to gas adsorption and density measurements under a wide range of extreme and real-world conditions.

The measured adsorption all showed a maximum value and then decreased with increasing pressure due to the difference between excess adsorption and absolute adsorption in a supercritical adsorption system. The excess adsorption (n_e), i.e., the difference between the absolute adsorption (n_a) and the amount that would be present in the same volume at the density of the gas in the bulk phase (Pini 2014; Hu and Mischo 2022), can be written as follows according to the Gibbs definition:

where ρ_g is the density of the bulk gas, ρ_a is the average density of the adsorbed phase, V_a is the average volume of the adsorbed phase. Neither the ρ_a nor the V_a in Eq. (1) can be measured directly by experiment. As the equilibrium pressure increases, a phenomenon occurs where the rate of increase in the product of V_a and ρ_g surpasses the rate of increase in the number of absolute adsorption (n_a). Under such conditions, the excess adsorption reaches a maximum value before eventually declining.

The Langmuir- or D–R-based models were employed predominantly to characterize the experimental adsorption results. The Langmuir model is a widely used and popular tool for studying gas adsorption due to its simplicity and practicality in engineering applications. This model assumes monolayer adsorption and has proven valuable for many engineering purposes(Yang et al. 2017, 2015; Xiong et al. 2017; Tian et al. 2016; Ji et al. 2012; Xia et al. 2017; Shabani et al. 2018).

However, a key limitation of the Langmuir model is that it assumes that the pore surface of porous media is homogeneous and uniform, which does not accurately reflect the heterogeneous nature of most natural reservoirs. Furthermore, the Langmuir model fails to capture the change in adsorption heat that occurs as surface coverage increases, which is a crucial factor in many systems due to its inherent assumption of energy homogeneity (Yang et al. 2014). As a result, while the Langmuir model has been useful in certain contexts, more advanced models are necessary to account for the complexity of many natural systems. The limitations of the single-site Langmuir model for accurately describing the adsorption behavior of shale gas are widely recognized. Shale reservoirs consist of a complex mixture of kerogen and clay minerals, each with unique chemical and structural properties, leading to a highly heterogeneous pore surface (Chen et al. 2021b; Guan et al. 2020, 2022; Chiang et al. 2018; Liu and Zhu 2016). As a result, the adsorption energy at each location within the shale can vary significantly (Chen et al. 2023), rendering the single-site Langmuir model inadequate for accurately describing the adsorption behavior of shale gas in realistic systems. On the other hand, the relationship between adsorption characteristics and the temperature has not been thoroughly inspected in previous research. Typically, previous studies have independently fit adsorption isotherms at different temperatures and used temperature-independent parameters(Li et al. 2017, 2018b; Feng et al. 2020; Tian et al. 2016; Yang et al. 2015; Zhang et al. 2012; Song et al. 2018; Zou et al. 2017; Gasparik et al. 2014; Zhou et al. 2017). However, this method may overlook the potential link between adsorption characteristics and the temperature, which may result in reduced predictive accuracy. To predict the trend of adsorption characteristics more accurately with changing temperature, it is necessary to simultaneously fit both adsorption isotherms and the effect of temperature on adsorption characteristics, which can improve predictive accuracy and provide a better understanding of the physical and thermodynamic mechanisms involved in the adsorption process.

The conventional Langmuir model has been widely recognized for its success in describing high-pressure methane storage in gas-shales systems. However, a significant challenge remains in accurately estimating the density of the adsorbed phase. Up to now, direct measurement of the adsorbed phase density is not possible. The relative magnitude of the bulk phase density to the adsorbed gas density plays a crucial role in determining the shape of measured isotherms. The most common approach to address this issue is to treat the density of the adsorbed phase as an unknown fitting parameter that can be obtained by fitting it directly to the measured excess adsorption (Hu and Mischo 2020b). However, this approach yields densities of the adsorbed phase that fall within a wide range, sometimes up to 1027 kg/m³ for methane adsorption on dry shale (Tian et al. 2016) and up to 1011.2 kg/m³ for methane adsorption on wet shale (Krooss et al. 2017). These values are significantly higher than the density of liquid methane (422.36 kg/m³) at boiling point at atmospheric pressure, rendering the adsorbed phase density physically meaningless. An alternative approximation is to treat the density of adsorbed methane as a constant, such as the density of liquid methane at atmospheric boiling point (422.36 kg/m³) or the van der Waals density (373 kg/m³) (Tian et al. 2016; Yang et al. 2015; Sakurovs et al. 2007; Chen et al. 2021a; Gensterblum et al. 2013; Gai et al. 2020). However, recent investigation found that this method is not valid for all cases and can sometimes lead to significant declines on absolute adsorption as the equilibrium pressure increases (Zhou et al. 2018).

A third approach to estimate the density of the adsorption phase is to consider that at high pressures, methane adsorption reaches saturation, and the density of the adsorbed phase remains relatively constant(Pini 2014; Moellmer et al. 2011). With this in mind, the excess adsorption during the high-pressure stage can be assumed to be linearly related to the density of the bulk phase, with the absolute slope of this linear relationship representing the volume of the adsorbed phase. The density of the adsorbed phase can be obtained by determining the value of the intersection of the fitted line with the bulk phase density axis, which corresponds to the density of the adsorbed phase at zero excess adsorbed volume (Chen et al. 2021a; Cai et al. 2018; Gensterblum et al. 2010). To apply this method, the excess adsorbed volume and bulk phase density within the high-pressure region are measured and a linear relationship is fitted to the data. The value of the horizontal axis intercept of the fitted line represents the density of the adsorbed phase, as shown in Fig. 1. This method provides a more physically meaningful estimate of the density of the adsorbed phase compared to the previous approaches. But recent investigation has demonstrated that this method would overestimate the density of the adsorbed phase due to the swelling induced by adsorption (Zhou et al. 2018; Clarkson and Haghshenas 2013).

The schematic diagram of linear fitting method to determine the density of the adsorbed phase

To characterize the effect of heterogeneous shale surfaces and different adsorption energies on gas adsorption, the multi-site Langmuir model was proposed. The multi-site Langmuir model not only successfully fits the experimental data, but also characterizes the variation of isosteric heat of adsorption with increasing surface coverage (Tun and Chyun 2021), due to the non-homogeneity of the adsorption sites. Lu et al. applied a dual-site Langmuir model to describe the pressure dependence and temperature dependence of methane adsorption on Devonian shales (Lu and Li 1995). Stadie et al. (2013) showed that the dual-site and three-site Langmuir models achieved similar results, but were much better than the single-site Langmuir model, by investigating the adsorption of methane on zeolite-templated carbon over a wide range of temperatures and at pressures of up to 9 MPa. Tang et al., used the dual-Langmuir model to depict the methane adsorption on shale and carbon dioxide adsorption on coals with a wide range of temperatures and pressures (Tang et al. 2016, 2017a, b, 2019). Li et al. developed a multi-site Langmuir model taking heterogeneity of surface energy distribution and the temperature dependence of adsorbed gas density into consideration (Li et al. 2018c, 2019). Gao et al. (2020) integrated an empirical equation for the maximum adsorption capacity and applied this model to depict isothermal adsorption data from the gravimetric method (Feng et al. 2020). Although the high-precision magnetic suspension balance is advantageous in gravimetric methods for measuring both gas phase density and excess adsorption simultaneously (Hu et al. 2022a; Yang et al. 2020; Hwang et al. 2019; Ottiger et al. 2008), the authors choose to use the Redlich–Kwong equation of state to calculate the gas phase density instead. Lin et al., introduced the compressibility factor to modify the dual-site Langmuir model to achieve an accurate modeling of molecular simulation data under ultra-high pressure (up to 50 MPa) conditions (Lin et al. 2020). Yue et al. applied the dual-site Langmuir model to fit carbon dioxide adsorption data on coal up to 12 MPa and extrapolate the predicted adsorption amount to higher pressures and temperatures (Yue et al. 2019).

A careful review of the literature reveals that the adsorption phase density is either artificially set to a constant value (373 kg/m³ or 422.36 kg/m³) (Stadie et al. 2013) or set to an unknown parameter to be optimized (including temperature-dependent models) (Gao et al. 2020; Tang et al. 2017a; Li et al. 2018c; Lin et al. 2020) in the above dual-site and multi-site Langmuir models. Namely, the density of the adsorbed phase does not change with pressure, and the pressure dependence on the density of the adsorbed phase in this liquid-like system was ignored. Recent molecular simulations (Xiong et al. 2017; Lin et al. 2020; Huang et al. 2018; Pang and Jin 2019; Wu et al. 2019; Heller and Zoback 2014; Jin and Firoozabadi 2013; Zhao et al. 2018; Liu et al. 2018; Hu et al. 2022b) and simplified local density (SLD) models(Hu and Mischo 2020c; Qi et al. 2019a; Zeng et al. 2021; Pang et al. 2022, 2020; Guo et al. 2017; Miao et al. 2022; Huang et al. 2022) have provided evidence that the density of the adsorbed phase will increase with pressure during the adsorption process. As a result, the fixed density model fails to accurately recover excess adsorption, absolute adsorption, and thermodynamic parameters, due to the observed deviation in the density of the adsorbed phase with increasing pressure during the adsorption process. Several scholars have incorporated pressure-dependent adsorption phase density into the supercritical Dubinin–Radushkevich (SD–R) model to effectively characterize the adsorption behavior of methane on shale (Kong et al. 2021). However, the effect of thermal expansion on the density of the adsorbed phase has been overlooked in this model. Thermodynamic parameters, specifically the enthalpy change of adsorption (also referred to as the isosteric heat of adsorption, typically represented by a positive number) and the entropy change, are significant for shale gas research. The examination of thermodynamic parameters can provide additional insights into the mechanisms and processes involved in gas adsorption within shale formations, while also enabling the characterization of the heterogeneous nature of shale surfaces (Yang et al. 2022a; Hu et al. 2021; Gao et al. 2023; Duan et al. 2022). Furthermore, a precise evaluation of the thermodynamics of gas adsorption can contribute to a comprehensive understanding of the transport mechanisms of shale gas within nanopores (Afagwu et al. 2022). This knowledge can also establish a robust theoretical foundation for the advancement of separation and purification facilities dedicated to shale gas (Tun and Chyun 2021; Hamid et al. 2023). Direct and indirect are two commonly employed methods for determining the thermodynamic parameters of adsorption. The direct method involves measuring the heat released during adsorption directly using a calorimeter. However, due to its complexity and high cost, only a few studies have used the direct method to determine the enthalpy change of adsorption. Therefore, the most prevalent approach for determining the enthalpy change of adsorption as a function of adsorption quantity (loading) is the indirect method, which uses adsorption isotherms (Nuhnen and Janiak 2020). This paper aims to overcome this aforementioned limitation by using a temperature- and pressure-dependent density of the adsorption phase model. This approach is used to develop a dual-site Langmuir model and assess its effectiveness in estimating both thermodynamic parameters and absolute adsorption.

Although the single-site Langmuir model and the dual-site Langmuir model have been successful in analyzing supercritical gas adsorption, it is important to note that the former was developed based on a homogeneous surface and is unable to represent changes in adsorption heat, while the latter assumes a constant density of the adsorbed phase under each or all temperatures. However, numerous calorimetric experiments and numerical simulations have demonstrated that the adsorption heat is dependent on the loading and surface heterogeneity.

2.1 Temperature and pressure-dependent density of adsorbed phase

As previously mentioned, directly measuring the density of the adsorbed phase is difficult. However, recent molecular simulation (He et al. 2019; Xiong et al. 2017; Pang and Jin 2019; Pang et al. 2020; Xu et al. 2022; Tian et al. 2017; Mosher et al. 2013; Jin and Firoozabadi 2016b), the simplified local density model (SLD) (Pang et al. 2020; Miao et al. 2022; Hu and Mischo 2020d; Qi et al. 2019b; Yang et al. 2022b), the density function theory (DFT) (Sermoud et al. 2022), and indirect experimental measurement (Liu et al. 2022) offer some insights into the density of the adsorbed phase. After briefly reviewing previously published studies, we can conclude the dependence of the density of the adsorbed phase on temperature and pressure as follows: The density of the adsorbed phase increases with pressure increases but declines with temperature increases. To represent the monotonically increasing density of adsorbed phase with pressure, we use a function Langmuir form to describe the relationship between density of adsorbed phase and pressure (Liu et al. 2022).

It should be noted that the parameters in Eq. (2) do not have the same physical meaning as those in the conventional Langmuir model. Equation (2) is merely an empirical formula, and further research is needed to better understand the relationship between the density of the adsorbed phase and pressure. Moreover, as the temperature increases, the adsorbed phase with a liquid-like appearance will exhibit thermal expansion. Therefore, it is possible to employ an approximate method to describe the relationship between the density of the adsorbed phase and temperature. Taking the effects of pressure and temperature under consideration, the density of adsorbed methane can be written as:

where ρ_a is the molar density of adsorbed methane, mol/m³; T is the temperature, K; T₀ is the initial temperature, K; β is the coefficient of thermal expansion, K⁻¹; P is the pressure, Pa; and a, b, β are parameters that need to be fitted. It is important to reiterate that, unlike the original meaning of the Langmuir function, its parameters are unlikely to possess a direct physical interpretation and such an interpretation is not necessary (Mertens 2009).

The dual-site Langmuir model can be written as (Tang et al. 2017a):

The absolute adsorption and the coverage ratio θ can be expressed as follows:

where, n_e is the Gibbs excess adsorption, mol/kg; n_max is maximum adsorption capacity of shale, mol/kg; ρ_g is the bulk methane density, mol/m³; α is the weighting coefficient assigned to second adsorption site, dimensionless; and P is the equilibrium pressure, Pa. The equilibrium constant for the ith adsorption site, denoted by K_i(T), can be expressed using the Arrhenius equation:

where, A_i is the preexponential factor of the ith adsorption site, 1/Pa; E_i is the adsorption energy of the ith adsorption site, J/mol; R is the gas constant, 8.314 J/K mol; and T is the temperature, K. The term exp(−E_i/RT) can be used to represent the probability that a methane molecule, which is adsorbed on a surface, has sufficient energy to overcome the attractive potential of the surface (Lowell et al. 2013).

2.2 Data acquisition and processing

The data used in this study were measured by our previous study using a gravimetric apparatus (ISOSORP-HP Static II, Rubotherm GmbH, Germany). The ISOSORP-HP instrument demonstrates a precision with regards to temperature, pressure, and mass measurements, offering an accuracy of 0.01 °C, 0.01 bar (1 kPa) and 0.01 mg, respectively. The instrument’s temperature range spans up to 150 °C, while the maximum pressure is 35 MPa. Supercritical methane adsorption was determined at various temperatures (40–100 °C) and at pressures of up to 30 MPa. Samples NY09, NY11, NY17, and NY21, which were drilled at different depths from the Longmaxi formation (lower Silurian) in the Changning–Weiyuan region of Sichuan Province, China, were ground into grains for the high-pressure methane adsorption. Further details on the samples, including mineral composition, and pore characteristics, can be found elsewhere (Hu and Mischo 2020b).

In this study, the adsorption of methane was measured using two measurement points: Measurement Point 1 (MP₁), which included the weight of the sample container and shale sample and was measured at vacuum (MP_1,0), as well as under experimental conditions (MP₁(ρ,T)); and Measurement Point 2 (MP₂), which included the weight of the titanium sinker, sample container, and shale sample, and was also measured at vacuum (MP_2,0) and under experimental conditions (MP₂(ρ,T)) (Lin et al. 2020). The density of bulk methane inside the chamber was calculated using the expression (Ottiger et al. 2008).

In Eq. (8), m_sk,0 and m_sk represent the weight of the titanium sinker under vacuum and under experimental conditions, respectively, while V_sk denotes the known volume of the titanium sinker. The excess adsorption (n_e), i.e., the difference between the absolute adsorption (n_abs) and the amount that would be present in the same volume at the density of the gas in the bulk phase, can be expressed as follows, according to the Gibbs definition:

where, M is the mole mass of methane; m_s is the weight of the shale sample; and V₀ is the combined volume of the shale sample and sample container. The values of the V₀ and m_s were determined using high pressure helium gravimetry (Hwang and Pini 2019):

2.3 The enthalpy change, and the heat of adsorption

The entropy change that occurs upon adsorption, denoted as ΔS_ads, can be determined as a function of the absolute adsorption uptake (n_a) using the isosteric method, which is similar in principle to the Clausius-Clapeyron relation (Stadie et al. 2013):

where, S_a is the entropy of adsorbed phase; S_g is the entropy of the bulk phase, V_g is the molar volume of the bulk phase, V_a is the molar volume of the adsorbed phase, which is the reciprocal of the molar density of the adsorbed phase:

where ΔS_ads is the entropy change during the adsorption processing. The heat of adsorption which takes real gas behavior and the volume of the adsorbed phase into account can be expressed as:

where Q_st is the isosteric heat of adsorption and ΔH_ads is the enthalpy change during the adsorption processing. It is noteworthy that the Eq. (13) incorporates both the actual gas behavior and the volume of the adsorbed phase.

Assuming that the volume of the adsorbed phase is negligible (i.e., V_a ≈ 0), an alternative expression for the isosteric heat of adsorption can be derived.

When the gas is assumed to behave as an ideal gas and the volume of the adsorbed phase is considered negligible, a different formulation for the isosteric heat of adsorption can be derived.

when the gas is treated as an ideal gas and the volume of the adsorbed phase is not negligible, a distinct expression for the isosteric heat of adsorption can be derived.

The \(\left( {\partial P/\partial T} \right)_{{n_{\text{a}} }}\) can be expressed as a chain rule:

However, in many literatures, Eq. (21) has been erroneously written in the following form (Gao et al. 2020; Tang et al. 2016, 2017a, 2019; Tang 2016):

Alternatively, the enthalpy change can be calculated by using the integral form of the Clausius–Clapeyron equation, which involves plotting the natural logarithm of pressure (ln P) corresponding to absolute adsorption against the reciprocal of temperature (1/T) (Mabuza et al. 2022):

where C is constant induced by the integration. It is worth noting that the Clausius–Clapeyron equation ignores the effect of non-ideality of gas and the volume of adsorbed phase. Table 2 summarizes the different forms of the parameters K_i(T) and \(\left( {\frac{\partial P}{{\partial T}}} \right)_{{n_{\text{a}} }}\) in both the literature and this study.

2.4 Model verification

In this study, we obtained data from Tang (Tang 2016) and applied an improved dual-site model with a pressure and temperature dependent density of adsorbed phase to fit the data. The fitting process was carried out using a universal global optimization (UGO) method in the software of 1stOpt, allowing for simultaneous optimization of all model parameters. The parameters are listed in Table 3 and the fitted curves and experimental data are display in Fig. 2. Table 3 and Fig. 2 illustrate the superior performance of the improved model over the density constant model. This can be attributed to the consideration of the adsorbed phase density’s variation with respect to pressure and temperature, and the kinetic model can better depict the thermodynamics and interaction between methane molecules and shale surface. The root mean square error (RMSE) (Fitzgerald et al. 2005) of the improved model is only one-fifth of that reported in the original paper (0.0147 mmol/g) (Tang et al. 2017a).

The comparison between the experimental data (solid symbols) and the fitted results (dashed lines) from Tang (2016) and the improved model fitted results (solid lines) from this study

3.1 Results of bulk methane density

Accurately estimating the bulk methane density is crucial in both the volumetric method and gravimetric method for isotherm experiments. However, previous studies have relied on equations of state to calculate the bulk methane density, which can introduce significant experimental error due to the difficulty in accurately predicting the bulk phase density with these commonly used equations. Therefore, finding more accurate methods to determine bulk methane density is essential for improving the accuracy and reliability of experimental results. As shown in Fig. 3, the bulk methane densities exhibit a nearly linear relationship with pressure, while increases in temperature results in a decrease in methane density. The bulk methane density measurements obtained through the gravimetric method are in excellent agreement with the data available on the NIST webbook (https://webhook.nist.gov). The root mean square error (RMSE) values for temperatures of 313.15 K, 333.15 K, 353.15 K, and 373.15 K are 0.052 kg/m³, 0.041 kg/m³, 0.034 kg/m³, and 0.028 kg/m³, respectively. These findings demonstrate the accuracy and reliability of our experimental method for determining bulk methane density.

The experimental determined bulk methane densities and collected data from the NIST webbook

3.2 Results of excess adsorption and model fitting

Figure 4 illustrates the experimental determined excess adsorption on samples NY09, NY11, NY17 and NY21. All the excess adsorption All isotherms show a trend of increasing to a maximum and then decreasing, which is typical of methane adsorption on shale. Because physisorption is an exothermic process that occurs spontaneously, increasing temperature inhibits adsorption. From Fig. 4 and Table 4, the improved dual-site Langmuir model proposed in this paper can accurately characterize excess adsorption, and it is worth noting that all parameters of Table 4 are temperature independent. This means that the model can be used to predict the excess adsorption at higher temperatures and pressures. The RMSE of this mode ranges from 0.0015 to 0.0035 mmol/g, which indicates that this model is far superior to the original dual-site Langmuir model, which assumes a constant value for the adsorption phase density at all temperatures. The present study employs an assumption regarding the density of the adsorption phase as a function of pressure and temperature dependence. As a result, the proposed model enables a more precise characterization of the mechanism underlying methane adsorption on shale.

The experimental measured excess adsorption (dots) and dual-site Langmuir model fitted results (lines)

The absolute adsorption evaluated from the proposed dual-site Langmuir model and the original single-site Langmuir model are depicted in Fig. 5. The data from Fig. 5 share the following characteristics: (1) there is a positive correlation between the pressure and the predicted absolute adsorption, with the latter increasing as the former rises; (2) an increase in temperature had an adverse impact on the total quantity of methane adsorbed within the shale, and (3) the commonly employed single-site Langmuir model, which assumes a constant adsorbed phase density, is prone to underestimating the actual absolute adsorption capacity, and comparable findings have been reported in previous investigations(Hu and Mischo 2022; Tang et al. 2019). The adsorption isotherm surface (Fig. 5) possesses a distinct advantage in that it can be used to estimate gas content under subsurface reservoir conditions at greater depths with high pressure and high temperature, since the adsorption isotherm constitutes a component of the adsorption surface.

The 3D surface of absolute adsorption in various shales subjected to increased pressures and temperatures. The predicted absolute adsorption surface is represented by a color gradient ranging from blue to red, while the absolute adsorption calculated by the single-site Langmuir model is depicted by a color gradient ranging from white to dark grey

3.3 Thermodynamic analysis

3.3.1 The isosteric heat of adsorption for different models

Figure 6 illustrates the variation in enthalpy of sample NY09 at 313.15 K during the progression of adsorption, corresponding to the augmentation in surface coverage. Figure 6 clearly demonstrates that the commonly used incorrect method results in an underestimation of the adsorption heat by approximately 2.5 kJ/mol at low surface coverage.

The enthalpy changes of sample NY09 at 313.15 K, which are determined using different models. The solid lines represent the enthalpy changes calculated from the model proposed in this study, while the dashed lines indicate the enthalpy changes evaluated from Tang's method. The purple dotted line represents the compressibility factor of methane at 313.15 K

The isosteric heat of adsorption (Q_st,IG, the blue solid line in Fig. 6) is found to remain constant as surface coverage increase, based on the assumption that the volume of the adsorbed phase is negligible and the gas behaves ideally. However, this approximation does not account for any changes in volume of the adsorbed phase that may occur during the surface coverage or deviations from ideal gas behavior. As a result, the calculated Q_st,IG may be overestimated to a large extent.

The trend of the isosteric heat of adsorption (Q_st,RG) is depicted by the green solid line in Fig. 6. This line is generated using a method that accounts for the real gas behavior, and it shows a pattern of decreasing values as the surface coverage increases, followed by a subsequent increase. In Fig. 6, the trend of the isosteric heat of adsorption (Q_st,RG) is observed to be identical to the trend of the compressibility factor (Z), which is represented by a purple dotted line. The compressibility factor is the ratio of ideal gas density to real gas density of methane at 313.15 K. This suggests that the variations in Q_st,RG are influenced by the compressibility factor of methane.

In Fig. 6, the pink solid line represents the isosteric heat of adsorption (Q_{st, IG+VA}) determined using Eq. (16), which accounts for the volume of the adsorbed phase. This line shows a decreasing trend as the surface coverage increases, which is observed at both low and high surface coverage. At high surface coverage, a significant decrease in Q_st,IG+VA is observed instead of an increase, suggesting that the volume of the adsorbed phase has a stronger influence than non-ideal gas behavior. This indicates that the effect of the volume of the adsorbed phase predominates over the effect of non-ideal gas behavior at high surface coverage.

The isosteric heat of adsorption (Q_{st, RG+VA}, represented by the red solid line in Fig. 6), as determined from Eq. (13) that considers both the real gas behavior and the volume of the adsorbed phase, exhibits a declining pattern as the surface coverage increases. Based on the trend of Q_{st, RG+VA}, it can be concluded that both non-ideal gas behavior and the volume of adsorbed phase have a significant impact on the isosteric heat of adsorption. Neglecting either the non-ideal gas behavior or the volume of adsorbed phase will lead to an overestimation of the isosteric heat of adsorption.

It is also worth noting that at very low surface coverage (very low equilibrium pressure) of the isosteric heat of adsorption, calculated by considering the volume of the adsorbed phase (Q_{st, RG+VA} and Q_{st, IG+VA}), is 1 kJ/mol lower than the isosteric heat of adsorption calculated without considering the volume of the adsorbed phase (Q_{st, RG} and Q_{st, IG}). Nonetheless, the contrast was not detected in conventional Langmuir models, including the single-site and dual-site Langmuir models, as they assume a constant density of the adsorbed phase. Assuming a constant density of the adsorbed phase leads to an overestimation of the density of the adsorbed phase and an underestimation of the volume of the adsorbed phase at low surface coverage. This is due to the fact that the absolute adsorption is exceptionally low during this stage (Hu and Mischo 2022).

3.3.2 The isosteric heat of adsorption at various temperatures and pressures

At lower pressures, methane molecules exhibit a preferential adsorption tendency towards sites with high potential energy, which results in a greater amount of heat being released during the adsorption process. Subsequently, as the equilibrium pressure further increases, methane molecules gradually occupy sites with lower potential energy, leading to a decrease in the heat of adsorption, as displayed in Fig. 7. This can be explained by the increased availability of adsorption sites at higher pressures, which allows the molecules to occupy weaker adsorption sites with lower potential energy.

The −ΔH_ads of various shale samples under different temperatures and pressures

The isosteric heat of adsorption for all samples demonstrates significant variation from the low to high surface coverage ranges. This result is attributed to the fact that over 83.3% of the adsorption sites in these samples possessed a binding energy of greater than or equal to 13.23 kJ/mol. Conversely, less adsorption sites exhibited a binding energy less than 4.55 kJ/mol. Consequently, the disparity in binding energies resulted in a more noticeable alteration in the heat of adsorption for all samples from low to high surface coverage range. Additionally, the sample exhibiting the lowest adsorption capacity, NY17, demonstrated the highest isosteric heat of adsorption, which could possibly be attributed to its strong energetic heterogeneity. This phenomenon has also been reported in previous literature (Stadie et al. 2013).

On the other hand, the degree of variation observed in the heat of adsorption as a function of surface coverage can serve as an indication of the heterogeneity of adsorption sites. This is due to the heightened sensitivity of the heat of adsorption to shale microstructure in comparison to the adsorption isotherm. Distinguishing between adsorption sites with distinct energy levels is more readily achievable through analysis of the fluctuation in the heat of adsorption than through examination of the adsorption isotherm. This phenomenon arises due to the fact that for the single-site Langmuir model, the heat of adsorption remains constant as the surface coverage increases (Li et al. 2019). In the single-site Langmuir adsorption model, it is postulated that the binding energy of gas molecules at each site is uniform, resulting in a consistent heat release during the adsorption process. The sample NY17 exhibits the most prominent degree of heterogeneity in terms of the binding energy associated with its adsorption sites. Consequently, as the surface coverage of the sample increases, a marked reduction in the corresponding isosteric heat of adsorption is observed.

At elevated surface coverages, all samples show a slower decline in the isosteric heat of adsorption. In Eq. (13), the term -T(∂P/∂T)n_a exhibits a monotonic increase as a function of surface coverage, while the term (V_g − V_a) demonstrates a steep monotonic decline at low surface coverage, followed by a gradual decrease as the surface coverage increases. The minimum value of isosteric heat will be reached when the rate of increase of the term -T(∂P/∂T)n_a is balanced by the rate of decrease of the term (V_g − V_a), leading to a subsequent slower decrease in the heat of adsorption.

A rise in the isosteric heat of adsorption at high surface coverage was observed in the dual-site Langmuir model, which presumes a constant volume of adsorbed phases (Tang et al. 2019). The anomalous behavior observed can be attributed to the assumption of a constant volume for the adsorbed phase, which leads to overestimation of its actual volume at low surface coverage and underestimation at elevated surface coverage. Hence, it fails to capture the effect of the volume of adsorbed phase on the isosteric heat of adsorption and the complicated mechanism behind supercritical adsorption.

Figure 7 illustrates that the isosteric heat of adsorption at low coverage increases with temperature, indicating a positive correlation. However, at high surface coverage, the isosteric heat of adsorption decreases with increasing temperature. In contrast, disparate results were obtained using a temperature-dependent, single-site Langmuir model that accounts for the real gas behavior and the volume of the adsorbed phase (Yang et al. 2018).

3.3.3 The isosteric excess heat of adsorption

In general, adsorption experiments only allow for the measurement of excess adsorption, making it challenging to accurately determine the absolute adsorption and volume of the adsorbed phase directly from experimental results. To address this issue, an alternative method involves using excess adsorption to evaluate the “isosteric excess heat” of adsorption, as opposed to the “isosteric absolute heat” of adsorption (Salem et al. 1998). The isosteric excess heat accounts for real gas behavior, while disregarding the impact of the volume of the adsorbed phase. Mathematically, it is expressed as follows:

Figure 8 illustrates the isosteric heat of sample NY09 under various temperatures, estimated using both absolute adsorption and excess adsorption. As shown in Fig. 8, the isosteric heat calculated from absolute adsorption decreases as the equilibrium pressure increases. In contrast, the isosteric excess increases sharply, reaching an infinitely large value near 10 MPa, before turning to an infinitely negative value, which is unrealistic for physisorption. The infinitely large value is attributable to the denominator of Eq. (25) being equal to zero at the maximum of the excess adsorption. Moreover, the absence of differentiation between absolute and excess adsorption only allows for a rough estimate of the isosteric heat at very low pressures (below 2 MPa for sample NY09).

The isosteric excess heat and the isosteric absolute heat of adsorption as a function of the equilibrium pressure for sample NY09. The solid lines indicate the isosteric absolute heat calculated from absolute adsorption, while the dashed lines indicate the isosteric excess heat calculated from excess adsorption

In this paper, an improved dual-site model was established based on the pressure-dependent and temperature-dependent density of the adsorbed phase. The experimental data from our previous study and the literature were used to validate the feasibility of this model. The isosteric heat of the adsorption was discussed in detail. Through the theoretical analysis of this study, our conclusions are made as follows:

1. (1)

   Our approach provides a more robust representation of the experimental excess adsorption data across a wide range of pressures and temperatures. This approach works by assuming the density of the adsorbed phase increases in a manner analogous to the Langmuir function, with increasing equilibrium pressure and decreases in a manner analogous to the liquid expansion model with increasing temperature.

2. (2)

   The improved dual-site Langmuir model not only successfully represents the experimental excess adsorption data, but also provides a robust approach for extrapolating adsorption temperatures and pressures beyond the range of existing experimental instruments. Furthermore, a comparative evaluation between the absolute adsorption estimated using the improved model and the traditional model, which assumes a constant density of the adsorbed phase for fitting, demonstrates that the traditional method significantly underestimates absolute adsorption under high pressure.

3. (3)

   This paper derives the correct form of the equation for the partial derivative (∂P/∂T)n_a using a chain rule. Comparison of the isosteric heat indicates that the widely used incorrect form of (∂P/∂T)n_a will lead to an underestimation of the isosteric heat across the entire range of the equilibrium pressure. On the other hand, neglecting the real gas behavior, the volume of the adsorbed phase, and the energy heterogeneity of the adsorption sites would result in an overestimation of the isosteric heat cross the entire range of the equilibrium pressure.

4. (4)

   The isosteric heat, which considers real gas behavior, volume of the adsorbed phase, and energy heterogeneity of adsorption sites, decreases as the equilibrium pressure increases. At very low equilibrium pressure, the isosteric heat slightly increases with increasing temperature, but at high equilibrium pressure, it decreases significantly with increasing temperature.

5. (5)

   The isosteric heat, which is determined from the experimental excess adsorption data, is only an approximate representation of the actual isosteric heat at pressures below 2 MPa. A singularity is observed in the isosteric excess heat curves at the point corresponding to the maximum excess adsorption. However, in the case of shales buried at great depths, it is necessary to differentiate between these two types of adsorption heat.