Characterization of Pore Heterogeneity in Lacustrine Shale Based on MIP, LTNA, NMR, and Multifractal Characteristics: A Case Study of Jurassic Dongyuemiao Member, China

1. Introduction

Shale oil and gas exploration originated in North America. By 2009, shale gas production accounted for over 12% of total natural gas production in North America[1]. In recent years, the development of shale oil, particularly using the Permian Basin marine shale as a reservoir, has also achieved remarkable success in North America. Consequently, the United States has transitioned from an energy importer to an energy exporter, achieving energy independence. China, like North America, possesses substantial shale oil and gas resources, especially in Paleozoic marine shale and Mesozoic-Cenozoic lacustrine shale formations [2,3,4]. After nearly a decade of scientific research, the development of marine shale gas in southern China has reached a significant scale [5,6], notably improving the country’s energy consumption structure. Recently, inspired by the success of marine shale gas development, the vast oil and gas potential in lacustrine shale has been increasingly recognized in the Ordos Basin, Bohai Bay Basin, and Sichuan Basin in China [7,8,9].

Shale, as a unique porous medium, comprises both organic and inorganic components, making it highly complex and irregular [10,11]. The heterogeneity of lacustrine shale is primarily reflected in two aspects: organic geochemical heterogeneity and inorganic mineral heterogeneity. The organic geochemical characteristics of different shales vary significantly, with total organic carbon (TOC) content ranging from 0.5% to 15% and vitrinite reflectance (Ro) ranging from 0.5% to 1.3%. The complexity of mineral composition is evident in the varying dominant mineral types of different lacustrine shales. For example, the Lucaogou Formation in the Jimusarsag composed of felsic minerals, with clay mineral content often below 8%. In contrast, the Ziliujing Formation in the Sichuan Basin features comparable amount of felsic minerals and clay minerals is comparable [12,13,14]. This compositional heterogeneity in lacustrine shale results in a more complex pore structure compared to marine shale [15,16]. The shale pore structure significantly influences the occurrence and permeability mechanisms of hydrocarbons [17,18]. Therefore, effectively characterizing pore structure heterogeneity is a critical factor in reservoir evaluation and enhancing hydrocarbon recovery [19,20].

Various qualitative and quantitative techniques have been used to characterize the pore structure of fine-grained sedimentary rocks [21,22]. Qualitative methods, such as field emission scanning electron microscopy (FE-SEM), transmission electron microscopy (TEM), and computed tomography (CT), are widely applied to study the morphology and connectivity of micro and nanoscale pore networks. Common quantitative methods for determining porosity and permeability include N₂/CO₂ physical adsorption (LTCA and LTNA), mercury intrusion capillary pressure (MIP), small angle/ultra small angle neutron scattering (SANS), and nuclear magnetic resonance (NMR) [23,24,25]. The pore structure of rocks such as coal and granite, in addition to fine-grained sedimentary rocks, can be qualitatively and quantitatively analyzed using the methods mentioned above.

The porosity and permeability of coal can be measured using methods such as mercury intrusion porosimetry, liquid nitrogen adsorption, and nuclear magnetic resonance. According to research, the porosity of coal reservoirs is generally less than 2%, which accounts for a relatively low proportion of the total porosity, while the porosity of bedrock is higher, accounting for about 80%. The permeability of coal samples ranges from 0.0002 × 10^-15 to 0.6525 × 10^-15 square meters, mostly exhibiting medium to low permeability [26]. Han et al. characterized the pore structure of coal with different degrees of metamorphism using mercury intrusion porosimetry and N₂ adsorption methods. They combined the Frenkel-Halsey-Hill and Menger fractal dimension models to accurately define the pore size range of experimental coal. Cai et al. measured the porosity distribution of coal rock using low-temperature nitrogen adsorption/desorption isotherms, evaluated the compressibility of different coal rank coals, and evaluated the pore structure of coal using MIP and N₂ adsorption isotherms at 77 K; Zhao et al. utilized existing pore throat characterization techniques in combination with scanning electron microscopy (SEM), low-temperature nitrogen adsorption (LTNA), high-pressure mercury injection (HPMI), and rate controlled mercury injection (RCMI) techniques to highlight the throat size and pore distribution characteristics of tight sandstone reservoirs; Longinos et al. characterized the pore structure of granite using nitrogen adsorption and nuclear magnetic resonance experiments to study the fracture behavior of heated granite after cooling with liquid nitrogen [27,28,29,30]. Although each method has specific application scenarios and limitations, combining multiple techniques can effectively characterize the full-scale pore system of shale. Fractal theory, as an effective powerful tool for evaluating the complexity and heterogeneity of porous media, is increasingly applied to analyze shale’s complex pore systems. By integrating fractal dimensions with macroscopic and microscopic pore structure characteristics parameters, the correlations between fractal dimensions and these parameters can be thoroughly examined. The Dongyuemiao shale of the Jurassic Ziliujing Formation in the Sichuan Basin is characterized by a high clay mineral content (~ 50%), well-developed shell interlayers, and advanced thermal evolution (~ 1.5%, in the condensate gas generation stage). This terrestrial shale has undergone prolonged diagenesis and hydrocarbon generation, with production wells exhibiting simultaneous oil and gas production at the wellhead [31,32,33]. The lacustrine shale of the Lianggaoshan Formation in the Sichuan Basin has high a TOC content (average 1.57%) and porosity (average 4.21%), along with significant oil and gas content. It is also features well-developed fine sand layers, making it more favorable for hydrocarbon accumulation than the Dongyuemiao section [34,35,36]. Using digital 3D rock models, some researchers have investigated the evolution of fractal dimensions porosity [37]. Wang et al. (2015) evaluated the microporous structure of lacustrine shale using gas adsorption date and fractal theory [38]. These studies demonstrate a strong correlation between fractal dimensions and factors such as mineral composition, specific surface area, and adsorption volume. However, a single fractal model is insufficient for understanding the heterogeneity of pore structures and fails to distinguish samples with identical fractal dimensions but differing pore size distributions [39]. The multifractal approach divides complex fractal structures into multiple regions with uniform singular intensities and generalized fractal dimensions. By analyzing these singular intensities and generalized fractal dimensions region by region, specific multifractal characteristics can be identified. This study focuses on lacustrine shale of the Dongyuemiao Member. Previous research on shale reservoir pore structures has primarily relied on gas adsorption and scanning electron microscopy. As a complementary method, multifractal analysis holds significant potential for characterizing pore structure heterogeneity. This study employs X-ray diffraction, scanning electron microscopy, and nuclear magnetic resonance. Using image processing and multifractal theory, it analyzes the pore structure characteristics and heterogeneity of different lithofacies, providing theoretical support for lacustrine shale reservoir characterization and sweet spot prediction.

