Multifractal Methods in Characterizing Pore Structure Heterogeneity During Hydrous Pyrolysis of Lacustrine Shale

Abstract

By using gas physisorption and multifractal theory, this study analyzes pore structure heterogeneity and influencing factors during thermal maturation of naturally immature but artificially matured shale from the Kongdian Formation after being subjected to hydrous pyrolysis from 250 °C to 425 °C. As thermal maturity increases, the transformation of organic matter, generation, retention, and expulsion of hydrocarbons, and formation of various pore types, lead to changes in pore structure heterogeneity. The entire process is divided into three stages: bitumen generation stage (250–300 °C), oil generation stage (325–375 °C), and oil cracking stage (400–425 °C). During the bitumen generation stage, retained hydrocarbons decrease total-pore and mesopore volumes. Fractal parameters ΔD indicative of pore connectivity shows little change, while Hurst exponent H values for pore structure heterogeneity drop significantly, indicating reduced pore connectivity due to bitumen clogging. During the peak oil generation stage, both ΔD and H values increase, indicating enhanced pore heterogeneity and connectivity due to the expulsion of retained hydrocarbons. In the oil cracking stage, ΔD increases significantly, and H value rises slowly, attributed to the generation of gaseous hydrocarbons further consuming retained hydrocarbons and organic matter, forming more small-diameter pores and increased pore heterogeneity. A strongly negative correlation between ΔD and retained hydrocarbon content, and a strongly positive correlation with gaseous hydrocarbon yield, highlight the dynamic interaction between hydrocarbon phases and pore structure evolution. This study overall provides valuable insights for petroleum generation, storage, and production.

1. Introduction

Shale, as a crucial source of unconventional hydrocarbons, has a complex and heterogeneous pore structure that significantly affects the occurrence and mobility of oil and gas [1,2]. Lacustrine shale, due to its unique depositional environment and organic matter content as well as its importance in efficient and sustainable shale oil development in China, is an ideal and important subject for studying the changes in pore structure through thermal maturation. The heterogeneity of pore structure not only affects the generation, retention, and expulsion of oil and gas, but also has significant implications for reservoir enhancement and well productivity prediction. In recent years, with the continuous advancement of shale oil and gas exploration and development, methods for studying shale pore characteristics have become increasingly diverse [3,4,5,6]. Various 2D/3D visualization methods like optical microscopy, scanning electron microscopy, micro-3D X-ray microscopy, and focused ion beam scanning electron microscopy can characterize the surface porosity, pore type and geometry, and pore connectivity of shale [7,8,9,10,11]. Additionally, methods such as gas physisorption (GP; with CO₂, N₂, Ar), mercury intrusion porosimetry (MIP), and nuclear magnetic resonance (NMR) can evaluate pore structure parameters of shale reservoirs, such as pore volume and pore/throat size distribution [12,13,14,15,16,17,18,19,20,21]. Porosity and permeability can be determined using techniques like helium gas injection, fluid immersion, MIP, and NMR [22,23].

Thermal maturity is an important factor affecting the evolution of shale pore structure, with different types of kerogen and maturity stages having varying impacts on pore structure. Loucks et al. (2012) found that organic matter pores of Barnett shale were positively correlated with kerogen transformation [24]. Kuila et al. (2014) characterized various shales using GP techniques and found that the organic matter in low-maturity shale lacks micropores (<2 nm in diameters) and small pores, while these pores are observed in high-maturity shale organic matter [25]. Mastalerz et al. (2013) used a combination of GP and MIP techniques to characterize shale pores, and their results indicate that maturity is one of the important factors affecting the architecture of shale pores [26]. Its impact is comparable to that of organic matter content and mineral composition, but this pattern may be influenced by the characteristics of samples used in their studies. Due to the combined effects of various factors on geological samples, such as organic matter content, mineral composition, type of organic matter, burial depth, and depositional environment, it is difficult to analyze the impact of maturity on pore evolution in isolation. Therefore, some scholars use laboratory simulation methods to obtain different thermal maturity subsamples from the same initially immature sample to explore the impact of maturity on pore evolution. The main thermal simulation methods currently include gold tube thermal simulation, high-temperature and high-pressure closed system, and semi-open-system thermal simulation [27].

Compared to naturally evolved samples, the main advantage of simulated pyrolysis tests lies in its abilities to control a single variable when examining the pore development process. Although thermally matured samples obtained through high-temperature and high-pressure simulation tests cannot fully replicate geological conditions, this method allows us to study the specific effects of thermal processes on pore evolution. Throughout thermal evolution, various aspects of the pore structure change, including both geometric (pore size distribution, PSD) and topological attributes (pore connectivity) [28,29,30]. Throughout the thermal maturation of shale, the development of pore structure is influenced by various factors, including the alteration of organic matter, the creation and release of hydrocarbons, mineral changes, and organic–inorganic interactions [31,32,33,34,35]. Previous studies have shown that thermal maturation positively affects the pore structure of shale, with different types of organic matter or mineral compositions exhibiting a diverse evolution of pore structure during maturation [36,37]. Oil-producing shale often shows an increase in pore volume and the formation of large-diameter pores during the generation and expulsion of oil and gas, while gas-producing shale tends to form smaller-diameter pores [38]. However, there are few studies addressing the connection between pore structure heterogeneity and hydrocarbon generation/expulsion throughout thermal maturation, impeding the detailed assessment of shale reservoirs and the forecasting of shale petroleum presence.

