Characteristics and Inﬂuencing Factors of Multi-Scale Pore Structure Heterogeneity of Lacustrine Shale in the Gaoyou Sag, Eastern China

: The success of shale oil exploration and production is highly dependent on the heterogeneous nature of the reservoir pore structure. Despite this, there remains limited research on the heterogeneity characteristics of pores at different scales in lacustrine shale oil reservoirs and the factors that impact them. This study aims to quantitatively characterize the multi-scale pore heterogeneity differences of the lacustrine shale


Introduction
In recent years, China has made significant progress in shale oil production with the increase in shale oil exploration. In 2021, China's shale oil output amounted to 240 × 10 4 t, up by 60% from 2019. In addition, the recoverable resources of shale oil across the country were 7.4-37.2 billion tons, as estimated by the National Energy Administration

Geological Setting and Samples
The Subei Basin is a Mesozoic and Cenozoic sedimentary basin founded on the Cretaceous basement. [31,32]. The Subei Basin is located in eastern China and encompasses approximately 3.5 × 10 4 km 2 . With a size of roughly 3000 km 2 , the Gaoyou Sag is situated in the southern part of the Subei Basin ( Figure 1). In the Gaoyou Sag, there is a sequential deposition of five formations, namely, the Cretaceous Taizhou Formation, the Paleogene Funing Formation, the Dainan Formation, the Sanduo Formation, and the Neogene Yancheng Formation, from the bottom up ( Figure 2). The Funing Formation is the primary shale development horizon in the Gaoyou Sag. Gray and dark gray shale combined with siltstone and fine sandstone are produced in the Funing Formation, which was formed in a deep lake to semi-deep lake sedimentary environment [32,33]. The Funing Formation is typically more than 300 m thick, and because of its excellent organic geochemical conditions, it is regarded as a good source rock. With significant shale oil resource potential, several wells have been used to drill shale reservoirs in the Funing Formation [32,34]. A total of 17 samples of shale from the Funing Formation in the Gaoyou Sag were chosen for this investigation, and all samples were from Well F2 at a depth of 3460~3688 m.

TOC, Mineralogy, and Lithofacies Division
The content of the total organic carbon (TOC) is one of the key matrices used to assess the abundance of organic matter. Using a CS230HC (CS230HC, LECO, San Jose, CA, USA) carbon and sulfur analyzer, the TOC contents of core samples from the Funing Formation were determined in this investigation. The samples were crushed to powder samples of 0.1~0.5 g. (particle size of 200-mesh). Diluted hydrochloric acid was employed to eliminate the samples' inorganic carbon prior to the experiment. Following the procedure, the samples were dried at 60 to 80 • C. Finally, the samples were burned at a high temperature (930 • C), the CO 2 contents were determined, and the samples' TOC contents were computed in accordance with the results [35]. An MPV-I microphotometer was used to measure the vitrinite reflectance value (R o ) [36].
The mineral components of the samples were quantitatively examined using a Smart-Lab X-ray diffractometer (SmartLab, Rigaku Corporation, Tokyo, Japan). To prepare the samples for analysis, they were crushed into a 200-mesh powder. By examining the distinctive diffraction peak intensities of the X-ray diffraction patterns of individual mineral crystals, it was possible to identify the different minerals present in the samples [37,38].
Scanning electron microscopy can be used to observe the microscopic pores and fractures of shale [39][40][41]. In this study, the pore properties of the shale samples were qualitatively observed using a high-resolution field emission scanning electron microscope, the ZEISS Gemini SEM 500 (Gemini SEM 500, ZEISS, Oberkochen, Germany), with a resolution of up to 0.8 nm. The samples were cut into small pieces measuring 1.5 cm × 1 cm × 0.5 cm before the experiment and then polished using a HITACHI IM 4000 (IM 4000, HITACHI, Tokyo, Japan) argon ion polishing apparatus. Field-emission scanning electron microscopy observations of the results were made.
The lithofacies division method in this study was based on the three-terminal element method [42,43]. The specific method involved dividing the minerals in the shale into three categories, namely calcareous minerals (calcite and dolomite), siliceous minerals (quartz and feldspar), and clay minerals. With 50% as the dividing line of mineral content, it can be divided into four lithofacies types: mixed shale (the content of three types of minerals is less than 50%), clayey shale (the content of clay minerals is more than 50%), calcareous shale (the content of calcareous minerals is more than 50%), and siliceous shale (siliceous mineral content is more than 50%).

