The N₂ adsorption-desorption data of the seven samples are plotted in Figure 2, depicting the volume of gas adsorbed at STP vs. relative pressure (P/P ₀ ). All adsorption-desorption hysteresis loops are of type IV, suggesting the mesopores’ dominance in the samples. Most samples show type H3 loops (Figure 2), indicating the presence of significant plate-like particles in the samples, resulting in narrow slit-shaped pores (Thommes et al., 2015). The adsorption-desorption profiles in all the samples coincide at extremely low relative pressure, except in the specimen CB-CCH-5, which is due to the presence of wedge- and slit-shaped pores (Li et al., 2019). Slopes of the adsorption curves changed at different segments of P/P ₀ , suggesting micropore filling at extremely low P/P ₀ , followed by monolayer and multilayer adsorptions at relatively higher P/P ₀ (Figure 2). The knee-bend in the adsorption curves indicates the completion of monolayer adsorption at P/P ₀ ≈ 0.3 in most of the samples. At P/P ₀ ∼0.5, most samples show sudden closures of desorption to adsorption limbs, while CB-A and CB-C show gentler closures, and CB-D shows none. This phenomenon is attributed to the tensile strength effect in which the desorption curves coincide with the adsorption curves (Thommes et al., 2015). The sudden collapsing of the hemispheric meniscus at P/P ₀ ∼0.5 represents the presence of pores smaller than 4 nm (Cao et al., 2015). This can also be verified from the mesopore modal width of CB-A and CB-C, which is around 14 nm, and CB-D around 36 nm (Table 2). All other samples have mesopore modal widths <4 nm (Table 2).

Micro-imaging (using FESEM: field emission scanning electron microscopy) also shows intragranular elongated and semi-circular slit-type pores are more prevalent in clay minerals (Figure 3A). Pyrites feature a few heterogeneous inter-crystalline pores, whereas intragranular phyllosilicates have complex, elongated, tapering pores (Figures 3B–D).

We used low-pressure CO₂ adsorption data to characterize <2 nm pores. The slopes of all adsorption curves are relatively higher at the beginning of adsorption, and it gradually decreases with increasing P/P ₀ , suggesting a higher initial adsorption rate, followed by a slower adsorption rate (Figure 4A). Analysis of the data reveals that micropore volume ranges from 0.007 cc/g (CB-J and CB-CCH-5) to 0.027 cc/g (CB-1–2), and Micropore modal width varies in a narrow range of 0.573–0.822 nm (Table 2).