Along with an integrated approach of combining various experimental techniques for a quantitative characterization of pore properties, the fractal theory has become an important tool for investigating the complexity of pore structures in coal, tight sandstone, and shale across different observation levels [39,40,41,42]. Although certain research employs a single fractal approach to assess the pore evolution traits of shale during thermal maturation, fractal dimension parameters are primarily used to comprehend the surface irregularities and structural complexities of geological samples, which is insufficient for a comprehensive quantitative assessment of pore network heterogeneity [43,44,45]. Through multifractal theory, the internal structure of pores can be characterized in great detail [46,47], as it segments the pore network into multiple units with varying singularities and quantitatively assesses the overall internal fine structure of porous media based on the fractal properties of different regions [48,49,50]. This research procured several artificially matured shale samples via hydrous pyrolysis tests and measured pore structure heterogeneity through NGP tests and multifractal analyses. Unlike the studies by Glema et al. (2010) and Manzari et al. (2014) that focus on anisotropy (plastic or damage) [51,52], this study focuses on using multifractal theory to analyze the heterogeneity of pore size distribution in shale at different pyrolysis temperatures, and explores its relationship with residual and expelled oil. Using originally immature lacustrine shale samples and their artificially matured sub-samples, a thermal maturation model was developed to explore the connections between organic matter conversion, hydrocarbon production–retention–release, pore structure alterations, and associated heterogeneity. This research connects variations in pore structure heterogeneity to the conversion of organic matter into hydrocarbons, which is crucial for identifying the key stages and priorities of petroleum accumulation. This not only aids in refining the exploration and development approaches for shale oil and gas assets, but also offers foundational support for the assessment of unconventional oil and gas reserves.

2. Samples and Methods

2.1. Samples

The samples used for pyrolysis were taken from the organic-rich laminated felsic shale of the Second Member of the Kongdian Formation in the Cangdong Sag of the Bohai Bay Basin (Figure 1). The total organic carbon (TOC) content is 2.27%, the vitrinite reflectance is 0.52%, the hydrogen index (HI) is 679 mg/g, and the kerogen is of type II1. The study area is characterized by typical continental sedimentary mudstone and shale, composed of interbedded laminated and massive shales. The organic matter is mainly of types I1 and II1, with TOC content ranging from 2% to 3%. The shale is buried at depths of approximately 2500 to 4300 m, with the vitrinite reflectance mainly distributed in the low to mature stages. Six small cores, measuring 2.5 cm in diameter and 5 cm in height, were extracted from the original samples for hydrous pyrolyses.

2.2. Experiments and Data Interpretation

2.2.1. XRD and Geochemical Analyses

Total organic carbon (TOC) analysis was conducted using a LECO CS744 instrument (LECO Company, St. Joseph, MO, USA). Rock pyrolysis analysis was performed using a Rock-Eval 7 Analyzer (VINCI Company, Paris, France), following standard procedures. The measured pyrolysis parameters included volatile hydrocarbon content (S₁, mg hydrocarbons (HC)/g rock), remaining hydrocarbon generation potential (S₂, mg HC/g rock), and the temperature of maximum pyrolysis yield (T_max). Quantitative analysis of mineral composition by X-ray diffraction (XRD) was carried out using a Bruker X-ray diffractometer (Bruker Company, Karlsruhe, Germany).

2.2.2. Hydrous Pyrolysis Experiments

This experiment was conducted at the Lanzhou Key Laboratory of Petroleum Resources Research, Chinese Academy of Sciences, using the WYMN-3 semi-open high-temperature and high-pressure thermal simulator. The experiment was designed with seven thermal simulation temperature points to obtain shale samples with different thermal maturities as well as their gaseous and liquid hydrocarbon products. Considering the impact of sample specifications on hydrocarbon generation and pore evolution, cylindrical samples with a diameter of 2.5 cm and a length of 5 cm were used in the experiment. Cylindrical samples better approximate the properties of formation samples, whereas granular samples may affect the evaluation of hydrocarbon generation and retention rates. To ensure the uniformity of the samples as much as possible, all cylindrical samples were drilled parallel to the bedding direction. The seven sets of semi-open high-temperature and high-pressure thermal simulation experiments were designed not only with temperature settings but also with corresponding fluid and lithostatic pressures calculated based on the burial history. The temperatures for the seven sets of samples were 250 °C, 300 °C, 325 °C, 350 °C, 375 °C, 400 °C, and 425 °C, with each experiment lasting 48 h and a heating rate of 10 °C/min. The lithostatic pressure was set between 45 MPa and 110 MPa, and the fluid pressure ranged from 18 MPa to 44 MPa, to enhance geological relevance and simulate natural processes at the Bohai Bay Basin. Following the completion of the experiments, the liquid products and solid residues were gathered for subsequent analyses, including the GP tests for PSD.

2.2.3. Low-Pressure CO₂ and N₂ Physisorption

Using the Micromeritics ASAP 2460 instrument (Micromeritics Company, Norcross, GA, USA), low-pressure CO₂ gas physisorption (CGP) at 273.15 K was conducted to determine the structure of micropores (<2 nm), and low-pressure N₂ gas physisorption (NGP) at 77 K was performed to analyze the structure of mesopores (2–50 nm) and macropores (>50 nm). Prior to the analyses, the samples were crushed to 35–80 mesh and degassed at approximately 110 °C for about 24 h to remove impure gases in the connected pore space. Subsequently, adsorption and desorption experiments were conducted under the relative pressures of N₂ or CO₂. The specific surface area (SSA) and the volume of mesopores and macropores were calculated using the Brunauer–Emmett–Teller (BET) model [52] and the Barrett–Joyner–Halenda (BJH) model [52], respectively.

2.3. Multifractal Theory

A multifractal analysis helps reveal the distribution patterns of pores at different scales, providing a deep characterization of the pore structure. Understanding the multifractal characteristics of shale is crucial. These characteristics can reveal changes in pore structure at different scales, aiding in accurate predictions of reservoir performance. This is essential for assessing the storage capacity and fluid flow paths, and improving extraction strategies of oil and gas reservoirs. To better understand the variability of pore structures, the multifractal analysis technique was used to describe the complexity of pore structures through two corresponding mathematical descriptions: the multifractal singularity spectrum a ∼ f(a) and the generalized dimension spectrum q ∼ D(q). In this work, we primarily used the q ∼ D(q) spectrum to assess the heterogeneity of PSD in GP experiments. The essential step of multifractal analyses is to divide the target object (length (L), range ([a, b])) into N(ε) boxes of scale ε, where (N(ε) = 2^k) (k = 0, 1, 2, …). In the i-th box, the mass probability function of the pore volume is defined as Pi(ε).