N 2 Adsorption and MIP
The low-temperature gas adsorption method is a common characterization method for measuring the specific surface area and pore size distribution of porous materials. N 2 is one of the most widely used gas detection media [44]. The ASAP 2460 specific surface area and pore size analyzer was used to test 60-mesh samples. Before starting the analysis, the samples were put into a degassing station and degassed at 110 • C for 10 h to remove the moisture and volatile substances. Then, the degassed samples were transferred to the analysis station, and an adsorption test was carried out at 77.3 K with high-purity N 2 as the adsorbate. The relative pressure (P/P 0 ) gradually increased from 0.001 to 0.998 (saturated vapor pressure of N 2 ).
The AutoPore IV 9505 (AUTOPORE IV 9505, Micromeritics, Norcross, USA) automatic mercury porosimeter was employed in this study. The mercury injection pressure provided by the instrument in the test stage reached 200 MPa, and the corresponding minimum test pore size was approximately 7 nm. The samples were formed into cylinders with a diameter of 2.5 cm and a length of 5 cm before the test. The samples were also continually dried for 48 h at 60 • C in order to successfully remove any internal contaminants. Pore structure factors, such as the pore volume and pore size distribution, were determined using the Washburn equation [45,46]: where D is the pore diameter, µm; P is the pressure of the mercury injection, MPa; γ is the surface tension of the mercury, generally 0.48 N/m; θ is the contact angle between the mercury and the solid dielectric material (θ = 140 • for a non-wetting phase with [47,48]).

Fractal Dimension
The fractal theory was first proposed by the French mathematician Mandelbert in 1973 [49]. Initially, it was utilized to quantitatively characterize the pore structure of porous media such as coal and shale reservoirs. Fractal theory is a method for studying complex, fragmented systems in nature that are both self-similar and self-inverting [50]. Fractal dimension values can be used to quantify a system's fractal properties. Numerous data and models, including the mercury intrusion method, the low-temperature gas adsorption method, the high-pressure mercury intrusion experiment, and the scanning electron microscope method, can be used to compute the fractal dimension [21,51,52].
In order to calculate the fractal dimension of the pore structure in porous solid media, the fractal Frenkel-Halsey-Hill (FHH) model, which is based on a low-temperature N 2 adsorption experiment, is frequently utilized. PFEIFER proposed the FHH model, which is a straightforward model with broad applicability [53]. The FHH model has been applied to various rock types, including sandstones, carbonates, and shales, and has been found to provide a good approximation of the porous structure of these rocks [54][55][56]. Its formula is: where V is the gas adsorption amount corresponding to the balanced pressure P, cm 3 /g; P is the balance pressure, MPa; P 0 is the saturated steam pressure of the adsorbed gas, MPa; C is a constant; D is the fractal dimension of the porous material, and the value of D is between 2 and 3 [57]. The closer the D value is to 3, the more complex the pore structure in the porous media is and the rougher the surface is. The closer the D value is to 2, the simpler the pore structure is. For the determination of the fractal dimension of the macroporous pore structure in porous solid media, the capillary model based on the high-pressure mercury intrusion experiment is widely used [58]. Based on the geometric principles and basic elements of the capillary model, which include the relationship between the pore size and the capillary pressure, as well as the relationship between the pore size and the mercury saturation, the subsequent model has been further developed: where S Hg is the mercury saturation, %; LgP c is the capillary pressure, MPa; LgP s is the inlet capillary pressure, MPa.