2. Geological Setting

The Sichuan Basin is a typical craton basin located in the western part of the Upper Yangtze Block (Figure 1a), covering an area of approximately 260,000 km². During the Early to Middle Jurassic period, most regions of the Sichuan Basin were dominated by shore-shallow lake, semi-deep lake, and deep lake facies, experiencing four lake transgressions during this period (Figure 1a). Organic-rich shale is primarily thought to have developed in semi-deep and deep lake facies under anoxic condition [40,41]. From bottom to top, four sets of organic-rich shales were formed in the Zhenzhuchong Member, Dongyuemiao Member, Da’anzhai Member, and Lianggaoshan 2^nd Member, respectively. The Dongyuemiao Member is further divided into three sub members: the 1^st, 2^nd, and 3^rd sub member. The 1^st sub member is further subdivided into four sections: 1^st, 2^nd, 3^rd, and 4^th Sections (Figure 1b). Except for the 2^nd sub member, the other intervals predominantly constit of shale. The organic matter in the shale is mainly Type II, with localized development of Type III. The vitrinite reflectance (Ro) is mostly greater than1.2%, indicating that shale has reached a high-maturity evolutionary stage and is predominantly generating by gas generation [42].

3. Method and Process

3.1. TOC Content, Mineral Component, and Lithofacies Classification

Samples with varying heterogeneity from layers 1-4 of the Dongyuemiao section were selected for experiments. Total Organic Carbon (TOC) is a key indicator for assessing organic matter abundance in shale. In this study, TOC content was measured using a LECO CS230 carbon-sulfur analyzer. Before analysis, shale samples were ground into a 200-mesh powder. Inorganic carbon was removed by treating the samples with 12.5% dilute hydrochloric acid, followed by rinsing with distilled water to eliminate any residual acid. The sample powder was data temperature of 60-80°C and then placed in the placed in the carbon-sulfur analyzer for complete combustion (at >930 °C). The TOC content was calculated by measuring CO₂ generated during combustion. Mineral composition was analyzed using a Smart Lab X-ray diffractometer. For this analysis, shale samples were ground into 200-mesh powder. Scanning Electron Microscopy (SEM) was used to observe microscopic pores and fractures in the shale. An FEI Quanta 650 FE GFE-SEM was used for qualitatively analysis. Before SEM imaging, shale samples were cut into 2cm×2cm×1cm rectangular blocks, polished with an IM 4000 argon ion polisher, and coated to enhance conductivity, improving the resolution of micropore images.