$$P_{\mathrm{i}}(\epsilon )={\mathrm{V}}_{\mathrm{i}}(\mathsf{\epsilon })/{\sum }_{i=1}^{N\left(\epsilon \right)}\mathrm{V}\mathrm{i}\left(\mathsf{\epsilon }\right)$$

where Pi(ɛ) is the mass probability function of the i-th box, and Vi(ε) is the pore volume of the i-th square. The distribution probability in the interval is satisfied in the case of scale variation:

$$p_i(\epsilon )~{\epsilon }^{-a_i}$$

a_i reflects the singularity strength representing the smoothness or regularity of data. For samples with multiple fractal features, N_a(ε) and ε are correlated as follows:

$$N_a(\epsilon )~{\epsilon }^{-f\left(a\right)},\epsilon \to 0$$

where f(a) is the multifractal spectrum as the fractal dimension of subsets with the same singularity index. The curve formed by a(q) and f(a) is called the multifractal singularity spectrum. Both a(q) and f(a) can be solved by the following equation:

$$a(q)\infty {\displaystyle \frac{{\displaystyle \sum _{i=1}^{N(\epsilon )}[u_i(q,\epsilon )\lg \epsilon ]}}{\lg \epsilon }}$$

$$f(a)\infty {\displaystyle \frac{{\displaystyle \sum _{i=1}^{N(\epsilon )}[u_i(q,\epsilon )\lg u_i(q,\epsilon )]}}{\lg \epsilon }}$$

where u_i(q, ε) is a family of probability measures, to be defined as:

$$u_i(q,\epsilon )={\displaystyle \frac{{p_i}^q(\epsilon )}{{\displaystyle \sum _{i=1}^{N(\epsilon )}{p_i}^q(\epsilon )}}}$$

where q represents the order of statistical matrix, taken as an integer between [−10, 10] with a step size of 1. When q >> 1, the information of high aggregation is amplified; when q << 1, the information of low aggregation is amplified. This work utilizes the approach and methodology included in the Excel 2007 spreadsheet calculations of Li et al. (2015) developed for MIP data [40].

For a given statistical moment order q, a function that satisfies $u(q,\epsilon )\infty {\epsilon }^{\tau \left(q\right)}$ is called a mass exponent function $\tau (q)$, which is a characteristic function of fractal behavior. In addition, the calculation formula for the generalized fractal dimension spectrum D(q) is:

$$D_q={\displaystyle \frac{1}{q-1}}\lim {\displaystyle \frac{\log \chi (q,\epsilon )}{\log \epsilon }}={\displaystyle \frac{\tau (q)}{q-1}}$$

3. Results and Discussion

3.1. Mineral Component Characteristics

Table 1. Mineralogy characteristics of the samples during hydrous pyrolysis.

Temp. (°C)	Whole Minerals (%)							Clay Minerals (%)
	Quartz	Potash Feldspar	Sodium Feldspar	Calcite	Dolomite	Zeolite	Clays	Illite	I/S	Chlorite
23	17.9	5.4	3.8	3.8	21.3	31.5	16.3	47	43	10
250	20.6	5.7	4.6	3.6	21.2	30.6	13.7	49	42	9
300	25.6	5.6	5.9	3.7	18.1	29.2	11.9	50	42	8
325	27.4	6.8	4.3	3.3	16.2	30.6	11.4	49	42	9
350	26.1	8.2	4.5	3.5	17.4	28.5	11.8	58	35	7
375	27.8	7.8	6.6	3.4	15.2	28.3	10.9	62	32	6
400	28.2	8.2	5.8	2.8	13.3	31.4	10.3	66	28	6
425	28.1	8.6	7.4	2.1	11.8	32.3	9.7	69	26	5

presents the mineralogical data of shale at different pyrolysis temperatures, including both whole-rock and clay mineral analyses. Overall, as the pyrolysis temperature increases, the content of quartz and feldspar minerals shows an upward trend, while the content of calcite, dolomite, and clay minerals shows a downward trend (Table 1). The changes in mineral content are influenced by multiple factors, among which the generation of organic acids may be a significant factor in the reduction of calcite and dolomite content. Since calcite and dolomite are carbonate minerals, they can react with organic acids, leading to a decrease in their content and the formation of dissolution pores, which can affect the pore structure.

Additionally, with the increase in pyrolysis temperature, illite and illite–smectite mixed-layer minerals also exhibit regular changes. The content of illite shows an increasing trend, while the content of illite–smectite mixed-layer minerals shows a decreasing trend. This change may be due to the transformation of illite–smectite mixed-layer minerals into illite, with the formation of organic acids providing potassium ions for the entire transformation process.

3.2. Organic Geochemistry

The original sample is in an immature stage with an R_o value of 0.56%. As shown in

Table 2. Geochemical parameters of shale samples after pyrolysis at different simulation temperatures.

Temp. (°C)	TOC (%)	S1 (mg/g)	S2 (mg/g)	S1 + S2 (mg/g)	HI (mg/g TOC)	OSI (mg/g)	PI	PC (wt.%)
23	2.27	0.76	15.39	16.15	679	33	0.05	1.39
250	2.13	0.51	12.6	13.11	591	24	0.04	1.16
300	2.08	0.58	12.12	12.7	584	28	0.05	1.12
325	2.3	1.05	12.65	13.7	550	46	0.08	1.21
350	1.68	1.2	6.92	8.12	411	71	0.15	0.73
375	1.06	0.62	2.51	3.13	235	58	0.2	0.29
400	0.61	0.2	0.11	0.31	19	33	0.64	0.05
425	0.55	0.19	0.13	0.32	7	16	0.71	0.03