Types and Morphology of Pores
The majority of the Funing Formation shale's pores are made of inorganic minerals ( Figure 3). The two primary types of inorganic mineral pores are clay mineral pores ( Figure 3d) and intragranular pores ( Figure 3a). With an average pore size of tens to hundreds of nanometers, the pores of the clay minerals are mainly elongated and triangular. Intragranular pores are formed mainly by the dissolution of soluble minerals such as feldspar and calcite, generally circular or oval in shape. Intragranular pores are typically less than 100 nm in size, isolated, and densely dispersed (Figure 3a). Intergranular pores ( Figure 3b) are primarily formed between brittle mineral particles such as pyrite and quartz. Typically, these pores have sizes in the hundreds of nanometers, although they may occasionally exceed micrometers. The Funing Formation shale exhibited fewer organic pores than marine shale, and the majority of the organic material observed in the electron microscope images was non-porous ( Figure 3c).

Pore Structure from N 2 Adsorption and MIP
The shale samples from the various Funing Formation lithofacies were all reverse-Sshaped ( Figure 4). There was a clear lagged loop between the adsorption and desorption isotherms. Lagged loops are produced when capillary condensation takes place and the relative pressure (P/P 0 ) surpasses 0.5, in accordance with the Kelvin equation [59,60].
The pore shape of the materials in porous media can be seen in the N 2 adsorption/desorption isotherm. The lagged loops of all samples were close to the H 2 type, which corresponds to the ink-bottle pore shape, according to the IUPAC categorization method [59]. The H 3 type, which corresponds to parallel plate-like pores, is also present in the lagged loops of MS and SS. The findings demonstrate that parallel plate pores as well as ink-bottle pores were found in the shale samples of these two lithofacies.
The features of shale's pore structure can be determined using MIP. Figure 5 depicts the mercury injection-ejection curves for the samples. Different lithofacies' mercury injectionejection curves had similar growth trend features, and they often began to rise quickly at pressures of around 30 MPa. However, the mercury withdrawal efficiency (MWE) of various lithofacies varied ( Figure 5 and Table 2). The average mercury withdrawal efficiency for CS was 25.76%, which was considerably lower than that for MS (45.94%) and SS (47.91%). This indicates that the connectivity of CS is inadequate. Calcite and dolomite, two soluble carbonate minerals found in high concentrations in CS, form discrete and shoddily connected ink-bottle corrosion pores (Figures 3a and 4c).   To acquire the pore development characteristics of shale at the micro and nano sizes, the method of combining N 2 adsorption and MIP was chosen due to the complexity of the pore structure of the Funing Formation shale. The pore size distribution obtained by the N 2 adsorption experiments was calculated based on the capillary condensation mechanism. According to the Kelvin equation, the larger the pore, the greater the relative pressure required for capillary condensation. The relative pressure at the time of capillary condensation in the macropores was close to that of the saturated vapor pressure, which was challenging to detect in the experiment. Therefore, to characterize the pore volume and specific surface area of pores with a pore size of less than 50 nm, the density functional theory (DFT) and BET models of the N 2 adsorption method were utilized. The information about the pore size was obtained in accordance with the quantity of mercury molecules entering the shale under high pressure. The size of pores that mercury molecules can permeate through reduces with an increase in pressure. However, excessive pressure can lead to generated fractures and pore compression, which can produce false experimental results. MIP is therefore typically employed to examine pores with a diameter range greater than 50 nm. Pore size distributions of <50 nm and >50 nm were obtained by N 2 adsorption and MIP, respectively, depending on the range of applications of the two techniques. A thorough analysis of the pore structure parameters (PV, SSA, and MWE) of various scales was performed ( Figure 6 and Table 2). Among the three lithofacies, the PV and SSA of the CS were the largest (Table 2), with ranges of 0.0093~0.0104 cm 3 /g and 2.2467~2.8540 m 2 /g and averages of 0.0099 cm 3 /g and 2.5504 m 2 /g, respectively. The PV and SSA of the MS were relatively small, ranging from 0.0053 to 0.0130 cm 3 /g and from 0.3833 to 3.6429 m 2 /g, respectively, with averages of 0.0087 cm 3 /g and 1.8264 m 2 /g. The PV and SSA of the SS were medium, with ranges of 0.0086~0.0112 cm 3 /g and 0.8741~3.1888 m 2 /g and averages of 0.0097 cm 3 /g and 2.0012 m 2 /g, respectively. In general, the average PV of the CS in the study sample was 1.13 and 1.02 times that of the MS and SS, while the average SSA of the CS was roughly 1.40 and 1.27 times that of the MS and SS. Figure 4 shows that the pore size distribution of the shale from the Funing Formation mainly contained three peaks, at <10 nm, 20~30 nm, and >10,000 nm. Figure 1 shows that the last peak was mainly caused by microcracks, whereas the preceding two peaks were mainly caused by corrosion pores. Small pores with a diameter of <50 nm, which make up the majority of the Funing Formation shale reservoir space, produced 89% of the PV and 99% of the SSA.