Lithofacies represent a comprehensive reflection of sedimentary rocks and environments, encapsulating the petrological and sedimentary characteristics of rocks, as well as the macroscopic heterogeneity features of shale, including lithology, mineral composition, and other attributes [43,44,45]. Extensive research has been conducted on shale lithofacies classification and characteristics. For example, Allix et al. (2010) classified shale lithofacies based on the content of clay minerals, carbonate rocks, and felsic minerals [46]. Wang et al. (2012) incorporated organic matter abundance as one of the criteria for lithofacies classification based on mineral composition [47]. Shi et al. (2020) further refined lithofacies classification by considering organic matter abundance, mineralogical facies, and sedimentary structures [48]. A prominent feature of lacustrine shale, particularly in the Dongyuemiao Member, is the presence of high-frequency, variable thickness siliceous and calcareous laminae. In this study, lithofacies classification was conducted by grouping minerals into three categories: calcareous minerals (calcite and dolomite), siliceous minerals (quartz and feldspar), and clay minerals. Using a 50%mineral content threshold, shale was classified as mixed shale (all three mineral contents < 50%), clayey shale (with clay mineral content > 50%), and silty shale. By combining this mineralogical classification with laminae types, lithofacies were effectively categorized.

The main methods for characterizing the micro pore structure of shale include low-temperature nitrogen adsorption, high-pressure mercury intrusion porosimetry, and nuclear magnetic resonance. These methods each have their own advantages and disadvantages, and are suitable for different research needs.

3.2. Characterization Method of Microscopic Pore Structure

3.2.1. Mercury Intrusion Porosimetry (MIP)

The Mercury Intrusion Porosimetry (MIP) method is based on the capillary model, which assumes that porous media consist of capillaries with varying diameters [49]. Mercury, as a non-wetting phase, does not wet the surface of rocks, whereas, air or mercury vapor within the in rocks acts as the wetting phase. During the process, mercury is injected into the rock pores under pressure, displacing the wetting phase. When the injection pressure equals the capillary force, the corresponding capillary radius represents the pore throats. The volume of mercury entering the pores reflects the volume of rock pores connected by the throat. By continuously varying the injection pressure, the pore size distribution curve and capillary pressure curve can be obtained. The calculation formula is as follows:

$$P_c={\displaystyle \frac{2\sigma \cos \theta }{r}}$$

(1)

In the formula, ${\mathrm{P}}_{\mathrm{c}}$ represents capillary pressure (unit: Mpa). $\sigma$ is the interfacial tension between mercury and air, with $\sigma$=480dyn/cm, $\theta$ is the wetting angle between mercury and rock. For $\theta$=140 °, $\cos \theta$=0.765, and r represents the pore radius (unit: μm). The relationship between pore radius (r) and capillary pressure (P_c ) is given by:

$$r={\displaystyle \frac{0.735}{P_c}}$$

(2)

MIP data in this study were obtained using the American Auto Pore IV 9500 fully automatic mercury intrusion porosimeter. Prior to testing, samples were dried to a constant weight at 105 °C. The experiment involved two main processes: pressurized mercury injection and depressurized mercury removal. The maximum experimental pressure reached 200mpa.

3.2.2. Low-Temperature N₂ Adsorption (LTNA)

The low-temperature gas adsorption method is a commonly used technique for characterizing the specific surface area and pore size distribution of porous materials [50]. Nitrogen (N₂) widely used as the detection gas for this purpose. The analysis was conducted using an ASAP 2460 specific surface area and porosity analyzer. Before the analysis, the sample was degassed at 110 °C for 10 hours to remove moisture and volatile substances. After degassing, it was transferred to the analysis station, adsorption experiments were performed using 99.99% pure nitrogen as the adsorbent at 77.3K. The relative pressure (P/P₀) was gradually increased from 0.01 to the saturated vapor pressure (P/P₀ =0.998).

3.2.3. Low-Field Nuclear Magnetic Resonance (NMR)

Saturated water NMR technology relies the nuclear magnetic resonance (NMR) relaxation behavior of spin hydrogen nuclei in rock fluids under a uniformly static magnetic field and radio frequency fields [51]. In this study, samples were saturated with distilled water for 48 hours before measurement. The main frequency 12MHz. NMR relaxation times include longitudinal relaxation time (T₁) and transverse relaxation time (T₂). Under low magnetic field conditions, T₂ measurement is significantly faster than T₁ while providing equivalent pore size distribution analysis. Therefore, T₂ is traditionally used for pore size distribution. The measured T₂ relaxation time comprises three components, expressed by the following formula:

$${\displaystyle \frac{1}{T_2}}={\displaystyle \frac{1}{T_{2b}}}+{\displaystyle \frac{1}{T_{2s}}}+{\displaystyle \frac{1}{T_{2d}}}$$

(3)

where T_2s represents surface relaxation time, T_2d represents diffusion relaxation time, and T_2b represents bulk relaxation time. Diffusion relaxation (T_2d) arises from the self diffusion of fluid molecules, caused by proton spin diffusion through strong internal field gradients. However, in a uniform magnetic field, the magnetic field gradient does not exist, so the diffusion relaxation time (T_2d) is not considered in this experiment. The bulk relaxation time (T_2b) of the fluid depends on its properties, including chemical composition, viscosity, etc. In this experiment, the saturated fluid is water. Since the bulk relaxation time of water is relatively long (usually between 2000-3000ms), and far exceeds the surface relaxation time of pores, 1/T_2b can also be ignored. Thus, the formula (3) simplifies to:

$${\displaystyle \frac{1}{T_2}}={\displaystyle \frac{1}{T_{2s}}}=\rho {\displaystyle \frac{s}{v}}$$

(4)

where $\rho$ represents the surface transverse relaxation strength (unit: μm/ms), which is influenced by the physical properties of the rock. S represents the pore surface area (unit: cm²). The ratio of pore surface area to pore volume (v) is closely related to pore size. Therefore, the formula (4) can be simplified as follows:

$${\displaystyle \frac{1}{T_2}}=\rho {\displaystyle \frac{F_s}{r}}$$

(5)

where F_s is the pore geometry factor. For spherical pores, F_s=3, and r is the pore radius (unit: nm). Based on this formula, it can be inferred that there is a linear relationship between T₂ the spectrum and pore size. Consequently, the formula (5) can be further simplified as:

$$r=C{T}_2$$

(6)

In this formula, C represents a constant, and its value depends on the relaxation properties of the pore surface. The NMR experimental instrument used was a low-field nuclear magnetic resonance device manufactured by Suzhou Niu mag Analytical Instrument Company. For the NMR tests, the shale samples were first processed into cylindrical plunger samples. These samples were then dried at 110 °C for 24 hours. After being evacuated for 12 hours, the shale plunger samples were saturated with aqua pura under a pressure of 25 MPa for 48 hours. After the saturation process, the plunger samples were, removed and left to stand in the saturated fluid for 12 hours before obtaining the NMR T₂ signals. The saturated fluid used was a 20,000 × 10^-6 potassium chloride solution at a saturation pressure of 25 MPa. The NMR test sequence employed was CPMG, with the following parameters: Main frequency: 12 MHz, Echo time (TE): 0.1ms, Number of echoes (NECH): 10,000, Number of scans (NS): 64, Waiting time (TW): 3000ms.

3.3. Multifractal Theory

The box-counting method is a commonly used technique for analyzing and obtaining the multifractal spectra of porous media [52]. It assumes that a set of boxes or sublayers of varying lengths divided the total aperture layer length (L) into equal-scale linear segment (ε), defined as ε=L.2^-k (where k is a positive integer) [53,54,55,56]. In non-uniform porous media, for a given scale ε, the probability of data distribution in the i^th subinterval is defined as $P_i\left(\epsilon \right)$ .

$$P_i\left(\epsilon \right)={\displaystyle \frac{m_i\left(\epsilon \right)}{{\Sigma }_{\mathrm{i}=1}^{N\left(\epsilon \right)}}}={\displaystyle \frac{m_i\left(\epsilon \right)}{M\left(\epsilon \right)}}$$

(7)

$${\Sigma }_{\mathrm{i}=1}^{N\left(\epsilon \right)}{P}_i\left(\epsilon \right)=1$$

(8)

where $\epsilon$ is an arbitrary scale (i.e., box size in box counting) at which the porous media is examined, $m_i\left(\epsilon \right)$ represents the number of pixels or mass in any box, i, at size $\epsilon$, $M\left(\epsilon \right)={\displaystyle \sum }_{i=1}^{N\left(\epsilon \right)}{m}_i\left(\epsilon \right)$ is the total mass or sum of pixels in all boxes for this ε , $N\left(\epsilon \right)$ is the total nonempty boxes for each ε ,P_i (ε) represents the probability of mass$m_i$ at i relative to the total mass for a box.

When the interval becomes sufficiently small, the relationship between P _i (ε) and ε can be expressed as:

$$P_i\left(\epsilon \right)\propto {\epsilon }^{-a_{\mathrm{i}}}$$

(9)

$a_{\mathrm{i}}$ represents the Lipschitz-Hölder index or singularity index, indicating the singularity intensity of the i^th subinterval at scale ε. For multifractals, the target object can be divided into $N\left(\epsilon \right)$ equally long subintervals for a given ε where $N\left(\epsilon \right)$ is expressed as:

$$N\left(\epsilon \right)\propto {\epsilon }^{-f\left(a\right)}$$

(10)

where $a$ measures the local scaling behavior of a measure in the system at scale ε, $f\left(a\right)$ is the fractal spectrum, or singularity spectrum, that characterizes the distribution of singularity strengths $a$ within the system. The shape of $f\left(a\right)$, typically a convex curve, provides insights into the system‘s heterogeneity :a narrow $f\left(a\right)$ indicates with low variability in singularity strengths, whereas a wide f($a$) signifies high variability.

A complex fractal can be divided into subsets represented by different α values. where $f\left(a\right)$ describes the fractal characteristics of these subsets. The fractal spectrum, represented by the $a$~$f\left(a\right)$ relationship, serves as a fundamental descriptor of the local features of multifractals. The formula for defining the moment function is as:

$$M\left(\epsilon \right)={\displaystyle \sum }_{i=1}^{N\left(\epsilon \right)}{m}_i\left(\epsilon \right)$$

(11)