and Figure 2, the initial TOC content of the samples is 2.27%. With the increase in thermal simulation temperature, the organic matter is consumed to 0.55%. As the hydrocarbon generation process gradually enters the mature and highly mature stages, the kerogen and hydrocarbons remaining in the samples are gradually consumed, and the hydrocarbon generation potential becomes increasingly smaller. Free hydrocarbons show a trend of first increasing and then decreasing, while pyrolytic hydrocarbons show a decreasing trend. These parameters are also controlled by hydrocarbon generation and expulsion.

When the thermal simulation temperature reaches 425 °C, the corresponding hydrocarbon generation stage is the extensive gas generation stage. At this stage, the contents of free hydrocarbons and pyrolytic hydrocarbons decrease to below 0.1 mg/g, indicating that the hydrocarbon generation process is very thorough. The production index increases rapidly at 350 °C, indicating an enhanced degree of conversion to hydrocarbons. At a thermal simulation temperature of 425 °C, the production index reaches a maximum value of 0.71. The HI is highest in the unheated state, gradually decreasing with the degree of evolution, and rapidly decreasing at 375 °C, reaching a minimum value at 425 °C. The organic carbon content also gradually decreases with the degree of evolution, reaching a minimum value at 425 °C. The OSI (oil saturation index) showed a tendency to increase and then decrease, reaching a maximum at 350 °C. The PI (productivity index) gradually increases with the degree of evolution, and increases rapidly at 400 °C. The PC (effective carbon) gradually decreases with the degree of evolution, reaching a minimum value at 425 °C.

As shown in Figure 2, overall, when the thermal simulation temperature reaches 350 °C, the samples enter the extensive hydrocarbon generation stage, and the HI decreases rapidly. When the temperature reaches 425 °C, the hydrocarbon generation potential is almost zero, indicating that the hydrocarbon generation and expulsion process is very thorough.

3.3. Hydrocarbon Generation

3.3.1. Retained Hydrocarbons in Different Occurrence States

Comparing the pyrolysis data before and after extraction (

Table 3. Quantification of retained hydrocarbons in samples before and after pyrolysis.

Temp. (°C)	Ro (%)	Extracted		Unextracted			Bound Oil (mg/g)	Bound/ Free Oil	Retained Oil (mg/g)
		S1	S2	S1	S2	TOC (%)
		(mg/g)	(mg/g)	(mg/g)	(mg/g)
23	0.51	0.17	12.73	0.76	15.39	2.27	2.66	3.50	3.42
250	0.53	0.20	9.89	0.51	13.18	2.13	3.29	6.45	3.80
300	0.58	0.16	8.32	0.58	12.77	2.08	4.45	7.67	5.03
325	0.79	0.18	5.58	1.05	11.82	2.3	6.24	5.94	7.29
350	1.08	0.16	2.72	1.20	5.43	1.68	2.71	2.26	3.91
375	1.46	0.14	0.74	0.62	1.64	1.06	0.90	1.45	1.52
400	1.75	0.16	0.05	0.20	0.17	0.61	0.12	0.60	0.32
425	1.99	0.11	0.02	0.19	0.13	0.55	0.11	0.58	0.30

), it was found that the post-extraction pyrolysis S₁ was almost zero, and S₂ also decreased partially. This indicates that the retained hydrocarbons include not only S₁ but also this extractable portion of S₂, which constitutes a significant proportion of the retained hydrocarbons. The changes in retained and bound hydrocarbons in the mudstone and shale as the thermal evolution maturity increases are shown in Figure 3. The content of retained oil ranges from 0.30 to 7.29 mg/g. With the increase in thermal simulation temperature, it first increases and then decreases, reaching the highest content at 325 °C. The ratio of bound hydrocarbons to free hydrocarbons increases with thermal maturity, first increasing and then decreasing, reaching a maximum value of 7.67 at 300 °C and a minimum value of 0.30 at 425 °C. From the initial state to 300 °C, due to the low thermal maturity, most of the generated hydrocarbons are bound to the organic matter. As it enters the oil and gas generation stage, the amount of free hydrocarbons increases. The amount of bound oil reaches a maximum value of 6.24 mg/g at 325 °C, then continues to decrease, reaching a minimum value of 0.11 mg/g at 425 °C. After 350 °C, both bound oil and free oil decrease with increasing temperature, but the rate of decrease for bound oil is higher than that for free oil. This indicates that bound oil is less likely to be preserved at this stage, and free oil is more easily accumulated in the shale at higher maturity levels.

Retained Hydrocarbon = S_{1 Unextracted} + S_{2 Unextracted} − S_{2 Extracted}

(1)

Bound Oil = S_{2 Unextracted} − S_{2 Extracted}

(2)

The remaining hydrocarbon generation potential (S2) of the extracted samples decreases exponentially with increasing maturity, with over 90% of the hydrocarbon generation potential being consumed within the oil window (Figure 4). At a thermal simulation temperature of 375 °C, the remaining TOC content is only 46.6% of the original sample, indicating that the consumption of TOC mainly occurs at relatively low maturity stages. When the simulation temperature is between 400 °C and 425 °C, the organic matter content decreases only from 0.61% to 0.55%, indicating that this stage primarily involves the cracking of liquid hydrocarbons into gaseous hydrocarbons.

3.3.2. Characteristics of Hydrocarbon Generation and Expulsion

During the high-temperature and high-pressure thermal simulation experiments in a semi-open system, organic-rich shales gradually reach the mature and highly mature stages, accompanied by the generation of oil and gas. As the amount of generated oil and gas increases and reaches the preset expulsion pressure, some of the gaseous and liquid products are expelled from the samples. In this study, the gaseous and liquid products generated during thermal evolution were collected and categorized into expelled oil, retained oil, and gaseous products based on their different states.

Table 4. Hydrocarbon generation yields during thermal maturation.