Fractal Dimension 4.4.1. Fractal Dimension Based on N 2 Adsorption
The fractal dimensions D 1 and D 2 were calculated using the nitrogen adsorption data and the FHH model (Equation (2)). Based on the N 2 adsorption data with relative pressure between 0 and 0.5, D 1 was calculated, with values ranging from 2.4948 to 2.5681 (Table 3 and Figure 7). D 1 refers to the fractal dimension of the pore surface and is mainly governed by van der Waals force, representing the degree of irregularity of the shale surface [61]. Based on the N 2 adsorption data with relative pressure between 0.5 and 1, D 2 was calculated, with values ranging from 2.5952 to 2.8513 (Table 3 and Figure 7). D 2 , as a result of capillary condensation, represents the degree of irregularity of the pore structure, which is called the fractal dimension of the pore structure [61]. The different lithofacies of the shale demonstrates that D 2 was greater than D 1 , showing that the pore structure is more intricate than the pore surface ( Figure 7).

Fractal Dimension Based on MIP
Shale from the Funing Formation features a complicated pore structure with multifractal characteristics (Figure 8 and Table 3). The two-segment fractal characteristics of the MIP fractal curve are evident ( Figure 8). As a result, the pore complexity values of various pore ranges matched to the fractal dimension were calculated by varying subsection pressure. In this study, D c1 and D c2 represent the optimal fractal properties of pores with widths between 50~100 nm and >100 nm, respectively. D c1 ranged from 2.5031 to 2.9158 (average 2.7096), and D c2 ranged from 2.9612 to 2.9985 (average 2.9937). In general, D c1 < D c2 , suggesting that, as pore size grows, the complexity of the pore structure also increases. Based on the idea of comprehensive fractal dimensions, a summary of the multisegment fractal dimensions calculated using the MIP data is provided, allowing for easier comparison of the differences among the fractal dimensions of various sizes. According to the amount of PV within each diameter range, the scattered fractal dimension was then combined to create a total fractal dimension D 3 , as given in Formula (4) [62]. D 3 is used to characterize the complexity of reservoir spaces with diameters greater than 50 nm. The range of D 3 was 2.7634 to 2.9858, with an average value of 2.9194 and a median value of 2.9390, indicating that the macropores are highly heterogeneous, and the differences among the lithofacies are very small.
where D t is the comprehensive fractal dimension; D i is the fractal dimension of different pore diameters; V i is the proportion of pore volume in different pore sizes, %. Figure 7 shows the correlation between the fractal dimension and the pore structure parameters of the Funing Formation shale samples. As shown in Figure 9a,b, D 1 and D 2 were positively correlated with the PV and SSA. According to Table 2, 89% of the PV and 99% of the SSA were contributed to by pores under 50 nm in size. From Figure 10a,b, it can be seen that the PV and SSA had a significant negative correlation with the AVE PD. Therefore, the larger the SSA and PV, the larger the number of small pores, and the larger the D 1 and D 2 . The regression equation between the roughness of the shale rock pore surface (D 1 ) and the AVE PD had a slope of −0.0058 and an R 2 value of 0.361, indicating a limited impact of the AVE PD on D 1 . The roughness of the pore surface in shale is primarily influenced by factors such as the mineral composition and depositional environment, while factors such as the pore size, pore volume, and connectivity have a minimal impact [26,29,55,63]. There is little variation in the degree of roughness among pore surfaces with different pore sizes ( Table 3). The mineral composition and depositional environment of the 17 samples in this study were highly comparable, as demonstrated in Table 1. As a result, there was minimal variation in the D 1 values among the different samples, as shown in Table 3. The structural complexity of the small pores, as measured by D 2 , was found to be significantly influenced by the AVE PD, as indicated by the slope of −0.031 in the regression equation between the two parameters, with an R 2 value of 0.9467. This suggests that there is a substantial difference in the structure complexity of small pores with varying aperture sizes. The main source of the pore volume was the small pores, and the smaller the AVE PD (Figure 10b), the greater the number of small pores, and the more complex the structure of these small pores tended to be. As depicted in Figures 9 and 10, it can be observed that the differences in the D 1 and D 3 values among the different samples were relatively small, while the differences in the pore structure parameters were often attributed to differences in the development of small pores. As shown in Table 2, large pores contribute less to the pore volume and specific surface area, a result of which D 3 is irrelevant to the PV, SSA, and AVE PD. Furthermore, the correlation coefficient R 2 was smaller than 0.1. As for D 1 , D 2 , and D 3 , they were negatively correlated with the MWE. For them, the correlation coefficient R 2 values were 0.168, 0.233, and 0.3119, respectively (Figure 7d). The greater the pore network connectivity of shale, the higher the regularity of pores, and the lower the D 1 and D 2 . In addition, D 3 had the highest correlation coefficient with regard to the MWE, and the MWE shows a positive correlation with the AVE PD, as can be seen in Figure 10c. This suggests that regular macropore development has a considerable impact on the connectivity of the pore network. As shown in Figure 3a, nonporous organic matter accounts for a large majority of the organic matter in the shale samples of the Funing Formation. According to a previous study [64], the lower limit of maturity of organic pores in shale is 0.9%. The vitrinite reflectance of the sample ranged mainly between 0.7 and 0.9 (Table 1), and the maturity of organic matter remained low. This hindered the organic matter from generating a large amount of gas to produce pores. Instead, macropores and microfractures developed on a small scale (Figure 3c). On the one hand, with the increase in TOC, the space of development was reduced for small pores, which caused significant reductions in the D 1 and D 2 (Figure 11a). Due to the decline in the number of small pores, both the PV and SSA decreased (Figure 11b,c). On the other hand, the connectivity of the pore network was improved due to the small number of large pores and micro-fractures produced by organic matter. Therefore, the relationship between the MWE and the TOC content of organic matter was insignificantly positive (Figure 11d). According to the scanning electron microscope observations (Figure 3), the large pores mainly consisted of a large quantity of clay mineral pores and intergranular pores and a small number of organic matter pores, with the latter contributing little to the pore structure of the large pores. Therefore, the level of the structural complexity of the large pores (D 3 ) does not depend on the amount of organic matter. Consequently, there was no clear correlation between the macropore fractal dimension D 3 and the TOC (Figure 11a). In summary, the D 1 , D 2 , PV, and SSA decreased with an increase in the TOC, but the connectivity was improved slightly. Figure 11. Impact of TOC on fractal dimension (a)and pore structure parameters (b-d).