$$f\left(\alpha \right)={\displaystyle \frac{loglog{\Sigma }_{\mathrm{i}=1}^{N\left(\epsilon \right)}loglog{u}_i\left(\epsilon \right)}{loglog\epsilon }}$$

(12)

$$a\left(q\right)={\displaystyle \frac{loglog{{\displaystyle \sum }}_{\dot{i}=1}^{N\left(\epsilon \right)}{u}_i\left(\epsilon \right)log\log p_i\left(\epsilon \right)}{loglog\epsilon }}$$

(13)

Let introduce $u_i\left(q,\epsilon \right)$ denotes the uniform partition function, which is defined as (14):

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

(14)

where q represents the multifractal dimension of an object at different scales, ranging from -10 to 10 in unit steps. The probability distribution function $u\left(q,\epsilon \right)$is then given by:

$$u\left(q,\epsilon \right)={\Sigma }_{\mathrm{i}=1}^{N\left(\epsilon \right)}{p}_i{\left(\epsilon \right)}^q~{\epsilon }^{\tau \left(q\right)}$$

(15)

where${\epsilon }^{\tau \left(q\right)}$ represents the mass scale function, and τ(q) is the scaling exponent or mass exponent. The relationship between τ(q) and q reflects the system’s heterogeneity: it is linear for homogeneous fractals and nonlinear for multifractals. In formula (15), ${\epsilon }^{\tau \left(q\right)}$ is the mass scale function, express as:

$$\tau \left(q\right)=-\underset{\epsilon \to 0}{lim}{\displaystyle \frac{{\Sigma }_{\mathrm{i}=1}^{N\left(\epsilon \right)}{p}_i^q\left(\epsilon \right)}{\log \epsilon }}$$

(16)

Based on formulas (15) and (16), the singularity strength $a\left(q\right)$ and the mass index $\tau \left(q\right)$ can be transformed into the following formulas using the Legendre transformation, yielding the following expression:

$$\alpha \left(q\right)={\displaystyle \frac{\mathrm{d}\tau \left(q\right)}{\mathrm{d}q}}$$

(17)

$$f\left(\alpha \right)=q\alpha \left(q\right)-\tau \left(q\right)$$

(18)

Therefore, the generalized dimension ($D_q$) can be expressed by the following formula:

$$D_q={\displaystyle \frac{1}{q-1}}\underset{\epsilon \to 0}{lim}{\displaystyle \frac{loglog{{\displaystyle \sum }}_{i=1}^{N\left(\epsilon \right)}{p}_i^q\left(\epsilon \right)}{log\log \epsilon }}={\displaystyle \frac{\tau \left(q\right)}{q-1}}$$

(19)

For q>0 and q=1, D₁ can be calculated using the following formula:

$$D_1=\underset{\epsilon \to 0}{lim}{\displaystyle \frac{{\Sigma }_{\mathrm{i}=1}^{N\left(\epsilon \right)}{p}_{i\left(\epsilon \right)}\log p_i\left(\epsilon \right)}{\log \epsilon }}$$

(20)

$\mathsf{\tau }\left(\mathrm{q}\right)$ can also be calculated using the following formula:

$$\tau \left(q\right)=\left(q-1\right)D_q$$

(21)

In the multifractal model of the T₂ spectrum, T₂ is considered the total interval (I) of multifractal. Based on chosen the ε value, the total interval is divided into $N\left(\epsilon \right)=L∕\epsilon =2^k$ subintervals of equal length. To ensure each sub-interval contains at least one value, k is selected as 7 based on the T₂ length, resulting in 128 sub-intervals are divided. Each subinterval sub interval contains two values, yielding 256 points after T₂ interpolation.

4. Results

4.1. Lithofacies Type

The Total Organic Carbon (TOC) content of lacustrine shale in the Dongyuemiao Member ranges from 0.57% to 4.03%, with an average of 1.71% (Table 1). The primary minerals in lacustrine shale are quartz (10.2%-30.6%, average 23.3%) and clay minerals (25.4%-66.8%, average 53.1%), which together constitute 76.4% of the total mineral content. In this study, lacustrine shale samples are classified into four lithofacies types (Figure 2a): shell-laminae mixed shale (SLMS) (Figure 2b), silty-laminae clay shale (SLCS) (Figure 2c), clast-laminae clay shale (CLCS) (Figure 2d), and clay shale (CS) (Figure 2e). Among these, CS has the highest TOC content of CS is the highest, averaging 2.19%, while SLCS has the lowest TOC content, averaging 1.31%. X-ray diffraction (XRD) results reveal variations in mineral composition across the four lithofacies: SLCS has the highest clay content, averaging 64.15% (average 64.15%), while SLMS has the (average 41.5%), CS has the highest quartz content (average 27.2%), whereas SLMS has the lowest, averaging 21.7%.

4.2. Microscopic Pore Types and Morphology