Temp. (°C)	Retained Oil (mg/g)	Expelled Oil (mg/g)	Liquid HC (mg/g)	Expelled Gas (mg/g)	Total HC (mg/g)
250	3.81	0.79	4.61	0.22	4.83
300	5.03	1.11	6.14	0.30	6.43
325	7.29	1.50	8.78	0.48	9.27
350	3.91	6.00	9.91	0.82	10.7
375	1.52	8.25	9.77	1.11	10.9
400	0.32	8.02	8.34	3.02	11.4
425	0.30	6.63	6.93	4.49	11.4

lists the yields of gaseous and liquid products collected at different thermal simulation temperatures, and calculates the total liquid oil yield and total hydrocarbon yield. As the thermal simulation temperature increases, the retained oil yield first increases and then decreases (Figure 5).

At 250 °C, the samples begin to generate hydrocarbons [36], with a retained oil yield of 3.81 mg/g. As the temperature increases, the samples gradually enter the oil generation stage, reaching a maximum retained oil yield of 7.29 mg/g at 325 °C. Subsequently, as thermal maturity continues to increase, some hydrocarbons retained in the pores begin to be expelled, or light hydrocarbons start to be generated and expelled when thermal maturity reaches a certain level. The retained oil yield approaches zero at 425 °C. The yield of expelled oil also shows a trend of initial increase and subsequent decrease, but the temperature at which it reaches the expulsion peak is 375 °C, later than the peak value of retained oil. This stage is mainly due to the increasing degree of thermal evolution, with organic matter gradually reaching the hydrocarbon generation peak, resulting in the expelled oil gradually increasing and reaching a maximum value of 8.25 mg/g.

After the experimental temperature of 375 °C, as the temperature increases, the samples gradually enter the wet gas generation stage. During this stage, some hydrocarbons are expelled in the form of gaseous products, leading to a decrease in expelled oil. The yield of gaseous hydrocarbons shows a trend of gradually increasing with temperature. This is due to the increasing maturity of the samples, with kerogen gradually reaching the mature and highly mature stages, leading to an increase in the generation of gaseous hydrocarbons. The yield of gaseous hydrocarbons increases from an initial value of 0.22 mg/g to a final value of 4.49 mg/g, with a rapid increase observed at 400 °C, as the organic matter enters the gas generation stage. The total hydrocarbon yield shows a gradual increase with the thermal simulation temperature. After the samples reach the oil generation stage, the total hydrocarbon yield increases rapidly and then slowly increases in the later stages, reaching a maximum value at 425 °C.

3.4. Pore Structure Characteristics

3.4.1. CGP and Pore Structure

According to the IUPAC standards, pore sizes are categorized into micropores (diameter < 2 nm), mesopores (2–50 nm), and macropores (>50 nm) [36]. Shale micropores are typically analyzed using the CGP method, as CO₂ can more effectively access micropores due to its higher experimental temperature compared to N₂. CGP analyses can characterize pore structures smaller than 2 nm. The thermal simulation process of hydrocarbon generation involves the stages of liquid and gaseous hydrocarbon formation. Pore spaces smaller than 2 nm are crucial for studying the occurrence and migration of gaseous hydrocarbons. Therefore, this work uses low-temperature CGP to characterize the pore evolution during the bitumen generation, oil generation, and wet gas generation stages. The initial CGP amount is 0.494 cm³/g, reaching a minimum of 0.276 cm³/g at 325 °C, and then continuously increasing to a maximum of 1.273 cm³/g at 425 °C (Figure 6). The CGP pore volume distribution ranges from 0.24 mm³/g to 1.37 mm³/g, with an average value of 0.65 mm/g. It reaches a minimum at 325 °C and a maximum at 425 °C, showing an overall trend of first decreasing and then increasing. From the initial stage to 325 °C, the adsorption amount decreases with increasing temperature, possibly due to pore blockage caused by the generation of bitumen or heavy oil. From 350 °C to 375 °C, the samples enter a stage of significant oil generation, accompanied by the expulsion of oil and gas. During this stage, organic matter transforms into oil and gas, creating organic pores and increasing pore volume. From 400 °C to 425 °C, organic matter further matures thermally, entering the wet gas generation stage. Hydrocarbons retained in the pores further crack into gaseous hydrocarbons, leading to an increase in pore volume.

Figure 7 shows the CGP-derived PSD of solid residues from the hydrocarbon generation thermal simulation of the Second Member of Kongdian Formation at different thermal maturity stages. Overall, the PSD shows three peaks at 0.6 nm, 0.8 nm, and 1.3 nm. The peak values gradually decrease from the initial sample to 325 °C and increase from 325 °C to 425 °C. The main peak appears at 1.5 nm, with a significant proportion of pore volume in the 1–2 nm range. Figure 7 shows the distribution curves of specific surface area and pore size. The specific surface area curve also shows three peaks at 0.6 nm, 0.8 nm, and 1.3 nm. Compared to the pore volume distribution curve, the peaks at 0.6 nm and 0.8 nm are higher, indicating that micropores contribute more to the specific surface area than to the pore volume. The specific surface area ranges from 0.41 m²/g to 2.36 m²/g, with an average value of 1.19 m²/g. The specific surface area changes with the thermal simulation temperature, showing a trend of first decreasing and then increasing, reaching a minimum at 325 °C and then increasing to a maximum at 425 °C.

3.4.2. NGP Analyses and Pore Structure

Figure 8 shows the NGP adsorption and desorption profiles of the initial shale and the pyrolyzed samples. Based on the IUPAC hysteresis loop classification, the pyrolyzed samples are categorized as H3 type, indicating the prevalent presence of slit-shaped pores. As depicted in Figure 8, when the relative pressure is below 0.6, the adsorption capacity increases gradually, suggesting that the micropores and small mesopores in the samples contribute minimally to the pores measurable by the NGP method. The pore volume and SSA data of mesopores and macropores were obtained through NGP analysis, as their detectable pore size range (2 nm to 217 nm) is broader than that of CGP. As shown in