Impact of Lithofacies on the Heterogeneity of Pore Structure
Compared to the MS and SS lithofacies, the PV and SSA of the CS lithofacies were higher (Figure 12a), reaching 0.0093~0.0104 cm 3 /g and 2.2467~2.8540 m 2 /g, respectively ( Table 2). In addition, the MWE (25.76%) of the CS lithofacies was the lowest, reaching a significantly lower level than that of the MS (45.94%) and CS (47.91%), which indicates it has the worst connectivity. Since the content of soluble carbonate minerals (calcite and dolomite) in CS exceeds 50%, there were a large number of isolated and poorly connected dissolution small pores in development (Figure 3a). Therefore, the pore structure of the CS shale pores was highly heterogeneous, and the shale surface was rough. In addition, the D 1 and D 2 reached a high level, the average values of which were 2.5439 and 2.8178, respectively (Table 2 and Figure 13).  The PV and SSA of the SS were moderate, reaching 0.0086~0.0112 cm 3 /g and 0.8741~3.1888 m 2 /g, respectively. Similar to calcite and dolomite, SS contains a small proportion of feldspar (16.02%~22.69%), making it easy for isolated and poorly connected dissolution small pores to develop. In addition, the contents of rigid minerals, such as quartz and feldspar, in the SS shale reached above 50%. Since the rigid framework composed of quartz, feldspar, and other minerals promotes the preservation of primary intergranular pores, the SS shale contains a small number of small pores and well-preserved macropores. For this reason, the AVE PD of the SS shale was high, the connectivity was satisfactory (Figure 12), and the heterogeneity was moderate ( Figure 13).
The average PV and SSA of the MS were the smallest, reaching 0.0087 cm 3 /g and 1.8264 m 2 /g, respectively. By contrast, the AVE PD and MWE were higher ( Figure 12). MS contains plenty of clay minerals (>25%), which tend to aggregate with organic matter and pyrite to develop macropores and microfractures (Figure 3b,d). Therefore, the AVE PD of the MS was high and the pore connectivity was excellent (Figure 12b). Compared to siliceous minerals and carbonate minerals, clay minerals have greater plasticity but more susceptibility to damage caused by compaction, thus resulting in plate pore morphology ( Figure 3d). Therefore, the D 1 and D 2 of the MS were low, the averages of which were 2.5291 and 2.7450, respectively (Table 2 and Figure 13).
To sum up, the D 2 varied significantly among the different lithofacies, but the D 1 and D 3 barely differed. The shales of the different lithofacies show the order of D 3 > D 2 > D 1 . The pore structure is more complex than the pore surface, and the heterogeneity of large pores is greater compared to small pores. Among the three lithofacies, the CS had the largest PV, SSA, D 1 , and D 2 , indicating the development of a more complex pore structure network. This expands the space needed for shale oil occurrence. However, the connectivity of the CS lithofacies was the lowest among the three, which hinders shale oil production. Although the PV of the SS was slightly lower than that of the CS, its average pore diameter (AVE PD) and connectivity are significantly advantageous, making the SS an ideal shale reservoir.