FE-SEM images reveal numerous pores within or between mineral particles, including organic pores, inorganic pores, and micro-fractures (Figure 3). The organic matter in CLCS primarily consists of dense fragments of higher plants, with rare occurrences of solid bitumen. Inorganic pores are highly developed, predominantly comprising inter-particle (Inter-P) pores between clay platelets. These pores are large and exhibit good connectivity. A small number of intra-particle (Intra-P) pores within pyrite framboids are also observed. Micro-fractures are present between clay minerals. In SLMS, the organic matter is compact and organic pore development is minimal. Organic pores are mostly isolated, narrow, and elongated. The development of inorganic pores is significant, including inter-P pores between clay platelets and dissolution pores within carbonate particles. Micro-fractures primarily consist of bedding fractures, with additional micro fractures observed inside shells.

The organic matter content inside CS is high, primarily composed of vitrinite. However, organic pore development is limited, with most organic pores being isolated and circular. Inorganic pores are developed, mainly consisting of Inter-P pores between clay platelets. Bedding fractures dominate the micro-fractures, with a few microcracks observed at the edges of some organic matter. In SLCS, the degree of organic pore development is relatively high. Organic pores are mostly circular or partially. Inorganic pores arehighly developed, mainly comprising inter-P pores within clay minerals with larger pore sizes. Micro-fractures are predominantly bedding fractures.

4.3. Pore Size Distribution

4.3.1. MIP Characterization

The MIP primarily analyzes macropores [56]. The morphology of the MIP curve reflects the pore development and connectivity across different pore throat segments. While the MIP curves vary significantly between lithologies, their overall trend remain consistent (Figure 4). In terms of total mercury input, SLCS exhibits the highest value, indicating the most extensive pore development, followed by SLMS and CLCS, with CS having the lowest input. Initially, the steep curve slope and large mercury injection volume in CS suggest significant macropore development. As pressure increases, mercury progressively infiltrates smaller pores, stabilizing the curve without clear plateaus. When the injection pressure exceeds 10,000 MPa, the SLCS, SLMS, and CLCS curves show a sharp upward trend, while the CS curve remains relatively flat. The final mercury injection for CS is generally smaller than for other lithofacies. In the experiment, mercury saturation at the maximum pressure indicates the point at which mercury can no longer penetrate pores. During pressure release, mercury begins to escape, with the remaining mercury saturation representing the amount trapped in the sample. Mercury removal efficiency is calculated using the following formula:

$$w=\left(1-{\displaystyle \frac{s_r}{s_{max}}}\right)\times 100\%$$

where $w$ represents the mercury removal efficiency (unit: %). S_max is the mercury saturation at the maximum mercury injection rate (unit: %), and S_r is the residual mercury saturation (unit: %). Based on the mercury intrusion-mercury removal curves, the mercury removal efficiency is determined, which provides insights into pore connectivity. CLCS has the highest mercury removal efficiency, at 25.51%, suggesting it has the best pore connectivity among the four lithofacies. SLMS follows, with a mercury removal efficiency of 22.78%. SLCS and CS have the lowest mercury removal efficiencies, at 21.47% and 21.48%, respectively, with no significantly difference between them.

4.3.2. LTNA Characterization

LTNA is mainly used for analyzing micropores and mesopores. The LTNA curve shape reflects the partial pore size distribution within shale samples (Figure 5a). Curve for different lithofacies all exhibit an inverse "S" shape. The adsorption curve shows a slow upward trend in the low-pressure area, with a clear inflection point between relative pressures of 0.4-0.6, followed by a rapid increase in the high-pressure area. Due to variations in pore structure and morphology, the slope of the adsorption curve increases and the saturation adsorption plateau is not well-defined at high relative pressure (Figure 5b) This indicates that the desorption and adsorption curves do not overlap, resulting in a significant hysteresis loop (Figure 5c). According to the Kelvin equation, a hysteresis loop forms when capillary condensation occurs and the relative pressure (P/P₀) exceeds 0.5. Based on the IUPAC classification method, the desorption curves for all four lithofacies display distinct inflection points within the relative pressure range of 0.4 to 0.6. These characteristics correspond to the H3 hysteresis loop (Figure 5d), suggesting that the pores in all four lithofacies are predominantly narrow slit-shaped and wedge-shaped pores.

4.3.3. NMR Characterization

In this study, four types of T₂ curves were identified, each directly reflecting different pore size distributions [58]. The T₂ spectrum curve of CS exhibits a three-peak shape (Figure 6a), The first peak, on the far left, represents the pore size distribution of micropores, corresponding to relaxation times between 0.03 and 0.1ms. The second peak, with relaxation times between 0.3~5ms, represents micropores and mesopores. The third peak, with relaxation times ranging from 15 to 60ms, corresponding to the pore size distribution of macropores. Among the three peaks, the second peak is the most prominent, indicating a higher abundance of micropores and mesopores in CS. The T₂ curve of CLCS also displays a three-peak shape (Figure 6b). Compared to CS, CLCS shows better micropore development but less macropore development. The T₂ curve of SLMS exhibits a bimodal shape (Figure 6c). The first peak, with relaxation times ranging from 0.03 to 8ms, reflects the pore sizes of micropores to mesopores. The second peak, with relaxation times between 50-150ms, represents the pore size distribution of macropores. Compared CS and CLCS, SLMS shows reduced development of micropores, mesopores, and macropores. The T₂ curve of SLCS also exhibits a bimodal shape (Figure 6d). However, compared to SLMS, the development of micropores, mesopores, and macropores is further diminished. Compared to the "three-peak" pore size distribution, the bimodal micropore structure shows significantly reduced complexity.