We used adsorption arms of N₂ and CO₂ isotherms for the PSD analysis, as desorption arms limit the accuracy of the results in the 4–5 nm pore size region due to the tensile strength effect (Thommes et al., 2015). NLDFT model for CO₂ as adsorbate at 273 K has been used for estimating PSDs of micropores (<1.5 nm) and slit pores, whereas NLDFT equilibrium model-N₂ at 77 K evaluated the PSD of mesopore (1.5–35 nm). The PSD is limited to 35 nm as the NLDFT kernel for slit-type pores is limited to 35 nm. The PSD curves of all samples exhibit multimodal characteristics with several volumetric maxima (Supplementary Figure S4A, Figure 5). For a better understanding and quantitative assessment of the different pore architectures, we have applied the deconvolution method (Ulm et al., 2007) to determine the mean size and standard deviation of each pore family from their corresponding continuous pore size distribution curve (Figure 5). The distinct peak from PSD indicates the distinct pore size family. Normal (Gaussian) distribution describes the events regardless of whatever probability distribution describes the individual experiments. Here, we assume the pores could be divided into j=1 to n pore size families with sufficient contrast in pore size distributions. The jth pore family occupies a volume fraction f _j of the total porosity. The theoretical probability density function (PDF; P J x i , x ¯ j , σ j ) of a particular phase, which is assumed to fit a normal distribution, is defined as: P J x i , x ¯ j , σ j = 1 √ 2 π σ j 2 exp ⁡ ⁡ − x i − x ¯ j 2 2 σ j 2 (7) where x ¯ j and σ j are the mean value and the standard deviation of pore size distributions ( x i ) of phase j. Minimizing the difference between the data from the weighted model-phase probability distribution function (PDF) and the experimental PDF using the following equations, we can find the unknowns {f _J , x ¯ j , σ j }: min ∑ m ∑ i = 1 N ∑ j n f j P j   ( x i , x ¯ j , σ j − P x x i 2 (8) ∑ j = 1 n f j = 1 (9)

In Eq. 8, P x x i is the measured value of the normalized frequency of the pore size x i and m are the number of intervals (bins).

The deconvoluted pore size results of all the samples are presented as normal distribution curves with different colors corresponding to a particular pore family (Figure 5). The red dashed line is the fitted line to the solid green color experimental results. The fitting coefficients of all the curves are >0.889, suggesting a good fit between experimental and modeled data. The results show 3 to 4 families in the micropore region and 6 to 9 families in the mesopore region. The modal width of the micropore and mesopore families are ∼0.53, 0.62, 0.80, 1.12 nm, and ∼3.5, 5.1, 8.0, 14.6, 17.5, and 29.0 nm, respectively (Table 3).

We have analyzed the linear isotherm data and prepared the FHH plots (Supplementary Figures S5A,B) to understand the fractal dimensions (D _f ) of the pores (Table 4). There are two distinct linear segments in the FHH plot of all samples; one at 0<P/P ₀ < 0.5 (Region 1) and the other one at 0.5 < P/P ₀ <1 (Region 2) (Supplementary Figures S5A,B). We can verify this transitional change of fractal characteristics using the closures of the hysteresis loops at P/P ₀ ≈ 0.5 for all isotherm curves (Supplementary Figure S5A). Two fractal dimensions, D _f1 and D _f2, from the two linear segments of regions 1 and 2 (Supplementary Figure S5B) are calculated using Eq. 4 (Table 4). D _f1 are relatively low due to greater negative slopes, varying from 2.06 (CB-D) to 2.60 (CB-C). D _f2 are higher because of lower negative slope and lie within a narrow range of 2.63 (CB-CCH-5) to 2.67 (CB-C-B and CB-J) (Table 4).

We have applied multivariate PLS regression to find the dependency measure between independent and dependent variables in shale samples. We considered micro-, meso-, total pore volume, SSA, CO₂ uptake capacity, D _f1 , D _f2 , and D _s as dependent variables, whereas the mineral composition and TOC as selected to be independent variables. The PLS regression predicts the relationships between a set of dependent variables [Y] from a set of independent variables [X] when the number of dependent and independent variables is different (Geladi and Kowalski, 1986; Haenlein and Kaplan, 2004). This technique generalizes and combines features from principal component analysis and multiple regression model, defined as Y n × m = X n × p β p × m + e n × m (10) where β p × m is the regression coefficient matrix, e n × m is the error term, n is the number of observations, m is number of response variables, and p is the number of predicted variables. Multivariate PLS was applied to avoid the singular influence of one independent variable on predicting the response of the dependent variable. Dependency measures can be positive and negative depending upon the contribution of individual independent variables to predict the response model accurately. Independent variables showing positive and negative dependency measures suggest direct and inverse correlations. In the following subsections, we first evaluate and then discuss the dependency measure of pore parameters in our samples.

The PLS fitting results between independent and dependent variables show the TOC has a positive dependency measure with micro-, meso- and total pore volumes (Figure 7A, Supplementary Table S2). Quartz, Fe-bearing minerals, and clay have positive dependency measures with micropore volume, while negative dependency measures with mesopore and total pore volume (Figure 7A, Supplementary Table S2). Similar to the mesopore and total pore volume, the SSA, D _f2 , and D _s have a positive dependency measure with TOC and negative with quartz, Fe-bearing, and clay minerals (Figures 7B,C). CO₂ uptake capacity shows a positive dependency measure, while D _f1 has a negative dependency measure with all the independent variables (Figure 7C). In summary, TOC in Cambay shale strongly enhances the storage capacity by increasing SSA and pore volumes (micro-, meso-, total pore volume), eventually enabling higher sorption and free gas storage capacity. An increase in TOC also enhances pore structure complexity and roughness. CO₂ storage capacity shows a positive dependency measure with all the independent variables (Figure 7B). The Cambay shale, therefore, is a good site for CO₂ sequestration.

Deconvolution of the continuous pore size distribution by CO₂ NLDFT at 273K and slit pore, N₂ NLDFT at 77K model divides the pores into 13 pore families (Figure 5; Table 3). Families 1 to 4 belong to the micropore region, while families 5 to 13 belong to the mesopore region (Figures 8A-D; Table 3). We also checked the dependency measure of these pore families with the mineral composition and TOC (Figures 8C,D, Supplementary Table S3). In this analysis, the mineral composition and TOC are taken as independent variables, and the volume percent of each pore family as dependent variables (Supplementary Table S3).

Mesopore constitutes 95% of the total pore volume fraction, of which families 12 and 13 occupy more than 50% (Figures 8A,B and Supplementary Table S3). Family 12 and 13 show a strong positive dependency on quartz, Fe-bearing minerals, and clay (Figures 8C,D, and Supplementary Table S3). This indicates a higher concentration of quartz, Fe-bearing minerals (pyrite, siderite), and clay tends to increase the volume fraction of relatively larger pores (w ∼17.5 nm and 29 nm). While family 9, 10, and 11 show only positive dependency measures with TOC and negative dependency measures with other independent variables (Figures 8C,D and Supplementary Table S3). In the micropore region (families 1–4), families 3 and 4 constitute significant volume fractions (∼40–90% of micropore volume). They display a positive dependency measure with TOC and Fe-bearing minerals and a negative dependency measure with quartz and clay.

The micropore and early mesopore regions (up to a pore width of 3.60 nm) are primarily governed by TOC (Yang et al., 2016b). In contrast, the mesopore region (3.6–35.0 nm) is regulated by minerals such as Quartz, Fe-bearing minerals, and clay content (Figure 8E). The volume fraction of mesopore increases with hard minerals, including quartz, pyrite, feldspar, and dolomite. These minerals play a significant role in determining the size distribution and abundance of meso- and macropores (Liu et al., 2017).

Variation of pore parameters (SSA, micropore volume, mesopore volume, total pore volume, avg. pore width, mesopore modal pore width, and micropore modal pore width), clay content, TOC, EVRo, production index (PI), and CO₂ uptake capacity with depth is illustrated in Figure 9. EVR _o and PI show an increasing trend, while CO₂ uptake capacity and micropore volume decrease with depth (Figure 9). An increase in pressure and temperature with depth increases the vitrinite reflectance value (VRo) of kerogen; therefore, it increases the production index (PI). The pore parameters such as SSA, mesopore volume, total pore volume, average pore width, mesopore modal pore width, and micropore modal width do not show any overall trend with depth (Figure 9).

The porosity of shale decays exponentially with depth due to compaction (Magara, 1980). With the increase in depth, the diameter of macropores (>50 nm) formed in the intergranular space reduces due to overburden stress. In contrast, the modification of meso-and micropores that primarily occupy the intragranular space or the surface of the organic and inorganic matter of shale is negligible. Overall, the mesopores hardly altered with depth. We, therefore, conclude that depth has no or less significant role in controlling the mesopore architecture of the Cambay shales.

Fractal dimension of total pores D _s increases from 5.84% (CB-C) to 17.31% (CB-D) compared to the average surface fractal value (D _f ) of the accessible pore (Table 4). The lower surface roughness of accessible pores might be due to the smoothening of surface contact by fluids in the physisorption process or the fluid present in the source rock.

Dynamic processes like dissolution, precipitation, and diagenetic condition control the alteration of the pore surface fractal dimension in sedimentary rock-like shales (Aharonov and Rothman, 1996; Sen et al., 2002). Fractal dimension increases with an increase in diagenetic alteration. The alteration in fractal space is governed by reaction-limited growth and transport-limited growth (Aharonov and Rothman, 1996). In sedimentary rocks, ‘reaction-limited kinetics’ governs the growth (Nagy and Lasaga, 1992). The fractal dimension of a self-affine pore interface is related to the height (h(x)) evolved over the lateral extent of the surface interface (S). A statistical parameter can define the interface width (W), as W S = h x − h ¯ 2 1 2 (11) here, h ¯ is the mean height over the surface interface. For the self-affine surface, the interface width ( W S is related by a power law with linear dimension (S) of the substrate (Vicsek, 1992). W S ∼ S α (12)

Here, exponent α has simple relation with fractal dimension as: α = 3 − D (13)

We have generated 100×100 lattice space using simple discrete particle models of interfaces, roughening by deposition and dissolution variation (Figures 10A–F) (Aharonov and Rothman, 1996). Steps fill the original square lattice in a checkerboard manner, i.e., every filled site is a step of height 1 and surrounded by nearest holes by height 0. At every iteration of the deposition phenomenon, a block of height 2 is assigned to fill up randomly at local minima, and for the dissolution phenomenon, the subtraction of a block of height 2 is randomly chosen from the local maxima. At each iteration, a deposition event occurs with a probability of P + and a dissolution event with occur probability of P − ( P + + P − = 1 ).

Interfaces obtained for different deposition probability and their corresponding fractal dimension are shown in Figures 10G,H. The data reveal that deposition and dissolution probabilities greatly control the pore surface roughness. As the deposition probability increases, the surface becomes smoother (Figure 10H), while increasing dissolution probability yields pores with rough surfaces. D _s of Cambay shale has three modal clusters: 2.66 (CB-CCH-5), ∼2.78 (CB-1-2, CB-C, CB-C-B, and CB-D), and ∼2.87 (CB-A and CB-J) (Table 4), suggesting variable deposition and dissolution environments in Cambay shales. The deposition environment was prevalent when the deepest sample (CB-CCH-5, D _s = 2.66) was forged P + =1). The samples CB-1-2, CB-C, CB-CB, and CB-D (D _s ∼2.78) experienced dominating deposition and low dissolution environment (0.7 ≤ P + ≤0.8 and 0.2 ≤ P − ≤0.3). The samples CB-A and CB-J (D _s ∼2.87) were created in an environment of equal deposition P + and dissolution P − . This finding is consistent with previous studies by (Pandey and Dave, 1998; De et al., 2020), which focused on the Cambay shale deposited during a marine transgressive phase in the Palaeocene to Lower Eocene period. The Cambay shale is categorized into two formations: the Older and Younger Cambay Shale, with an erosional unconformity serving as the boundary between them. This unconformity promotes the dissolution process within the Cambay shale under marginal-marine to marine- depositional conditions.