Conclusions
(1) The Funing Formation shale is in the mature stage, with a TOC of about 60% in the samples exceeding 1%, and the source rock conditions are advantageous. The primary minerals of the Funing Formation include quartz and clay, followed by calcite and dolomite, but there are no obviously dominant minerals, with the average content falling below 30%. The lithofacies are comprised mainly of mixed shale (58.82%) and siliceous shale (29.41%); (2) The microscopic pores of the Funing Formation shale are dominated by inorganic mineral pores and fewer organic pores. Intragranular pores and clay mineral pores are two types of inorganic mineral pores that have developed widely. The MWE of the Funing Formation shale is usually lower than 50%. Among the three lithofacies, the CS had the worst MWE (25.76%), indicating its poor connectivity. Ink-bottle pores and plate pores represent the main pore morphology. The pore size distribution of the Funing Formation shale had three major peaks, at < 10 nm, 20~30 nm, and >10,000 nm respectively. In addition, 89% of the PV and 99% of the SSA are attributed to the pore sizes below 50 nm, both of which are the principal components of the Funing Formation shale reservoir space; (3) The fractal dimensions D 1 , D 2 , and D 3 were calculated respectively to characterize the roughness of the pore surface, the structural complexity of the small pores, and the structural complexity of the large pores (pore size > 50 nm). D 1 and D 2 are positively correlated with the SSA and PV but negatively correlated with the average pore size. D 3 is positively correlated with the connectivity and average pore size but irrelevant to the PV, SSA, and average pore size. With the increase in the TOC, the D 1 , D 2 , PV, and SSA decreased, but the connectivity improved slightly. The shales of the different lithofacies show the order of D 3 > D 2 > D 1 . The pore structure is more complex than the pore surface, and the heterogeneity of the large pores is greater compared to the small pores. Among the three lithofacies, the CS had the largest PV, SSA, D 1 , and D 2 , indicating the development of a more complex pore structure network, which expands the space required for shale oil occurrence. However, the connectivity of the CS was the worst among the three lithofacies, which is adverse to shale oil production. Although the PV of the SS was slightly smaller than that of the CS, the average pore size and connectivity are obviously more desirable, which makes the SS an ideal shale reservoir.

Data Availability Statement:
The data presented in this study are available upon request from the corresponding author. The data are not publicly available because some of the basic research involves confidentiality.