4.3.4. Full Pore Size Characterization of Shale

LTNA is advantageous for characterizing pore size distributions larger than 1nm, while MIP is better suited for analyzing mesoporous and macroporous structures. According to previous research, pores smaller than 2 nm are referred to as micropores, 2-50 nm as mesopores, and those larger than 50 nm as macropores [50]. Therefore, in this study, LTNA results were used to characterize micropores larger than 1 nm, and MIP data were employed for larger pore distributions that LTNA cannot detect. The combination of these two methods enables comprehensive full pore size characterization of shale (Figure 7a, 7b, 7c, 7d). As shown in Table 2, the statistical results indicate that the pore volume of the four lithofacies is predominantly composed of mesopores and macropores, with similar contributions from both.

4.4. Fractal Dimension Based on NMR Results

Wang et al. analyzed the pore structure of shale nanopores using fractal theory and multifractal theory, and calculated the fractal dimension and multifractal parameter capillary pressure under different fractal models. The results showed that shale nanopores have multifractal characteristics, and multifractal parameters can reflect the asymmetry of pore size, concentration, and pore size distribution. Shales with similar multifractal parameters have similar pore size distribution (Δ$a$, Δ f, $a$₀, $a$₁, $a$₂). The information dimension D₁ and correlation dimension D₂ are positively correlated with the size of shale nanopores. As the information dimension D₁ decreases, the contribution of pore size distribution becomes greater, and there is a strong negative correlation between the information dimension D₁ and the fractal dimension of the three-dimensional capillary model [5960]. Based on the T₂ spectral distribution and the multifractal model, the NMR multifractal characteristics under saturated water conditions were analyzed. The multifractal parameters for each lithofacies are presented in Table 3. The default setting for the multifractal moment order is a continuous interval between -10 and 10, assuming the statistical interval order is q=1. Figure 8 shows the generalized fractal dimension spectra ($D_q$~q) for different lithofacies, displaying an overall inverse S-shaped curve. The rate at which $D_q$ decreases with order (q) reflects the heterogeneity of the system. A faster decrease indicates stronger heterogeneity. Among the four lithofacies, CLCS exhibits the highest multifractal dimension difference (ΔD), signifying the strongest heterogeneity. As q increases, $D_q$ demonstrates two distinct trends:

When q>0, the multifractal parameters decrease significantly, highlighting the pore structure in high-probability regions.

When q<0, the multifractal parameters decrease more gradually, reflecting the pore structure in low probability regions.

Using q=0 as the boundary, (D_min-D₀) represents the heterogeneity of macropores. (D₀-D_min) indicates the heterogeneity of micropores and mesopores.

For multifractal spectra, the Δ$a$ value serves as an indicator of heterogeneity; larger Δ$a$ values signify more complex pore distributions and greater heterogeneity. Figure 9 shows the fractal dimension spectra for different lithofacies, with the peak values ($a$₀, f($a$₀ )) serving as dividing points: ($a$₀-$a$_min) is more suitable for evaluating the pore structure in low probability regions, such as macropores. Conversely, ($a$_max-$a$₀) better describes the pore structure in high probability regions, such as micropores. The fractal spectrum distribution f($a$) vs.$a$ in Figure 9 is convex (as expected), indicating the multifractal nature of porosity, and is generally asymmetric. On the left side of ($a$₀, f($a$₀)), the fractal spectrum distribution f($a$) increases sharply with increasing of $a$. On the right side, f($a$) gradually decreases as a increases. Among the four rock phases, CLCS exhibits the largest decrease in the generalized fractal dimension (Δ$a$), indicating a highly complex pore size distribution and significant heterogeneity within the pore space.

The f($a$) distribution in the CS and CLCS lithofacies is right-skewed, suggesting that large-scale, high-probability features (e.g., macropores) contribute more significantly to the system's heterogeneity. In contrast, the SLMS and SLCS lithofacies exhibits a more symmetrical f($a$), indicating no predominant micro-or macro-porosity but rather a more uniform distribution between them (Figure 8). The width of the f($a$) spectrum reflects the range of singularity strengths ($a$), which are associated with local scaling behaviors of the pore distribution: A narrow f($a$) implies uniform porosity with low heterogeneity, as observed in the CS-Well B (2509.7m) and CLCS-Well A (2937.66m) samples, while a wide f($a$) signifies a highly heterogeneous system with diverse pore sizes and distributions, as seen in the CS-Well B (2493.09m), CLCS-Well B (2507.62m), and Well A (2940.08m) samples.