Table 5. Pore volume and SSA results at different simulation temperatures by CGP and NGP.

Sample	CGP-Micropores		NGP-Mesopores		NGP-Macropores		(CGP + NGP)-Total
ID	Volume	SSA	Volume	SSA	Volume	SSA	Volume	SSA
	(10−3 cm3/g)	(m2/g)	(10−3 cm3/g)	(m2/g)	(10−3 cm3/g)	(m2/g)	(10−3 cm3/g)	(m2/g)
23 °C	0.456	1.307	21.9	3.54	7.14	0.173	29.59	5.02
250 °C	0.475	0.791	20.0	3.35	6.96	0.166	27.4	4.31
300 °C	0.311	0.538	16.3	2.46	7.99	0.181	24.6	3.18
325 °C	0.238	0.405	16.1	2.26	10.3	0.230	26.7	2.89
350 °C	0.344	0.587	14.4	1.94	12.0	0.253	26.8	2.78
375 °C	0.853	0.857	15.4	2.04	12.6	0.267	28.8	3.17
400 °C	1.176	1.64	15.4	1.85	16.2	0.342	32.8	3.83
425 °C	1.367	2.359	17.9	2.04	20.1	0.438	39.4	4.84

, the mesopore volume of the samples ranges from 14.445 to 21.868 10⁻³ cm³/g, and the specific surface area ranges from 1.84 to 3.54 m²/g. The macropore volume ranges from 6.96 to 20.1 × 10⁻³ cm³/g, and the specific surface area ranges from 0.166 to 0.438 m²/g.

The pore volume and specific surface area parameters of different pore types in the samples exhibit certain evolutionary characteristics with increasing simulated temperature. The total pore volume first decreases and then increases, reaching a minimum value of 24.6 × 10⁻³ cm³/g at 300 °C. This may be due to the bitumen generated by the shale before 300 °C blocking the pores, leading to a reduction in pore volume. As the simulated temperature increases, the organic matter further matures, generating and expelling more hydrocarbons, resulting in an increase in total pore volume. The mesopore volume of the samples decreases gradually and then increases with rising temperature. The mesopore volume is a significant contributor to the total pore volume of the shale samples, accounting for 45.5% to 74.2% of the total pore volume. The changes in mesopore volume exhibit a similar pattern to the total pore volume (Figure 9). The macropore volume of the samples gradually increases with rising temperature. This is due to the formation of new pores during the hydrocarbon generation process or the further expansion of existing small pores, leading to an increase in macropore volume. The total specific surface area and micropore specific surface area of the samples exhibit similar properties with temperature changes, showing a trend of first decreasing and then slowly increasing. The initial decrease in specific surface area is influenced by the reduction in micropores, while the later increase is related to the development of micropores and mesopores.

The PSD curves of the original samples and the thermally simulated samples are shown in Figure 10. The curves exhibit a bimodal pattern, with the main peak located at 40–100 nm and another lower peak at 3 nm. According to the graph, the peak at 3 nm in the unheated samples is significantly higher than that in the thermally simulated samples, indicating that the original samples have more small pore sizes. Overall, as the temperature increases, both peaks first decrease and then increase, with the right peak shifting towards larger pore sizes. When the temperature rises to 250 °C, the PSD curve changes significantly. The peak at 3 nm becomes less pronounced, and the curve in the 3–30 nm range drops sharply, indicating a rapid decrease in pores of this size, which is related to the generation of bitumen blocking the shale pores. The specific surface area and PSD curves show that the micropore specific surface area decreases sharply from the original sample to 250 °C. When the temperature rises to 300 °C, the PSD curve further declines, and the pores in the 3–40 nm range continue to decrease. Additionally, the specific surface area and PSD curves show that the specific surface area in the 3–20 nm range is also decreasing. The reason for the fewer pores of this size may still be related to the generation of bitumen blocking the pores. When the temperature rises to 325 °C, the PSD curve (Figure 9) shows that the pore volume in the 20–100 nm range increases to a certain extent, while the change in the left peak is not very obvious. This is due to the generation and expulsion of hydrocarbons at this stage, promoting the development of relatively larger pores. When the temperature rises to 350 °C, the PSD curve becomes unimodal, the left peak basically disappears, and the right peak continues to rise, with more pores larger than 30 nm developing. When the temperature continues to rise to 375 °C, the left peak re-forms and increases compared to 350 °C, the right peak shifts relatively to the right, and the pores in the 3–200 nm range increase compared to 350 °C. This stage is in the peak period of hydrocarbon generation and expulsion, which is conducive to pore formation. The specific surface area PSD curve shows an increase in the 3 nm peak, indicating the development of small pore sizes at this stage. When the temperature rises to 400 °C and 425 °C, the pore volume and PSD curves are still bimodal, with the left peak at 3 nm and the right peak at 100 nm. The curve peak at 100 nm rises significantly, and at this stage, the yield of gaseous hydrocarbons begins to increase significantly, accompanied by the expulsion of liquid hydrocarbons. The generation of gaseous hydrocarbons causes the remaining hydrocarbons in the pores to further pyrolyze, which is conducive to the increase in pore volume. Overall, at different simulated temperatures, the pore development shows the following patterns: from the initial state to 300 °C, the pores in the entire pore size range decrease; from 325 °C to 375 °C, pores larger than 20 nm develop significantly, while pores around 3 nm do not develop significantly; from 375 °C to 425 °C, pores at 3 nm increase, and pores larger than 20 nm further develop. The pore change patterns are consistent with the stages of hydrocarbon generation and expulsion: the initial to 300 °C is the bitumen generation stage, 325 °C to 350 °C is the oil generation stage, and 375 °C to 400 °C is the wet gas generation stage.

3.4.3. NMR Characteristics of Samples with Different Thermal Maturities