5. Discussion

5.1. Controlling Factors of Heterogeneity in Pore Structure

Multifractal theory reveals an intrinsic relationship between the heterogeneity of pore structure and the varying effects of, different components and organic matter content in rocks on pore architecture. Li et al. (2024) believe that the more brittle minerals and the less clay minerals in the mineral composition, the more complex the pore structure of shale and the stronger its heterogeneity [60]. The results show that the heterogeneity parameters of CS are significantly negatively correlated with quartz content (R²=0.93) (Figure 10a). This suggests that high quartz content promotes the preservation of pores, resulting in SLCS have a well-developed pore structure, lower heterogeneity, and consequently, a smaller fractal dimension. In contrast, the heterogeneity parameters of CLCS and SLMS exhibit a positive correlation with clay content (R²=0.57 and 0.74) (Figure 10b, 10c). Clay minerals are susceptible to alteration during diagenesis, and compaction often destroys the pore spaces within and between clay minerals. These processes lead to a more complex pore structure, increased heterogeneity, and higher fractal dimensions. The heterogeneity parameters of SLCS are significantly positively correlated with TOC (R²=0.54) (Figure 10d). Previous studies [60] suggest that higher TOC content leads to more developed internal micro- and mesopores, the higher the TOC content, resulting in a more uniform pore size distribution and reduced heterogeneity. However, this study argues that increased organic matter abundance shale permeability to some extent. However, the pore size distribution of micropores and mesopores developed within organic matter is more complex, contributing to increase heterogeneity, as reflected by the higher in fractal dimensions. In the future, principal component analysis（PAC）method can be used to further analyze the main controlling factors of heterogeneity in shale reservoir pore structure.

5.2. Multifractal Parameters Reflecting Heterogeneity of Porosity and Permeability

Physical parameters comprehensively reflect the shale pore system, with pore gap and permeability being the most critical factors influencing storage capacity and flow. These properties determine reservoir quality [61]. Establishing relationships between physical properties and multifractal parameters helps in understanding the heterogenous nature of shale reservoirs. Previous studies found that the multifractal dimension of shale is weakly correlated with porosity and permeability. The correlations between the porosity of CS, CLCS, SLCS and (D_min-D₀) are significant (R²=0.92, 0.55, and 1) (Figure 11a, 11b, 11d). (D_min-D₀) is effective in characterizing the heterogeneity of porosity in these three lithofacies (Figure 11c). A larger (D_min-D₀) value indicates porosity heterogeneity. In contrast, (D₀-D_max) is more suitable for characterizing the porosity heterogeneity in SLMS, where a high value represents strong heterogeneity, and a low value indicates the opposite. The permeability of the four lithofacies is most correlated with (D₀-D_max) and ΔD (R²=0.65, 0.60, 0.55 and 0.74). The permeability of CS, SLMS, and SLCS is positively correlated with both (D₀-D_max) (R²=0.62, 0.34, and 0.88) (Figures 12a, 12c, 12d), while higher values indicate stronger permeability heterogeneity. Conversely, the permeability of CLCS is negatively correlated with (D₀-D_max) (Figure 12b), with lower values indicating stronger permeability heterogeneity.

6. Conclusions

(1) Based on mineral content and laminae types, the Dongyuemiao lacustrine shale lithofacies can be classified into four types: shell-laminae mixed shale (SLMS), silty-laminae clay shale (SLCS), clast-laminae clay shale (CLCS), and clay shale (CS).

(2) The micro-reservoir and permeability spaces of lacustrine shale consist of inorganic pores, organic pores, and micro-fractures. Inorganic pores include inter-particle and intra-particle pores. Full pore-size distribution characterization reveals two main peaks in the pore-size distribution curves of all the four lithofacies, with pore sizes concentrated in the 2-10 nm and 50-80 nm ranges. Mesopores and macropores constitute the majority of the total pore volume in lacustrine shale.

(3) The heterogeneity of CS is primarily controlled by quartz content, with higher quartz content leading to weaker heterogeneity. In SLMS and CLCS, heterogeneity is controlled by clay mineral content, which greater clay mineral, development results in stronger heterogeneity. For SLCS, heterogeneity is mainly controlled by TOC content, with higher TOC levels increasing to increased heterogeneity. The heterogeneity of porosity in SLMS is best characterized by (D₀-D_max), while in the other three lithofacies, it is characterized by (D_min-D₀). Permeability of all four lithofacies is most strongly correlated with (D₀-D_max) and ΔD. For CS, SLMS, and SLCS, permeability is positively correlated with (D₀-D_max), where high values indicate strong permeability heterogeneity. Conversely, in CLCS, permeability is negatively correlated with (D₀-D_max), where low values represent strong permeability heterogeneity.