Nuclear magnetic resonance (NMR) experiments are widely used to characterize shale pore features due to their ease of operation, rapid testing, and non-destructive nature. In this study, the NMR signals of dry and water-saturated samples were tested, and the signals of the dry samples were subtracted to obtain the pore distribution characteristics of the shale. NMR utilizes the response characteristics of atomic nuclei to magnetic fields, with hydrogen nuclei (¹H) and carbon nuclei (¹³C) being commonly used. In this experiment, the samples were saturated with water, and the response characteristics of hydrogen nuclei (¹H) in the water were used to characterize the pore properties. The NMR signals generated by fluids in pores of different sizes correspond to different relaxation times. The relaxation time is positively correlated with the pore size, meaning that the larger the relaxation time, the larger the corresponding pore size. As shown in Figure 11, the T₂ spectra of the original and thermally simulated samples exhibit a bimodal distribution, with peaks at 1–10 ms and 50–90 ms. The NMR signals show a certain regularity with temperature changes. First, from the original sample to 325 °C, the main peak of the sample decreases, which may be related to pore blockage caused by bitumen generation, as previously discussed. From 325° C to 375 °C, the NMR signal intensity increases, accompanied by a right shift of the left peak and a slow rise of the right peak. These changes may be due to the significant generation of liquid hydrocarbons, leading to an increase in pores. From 400 °C to 425 °C, the two peaks of pore size further increase, and the corresponding pores also increase. Overall, the results obtained from NMR are consistent with the findings from NGP analyses.

3.5. Responses of Multifractal Parameters to Simulation Temperatures

3.5.1. Multifractal Characteristics

Organic-rich shales undergo the consumption of organic matter and the formation of hydrocarbons and organic acids during hydrous pyrolysis, thereby altering the physicochemical properties of the samples. These changes affect the pore structure, which undergoes corresponding variations at different pyrolysis stages, further leading to fluctuations in fractal dimensions. The results of the NGP analyses (

Table 6. Multifractal parameters of the original sample and artificially matured shale samples pyrolyzed at 250–425 °C.

Temperature (°C)	D1	D2	D−10	D10	a−10	a10	a0	amin−amax	Dmin−Dmax	H
23	0.5513	0.3486	1.8528	0.2077	2.0120	0.1869	1.4809	1.8251	1.6451	0.6743
250	0.4504	0.2350	1.5860	0.1341	1.6908	0.1207	1.4952	1.5701	1.4519	0.6175
300	0.4188	0.1996	1.4858	0.1122	1.6157	0.1010	1.4655	1.5147	1.3736	0.5998
325	0.3838	0.1722	1.3488	0.0964	1.4642	0.0868	1.4620	1.3774	1.2523	0.5861
350	0.4071	0.1861	1.3768	0.1039	1.4907	0.0935	1.4484	1.3971	1.2729	0.5930
375	0.4009	0.1887	1.4926	0.1058	1.5918	0.0952	1.5104	1.4965	1.3867	0.5943
400	0.3948	0.1907	1.7531	0.1072	1.9210	0.0965	1.5759	1.8245	1.6459	0.5953
425	0.3999	0.1959	1.8211	0.1101	1.9957	0.0991	1.5905	1.8966	1.7111	0.5979

; Figure 12) indicate that the q ∼ D(q) spectra of the samples before and after hydrous pyrolysis exhibit an inverse S-shape. The generalized dimension D(q) decreases monotonically within the q (−10~10) range, indicating that the shale pores have distinct multifractal characteristics. The pore structure exhibits varying complexities at different pyrolysis temperatures. The multifractal q ∼ D(q) characteristic parameters include the capacity dimension (D₀), information dimension (D₁), correlation dimension (D₂), and the dimensions of the minimum/maximum probability subsets (D_min and D_max). The Hurst exponent (H), defined as (D₂ + 1)/2, can be used to describe pore connectivity, with higher H values indicating better pore connectivity. ΔD, defined as the difference between D_min and D_max, is an important parameter for describing pore structure heterogeneity. The D₁ distribution ranges from 0.384 to 0.551. ΔD decreases initially and then increases with rising pyrolysis temperature, ranging from 1.25 to 1.71. ΔD reaches its minimum and maximum values at 325 °C and 425 °C, respectively, indicating that as the pyrolysis temperature increases, the pores become more uniform initially and then more concentrated, with the pore structure being most heterogeneous at 425 °C. During pyrolysis, the H value shows a trend of first decreasing and then increasing, ranging from 0.586 to 0.674, with the lowest value at 325 °C. This indicates that pore connectivity is poorest at 325 °C, and subsequently improves.

3.5.2. Relationship Between Pore Structure and Multifractal Dimension During Pyrolysis

With the increase in thermal maturity, the transformation of organic matter, the generation, retention, and expulsion of hydrocarbons, as well as the diversity in pore formation, collectively drive the changes in pore structure heterogeneity [53,54,55]. To understand the relationships between geochemical characteristics, thermal maturity, pore structure, and its heterogeneity, Pearson correlation analysis can be applied (Figure 13). The correlation between ΔD and SSA, mesopore, and macropore volumes is relatively weak, but it shows a strong correlation with the total volume. This indicates that during the hydrothermal pyrolysis experiment, the total pore volume significantly affects the heterogeneity of the pore structure. Parameters H and D₁ show a strong correlation with mesopore volume, possibly because hydrocarbons are not effectively expelled during the initial stage of pyrolysis, occupying smaller pores. During the extensive hydrocarbon expulsion stage, residual hydrocarbons in the mesopores are expelled or cracked to generate gas, further increasing the mesopore volume. The correlation between macropore volume and the H value is relatively weak and shows a certain negative impact, indicating that the influence of macropore volume on connectivity is not as strong as that of mesopores.

3.5.3. Pore Heterogeneity Variation and Multifractal Dimension During Hydrous Pyrolysis

During hydrous pyrolysis, the physicochemical properties of shale samples undergo significant changes, directly affecting the fractal characteristics of the pore structure. Simultaneously, the generation and expulsion of hydrocarbons also have a significant impact on the complexity of the pore structure. Based on the characteristics of hydrocarbon generation and expulsion during the hydrous pyrolysis process and organic geochemical data, the entire evolution process can be divided into three stages: bitumen generation stage (~300 °C), oil generation stage (325–375 °C), and oil cracking stage (400–425 °C). By comparing the PSD curves before and after extraction, the spatial distribution of retained hydrocarbons can be determined. The retention pattern of hydrocarbons is closely related to the thermal maturity stage. Firstly, during the heavy hydrocarbon generation stage (250–325 °C) [56,57,58], more hydrocarbons are retained in the pores. Figure 14 shows that the retention space of hydrocarbons gradually increases, and the newly formed pores exhibit a trend of increasing pore size. At a simulated temperature of 325 °C, the pore size range for retained hydrocarbons is 3–200 nm. When the temperature rises to 350 °C, the sample enters the oil generation stage, and a large amount of hydrocarbons is expelled from the sample, especially those retained in the 40–200 nm pores. At this temperature stage, more hydrocarbons are retained in pores smaller than 40 nm. Under conditions of 400–425 °C, the amount of retained hydrocarbons in the sample is only about 0.3 mg/g, with an expulsion efficiency exceeding 95%. The PSD curves of the samples before and after extraction are almost completely overlapping.

The correlation between ΔD and the amount of hydrocarbon generation shows a strong negative correlation with retained hydrocarbons and a strong positive correlation with gaseous hydrocarbons (Figure 15). This may be due to retained hydrocarbons blocking some small pores, making the pore distribution more uniform. The generation of gaseous hydrocarbons reopens the blocked small pores and creates new small pores, leading to a more uneven pore distribution. During the bitumen generation stage (~300 °C), TOC and HI show little change, the total oil generation is low, and retained hydrocarbons dominate. The pore structure shows a decrease in total pore volume and mesopore volume, with little change in macropore volume. The multifractal parameter ΔD shows little change, while the H value decreases significantly, possibly due to bitumen generation blocking some pore channels, reducing pore connectivity.

In the early oil generation stage (300–325 °C), ΔD decreases significantly. The pore structure transitions from high heterogeneity to low heterogeneity due to the reduction of small pores and the increase of large pores. Small pores are occupied by retained hydrocarbons, and some pores enlarge due to the consumption of organic matter during hydrocarbon generation. Meanwhile, retained hydrocarbons block pore channels, leading to poor connectivity. During the peak oil generation stage (350–375 °C), both ΔD and H values increase, indicating increased pore heterogeneity and connectivity. This is attributed to the expulsion of retained hydrocarbons, which releases small pores and opens pore channels. In the oil cracking stage (400–425 °C), ΔD increases significantly, and the H value increases slowly. At this stage, the content of retained hydrocarbons is very low, while the production of gaseous hydrocarbons increases significantly (Figure 16). The reduction in retained hydrocarbons can be attributed to oil cracking and gas generation. The expulsion of liquid hydrocarbons during the high oil generation stage significantly increases pore connectivity, which further increases during the gas generation stage of oil cracking. During the gas generation stage, a large number of small pores develop, leading to a significant increase in pore heterogeneity.

4. Conclusions

By combining multifractal theory with gas physisorption results, this work studies the evolution of pore architecture and heterogeneity in naturally immature and artificially matured shale samples. Thermal maturity influences the transformation of organic material, the processes of hydrocarbon generation, retention, and expulsion, as well as variations in pore architecture—key factors in petroleum formation, storage, and extraction. This study conducted semi-closed-system hydrous pyrolysis experiments on low-maturity shale samples from the Kongdian Formation in the Cangdong Sag of the Bohai Bay Basin, within the temperature range of 250 to 425 °C, revealing the relationships between organic matter thermal evolution, generated oil content, pore structure characteristics, and multifractal parameters.

The contents of quartz, feldspar, and illite minerals increase with the rise in thermal simulation temperature, while the contents of calcite, dolomite, and clay minerals show an overall decreasing trend. The TOC content and hydrocarbon generation potential decrease with an increasing temperature, while pore volume first decreases and then increases. Overall, the distribution of retained hydrocarbons is mainly controlled by the thermal evolution stage. During the heavy hydrocarbon generation stage, the generated hydrocarbons are retained within the pores, and the range of retained pore sizes gradually increases. After entering the oil generation stage, the expulsion efficiency of retained hydrocarbons increases, and hydrocarbons retained in larger pores gradually disappear. Upon entering the gas generation stage, the amount of retained hydrocarbons is extremely low, and from the PSD curves before and after extraction, almost no hydrocarbons are present.

Multifractal analyses indicate a significant correlation between ΔD and the generation of hydrocarbons. Specifically, ΔD is strongly negatively correlated with retained hydrocarbons and strongly positively correlated with gaseous hydrocarbon generation. This may be due to retained hydrocarbons blocking some small pores, resulting in a more uniform pore distribution. Conversely, the generation of gaseous hydrocarbons reopens these blocked small pores and creates new small pores, leading to a more uneven pore distribution.

While this experimental study has shortcomings and restrictions in the limited sample size using tested, laboratory systems not able to simulate natural geological conditions, and specific pressures and temperatures chosen not encompassing the full range of conditions experienced by shales in different geological settings, our findings can help optimize exploration strategies and enhance the efficiency of hydrocarbon extraction from shale formations by providing a deep understanding of the thermal maturity process and its impact on hydrocarbon generation and migration. These detailed analyses of pore structure heterogeneity can improve reservoir characterization, leading to an accurate modeling and prediction of shale reservoir performance. Furthermore, utilizing reopened and newly formed small pores in shale can optimize hydraulic fracturing techniques, enhancing hydrocarbon recovery rates.