Pore Characteristics of Lacustrine Shale Oil Reservoir in the Cretaceous Qingshankou Formation of the Songliao Basin, NE China

: Shale oil is hosted in nanopores of organic-rich shales, so pore characteristics are signiﬁcant for shale oil accumulation. Here we analyzed pore characteristics of 39 lacustrine shale samples of the Late Cretaceous Qingshankou Formation (K 2 qn) in the Songliao Basin, which is one of the main shale oil resource basins in China, using ﬁeld emission-scanning electron microscopy (FE-SEM), and low-pressure nitrogen adsorption. We accomplished fractal analysis, correlation analysis using correlation matrix and multidimensional scaling (MDS), and prediction of fractal dimensions, which is the ﬁrst time to predict pore fractal dimensions of shales. Interparticle pores are highly developed in K 2 qn. These shales have mesoporous nature and slit-shaped pores. Compared with the second and third members (K 2 qn 2,3 ), the ﬁrst member of the Qingshankou Formation (K 2 qn 1 ) has a larger average pore diameter, much smaller surface area, fewer micropores, simpler pore structure and surface indicated by smaller fractal dimensions. In terms of pore characteristics, K 2 qn 1 is better than K 2 qn 2,3 as a shale oil reservoir. When compared with marine Bakken Formation shales, lacustrine shales of the Qingshankou Formation have similar complexity of pore structure, but much rougher pore surface. This research can lead to an improved understanding of the pore system of lacustrine shales.


Introduction
With the increasing energy demands of the world, the transition from conventional petroleum to unconventional petroleum is inevitable [1][2][3][4][5][6]. Shale oil, as a kind of unconventional petroleum resource, is attracting growing attention [7][8][9][10]. Multistage-fractured horizontal drilling dramatically enhances the production performance of shale reservoirs [11]. Natural fracture spacing significantly affects shale reservoir performance [11]. It is more profitable to find oil instead of gas since oil has more energy density and it is easier to transport, store, and use [9].
A shale oil reservoir is a kind of self-contained source-reservoir system. The pores and micro-fractures in inorganic minerals and organic matter constitute a complicated pore-fracture system in shales [1,2,12,13]. Shale oil is stored in four types of pores, i.e., intraparticle (intraP) pores, direction with a diamond-shaped basin floor. It went through three major tectonic events, including rifting, thermal subsidence, and structural inversion [46]. During the Late Mesozoic, the interaction between the Pacific and Eurasian plates resulted in an extended rift-valley system in south Mongolia and northeast China [47,48]. From the Late Jurassic to the Cretaceous, the Songliao Basin, which was in northeast China, was developed as a rift basin in this extended system [46][47][48]. During thermal subsidence, lithospheric cooling and extension during the whole Cretaceous resulted in regional deformation and subsidence [46]. The Songliao Basin gently developed into a considerable lake with lake-level fluctuations and probable marine incursions [49]. Generally, the thickness of typical sediment packages is about 3000-4000 m. These sediment packages are composed of lacustrine, alluvial, and fluvial sediments [45]. Sediments in the SK-1 well cores were deposited during this stage. Structural inversion, which included uplift and folding, started with the deposition of the Upper Nenjiang Formation and culminated at the end of the Cretaceous [46]. Six structural units have been established in the Songliao Basin. They are Central Depression Zone, North Plunge Zone, Northeast Uplift Zone, Southeast Uplift Zone, Southwest Uplift Zone, and West Slope Zone, respectively (Figure 1b) [50].
Energies 2020, 13, x FOR PEER REVIEw 3 of 25 between the Pacific and Eurasian plates resulted in an extended rift-valley system in south Mongolia and northeast China [47,48]. From the Late Jurassic to the Cretaceous, the Songliao Basin, which was in northeast China, was developed as a rift basin in this extended system [46][47][48]. During thermal subsidence, lithospheric cooling and extension during the whole Cretaceous resulted in regional deformation and subsidence [46]. The Songliao Basin gently developed into a considerable lake with lake-level fluctuations and probable marine incursions [49]. Generally, the thickness of typical sediment packages is about 3000-4000 m. These sediment packages are composed of lacustrine, alluvial, and fluvial sediments [45]. Sediments in the SK-1 well cores were deposited during this stage. Structural inversion, which included uplift and folding, started with the deposition of the Upper Nenjiang Formation and culminated at the end of the Cretaceous [46]. Six [41]); (b) major structural zones in the Songliao Basin (modified after [41]); (c) structural cross-section across the central part of the Songliao Basin (modified after [46]); and (d) stratigraphic histogram of Qingshankou Formation in the SK-1s well (modified after [43]) and sampling depths.
The SK-1 drilling project includes two boreholes, the north borehole (SK-1n) and the south borehole (SK-1s). Figure 1b shows the locations of these two boreholes. These two boreholes are 60 km apart. The six stratigraphic formations in the SK-1 cores (Figure 1c) are primarily composed of mudstones, siltstones, and sandstones, deposited in various paleoenvironments from deep lakes to flood plains [45]. The Qingshankou Formation (1285.91-1782.93 m in SK-1s well) is primarily composed of black-gray mudstone interbedded with gray sandstone, siltstone, and dolomite. The black mudstone in the first member of the Qingshankou Formation (K2qn 1 ), which is also the most significant petroleum source rock in the Daqing Oilfield, was deposited in a deep lake environment [46]. The black-gray mudstone in the second and third members (K2qn 2,3 ) was deposited mainly in a semi-deep lake environment [50]. The Qingshankou Formation is Late Turonian to Late Coniacian in age [50].  [41]); (b) major structural zones in the Songliao Basin (modified after [41]); (c) structural cross-section across the central part of the Songliao Basin (modified after [46]); and (d) stratigraphic histogram of Qingshankou Formation in the SK-1s well (modified after [43]) and sampling depths.
The SK-1 drilling project includes two boreholes, the north borehole (SK-1n) and the south borehole (SK-1s). Figure 1b shows the locations of these two boreholes. These two boreholes are 60 km apart. The six stratigraphic formations in the SK-1 cores (Figure 1c) are primarily composed of mudstones, siltstones, and sandstones, deposited in various paleoenvironments from deep lakes to flood plains [45]. The Qingshankou Formation (1285.91-1782.93 m in SK-1s well) is primarily composed of black-gray mudstone interbedded with gray sandstone, siltstone, and dolomite. The black mudstone in the first member of the Qingshankou Formation (K 2 qn 1 ), which is also the most significant petroleum source rock in the Daqing Oilfield, was deposited in a deep lake environment [46]. The black-gray mudstone in the second and third members (K 2 qn 2,3 ) was deposited mainly in a semi-deep lake environment [50]. The Qingshankou Formation is Late Turonian to Late Coniacian in age [50].

Sampling
A total of 39 lacustrine shale samples were selected from K 2 qn in the SK-1s well. Among them, 20 samples were from K 2 qn 1 , and 19 samples were from K 2 qn 2,3 . Figure 1d shows the strata and sampling depths. Since these core samples were well preserved without weathering, the measured data is accurate. Each sample was divided into several pieces for conducting low-pressure nitrogen adsorption, and FE-SEM studies.

Field Emission-Scanning Electron Microscopy (FE-SEM)
We used an argon ion beam to polish one piece of each sample (about 10 mm in length, 10 mm in width, and 3-5 mm in height) to generate a flat and smooth surface by Leica RES102 Ion Milling System of Leica company from Germany at the Petroleum Geology Research and Laboratory Center. Every single sample was carbon-coated (13 nm thick). An Apreo FE-SEM of FEI company then imaged the samples. The maximum resolution of a single-pixel is 0.9 nm. Energy-dispersive spectrometer (EDS) measured the semi-quantitative analysis of mineral composition. We can directly observe pore types, pore morphology, and pore distribution [17,19].

Low-pressure Nitrogen Adsorption
Low-pressure nitrogen adsorption was conducted with a Micromeritics ASAP 2020HD88 Surface Area and Porosity Analyzer of Micromeritics company from the U.S. at Beijing Center for Physical and Chemical Analysis, China. Equivalent surface areas were determined with the Brunauer-Emmett-Teller (BET) method. Pore volumes were calculated using the Density Functional Theory (DFT) model. Average pore diameters were calculated with the adsorption branch by the Barrette-Joyner-Halenda (BJH) model. There is a comprehensive description of the different theories mentioned above [51]. In this study, with the DFT model, micropores with pore size from 1.6 nm to 2 nm, mesopores with pore size from 2 nm to 50 nm, and macropores with pore size from 50 nm to 148 nm were detected [52].

Pore Morphology from FE-SEM
Four types of pores were recognized in this study, including organic-matter pores, interparticle (interP) pores, intraparticle (intraP) pores, and micro-fractures.

Organic-Matter Pores
Organic-matter pores in the Qingshankou Formation are moderately developed and show high heterogeneity. Some organic-matter pores are elongated, with maximum pore size up to 10 µm long and 2.5 µm wide (Figure 2a-c). Some organic-matter pores are oval or round, with maximum pore size up to 15 µm in diameter (Figure 2d). There are also some honeycomb-like organic-matter pores, with a wide range of pore sizes between 0.1 µm and 7 µm (Figure 2e,f). Some irregular band-like organic matter with clear and sharp boundaries does not have visible pores (Figure 2g). Nevertheless, when band-like organic matter concentrates together and is in contact with each other, a small number of interparticle pores form with pore size up to 5 µm long and 1 µm wide (Figure 2h). Some gridded-like pores are found within organic matter, which coexists with clay minerals, calcite grains, and pyrite particles (Figure 2i). Some pores are also found within organic matter, and are surrounded by pyrite circle (Figure 2j) or coexist with pyrite particles within pyrite circle (Figure 2k). Some pores can be observed within the organic matter in pyrite framboids (Figure 2l).

Interparticle Pores
Interparticle pores are highly developed and are the most significant pore type in the Qingshankou Formation. Some interP pores can be found between organic matter and inorganic minerals (Figure 3a). Some grain-edge interP pores exist among pyrite particles, clay minerals, and other minerals (Figure 3b). Some interP pores are present around rigid grains (Figure 3c and d). There are also many interP pores between clay platelets (Figure 3e and f).

Interparticle Pores
Interparticle pores are highly developed and are the most significant pore type in the Qingshankou Formation. Some interP pores can be found between organic matter and inorganic minerals (Figure 3a). Some grain-edge interP pores exist among pyrite particles, clay minerals, and other minerals (Figure 3b). Some interP pores are present around rigid grains (Figure 3c,d). There are also many interP pores between clay platelets (Figure 3e,f).

Intraparticle Pores
Intraparticle pores are found within particles. Some of them are intercrystalline pores, but intragranular pores also exist. Intercrystalline intraP pores within pyrite framboids can be easily observed in the Qingshankou Formation, with pore size up to 0.5 μm (Figure 4a). Cleavage-sheet and cleavage-wedge intraP pores are also common within clay aggregates (Figure 4b and c). There are also some intragranular pores in dolomite grains, with a maximum pore size up to 0.5 μm (Figure 4b and c). Lots of dissolution-rim pores are observed in dolomite grains (Figure 4d). Some intragranular pores exist in K-feldspar grains, with pore size up to 0.8 μm (Figure 4e). Some intragranular pores are observed within plagioclase grains (Figure 4f).

Micro-fractures
Micro-fractures, which are a crucial part of an active pore network, can be found in the Qingshankou Formation. Some micro-fractures with a length of more than 400 μm and a width of 8

Intraparticle Pores
Intraparticle pores are found within particles. Some of them are intercrystalline pores, but intragranular pores also exist. Intercrystalline intraP pores within pyrite framboids can be easily observed in the Qingshankou Formation, with pore size up to 0.5 µm (Figure 4a). Cleavage-sheet and cleavage-wedge intraP pores are also common within clay aggregates (Figure 4b,c). There are also some intragranular pores in dolomite grains, with a maximum pore size up to 0.5 µm (Figure 4b

Intraparticle Pores
Intraparticle pores are found within particles. Some of them are intercrystalline pores, but intragranular pores also exist. Intercrystalline intraP pores within pyrite framboids can be easily observed in the Qingshankou Formation, with pore size up to 0.5 μm (Figure 4a). Cleavage-sheet and cleavage-wedge intraP pores are also common within clay aggregates (Figure 4b and c). There are also some intragranular pores in dolomite grains, with a maximum pore size up to 0.5 μm (Figure 4b

Micro-fractures
Micro-fractures, which are a crucial part of an active pore network, can be found in the Qingshankou Formation. Some micro-fractures with a length of more than 400 μm and a width of 8

Micro-Fractures
Micro-fractures, which are a crucial part of an active pore network, can be found in the Qingshankou Formation. Some micro-fractures with a length of more than 400 µm and a width of 8 µm that are filled with organic matter are also observed ( Figure 5a). Some micro-fractures are not filled with organic matter (Figure 5a-c).
Energies 2020, 13, x FOR PEER REVIEw 7 of 25 μm that are filled with organic matter are also observed ( Figure 5a). Some micro-fractures are not filled with organic matter (Figure 5a-c).

Nature of Isotherms and Hysteresis Loops
The nitrogen adsorption/desorption isotherms for samples in K2qn 1 and K2qn 2,3 are shown in Figures 6a and b, respectively. Adsorption/desorption isotherms show obvious hysteresis loops, which are visible when relative pressure is about 0.45 ( Figure 6). This phenomenon is called the 'forced closure' of the desorption branch due to the Tensile Strength Effect [53]. Various interpretations have been proposed for hysteresis [54][55][56], with a possible interpretation being capillary condensation in mesoporous structures [56]. All the studied samples are of Type IV [56]. Type IV isotherm is characterized by its hysteresis loop, which is related to capillary condensation taking place in mesopores, and the limiting uptake over a range of high relative pressure. Due to monolayer-multilayer adsorption, the first part of the Type IV isotherm follows the same path as the corresponding part of a Type II isotherm. This indicates that mesopores are the dominating pores of the samples in both K2qn 1 and K2qn 2,3 . Most of the isotherms in K2qn 1 and K2qn 2,3 rise rapidly or even vertically when relative pressure is close to 1 (Figure 6a and b), which indicates these samples might have macropores or fractures [29,56].

Nature of Isotherms and Hysteresis Loops
The nitrogen adsorption/desorption isotherms for samples in K 2 qn 1 and K 2 qn 2,3 are shown in Figure 6a,b, respectively. Adsorption/desorption isotherms show obvious hysteresis loops, which are visible when relative pressure is about 0.45 ( Figure 6). This phenomenon is called the 'forced closure' of the desorption branch due to the Tensile Strength Effect [53]. Various interpretations have been proposed for hysteresis [54][55][56], with a possible interpretation being capillary condensation in mesoporous structures [56]. All the studied samples are of Type IV [56]. Type IV isotherm is characterized by its hysteresis loop, which is related to capillary condensation taking place in mesopores, and the limiting uptake over a range of high relative pressure. Due to monolayer-multilayer adsorption, the first part of the Type IV isotherm follows the same path as the corresponding part of a Type II isotherm. This indicates that mesopores are the dominating pores of the samples in both K 2 qn 1 and K 2 qn 2,3 . Most of the isotherms in K 2 qn 1 and K 2 qn 2,3 rise rapidly or even vertically when relative pressure is close to 1 (Figure 6a,b), which indicates these samples might have macropores or fractures [29,56].

Nature of Isotherms and Hysteresis Loops
The nitrogen adsorption/desorption isotherms for samples in K2qn 1 and K2qn 2,3 are shown in Figures 6a and b, respectively. Adsorption/desorption isotherms show obvious hysteresis loops, which are visible when relative pressure is about 0.45 ( Figure 6). This phenomenon is called the 'forced closure' of the desorption branch due to the Tensile Strength Effect [53]. Various interpretations have been proposed for hysteresis [54][55][56], with a possible interpretation being capillary condensation in mesoporous structures [56]. All the studied samples are of Type IV [56]. Type IV isotherm is characterized by its hysteresis loop, which is related to capillary condensation taking place in mesopores, and the limiting uptake over a range of high relative pressure. Due to monolayer-multilayer adsorption, the first part of the Type IV isotherm follows the same path as the corresponding part of a Type II isotherm. This indicates that mesopores are the dominating pores of the samples in both K2qn 1 and K2qn 2,3 . Most of the isotherms in K2qn 1 and K2qn 2,3 rise rapidly or even vertically when relative pressure is close to 1 (Figure 6a and b), which indicates these samples might have macropores or fractures [29,56].  Though these samples display isotherms with hysteresis patterns, these isotherms show some differences in shape. Some isotherms show steep slopes when relative pressure is close to 1, and others show moderate or gentler slopes. The slope of the isotherms indicates the development of macropores [35]. The parameter ∆V G , which is the difference in volumes of gas adsorbed at the last two highest relative pressures recorded, was used to calculate the slope of the isotherms [57]. However, the slope should be calculated with the difference in volumes of gas adsorbed at the last two highest relative pressures recorded divided by the difference in the last two highest relative pressures recorded (∆V G /∆(P/P 0 )) in Figure 6. Therefore, we used the parameter ∆V G /∆(P/P 0 ) instead of the parameter ∆V G to represent the slope of the isotherms in this study.
All the ∆V G /∆(P/P 0 ) values are shown in Table A1. Some samples, such as SK1S04, SK1S05, SK1S12, and SK1S38, show high ∆V G /∆(P/P 0 ) values, and correspondingly show steep isotherms (which was interpreted to be indicative of macropores [35]) at the highest two relative pressure steps ( Figure 6). However, the presence of apparent hysteresis of these shales also represents the existence of mesopores. Some samples, such as SK1S28 and SK1S37, show low ∆V G /∆(P/P 0 ) values, as the slopes are gentle. Along with the hysteresis patterns, it indicates that these samples are mainly mesoporous. Other samples, such as SK1S01, SK1S11, SK1S18, SK1S26, and SK1S39, are marked by the existence of both mesopores (visible hysteresis pattern) and macropores (lower in concentration relative to samples with high ∆V G /∆(P/P 0 ) values).
According to the types of hysteresis loops [56], the hysteresis patterns shown by the samples are mainly of two types: (1) H3, such as SK1S01, SK1S04, SK1S05, SK1S11, SK1S12, SK1S18, SK1S26, SK1S38, and SK1S39; and (2) the combination of H3-H4, such as SK1S28 and SK1S37. The H3 hysteresis pattern, which does not show any limiting adsorption at high relative pressure, is observed with aggregates of plate-like particles producing slit-shaped pores [56]. The H4 hysteresis pattern is typically related to porous materials containing narrow slit-shaped pores [56]. However, explaining the pore shapes with hysteresis patterns might be not accurate [29], since the actual pore shape might be a mix of various pore types, which have also been observed in the FE-SEM images.

Pore Structure Parameters
The quantitative pore structure parameters of the studied samples from low-pressure nitrogen adsorption are shown in Table A1. The statistic values (minimum value, average value, and maximum value) of pore structure parameters in K 2 qn 1 and K 2 qn 2,3 are shown in Figure 7. When comparing pore structure parameters and fractal dimensions of two different members (K 2 qn 1 and K 2 qn 2,3 ), the average value of each member can represent the overall level of each member. Therefore, the average value has been applied here. The average pore diameter in K 2 qn varies from 7.34 nm to 18.22 nm, with an average of 11.89 nm. The average pore diameter in K 2 qn 2,3 (9.42 nm on average) is smaller than that in K 2 qn 1 (14.24 nm on average) (Figure 7a). BET specific surface area (BET SSA) in K 2 qn exhibits a wide range of values from 3.04 m 2 /g to 64.57 m 2 /g, with an average of 22.70 m 2 /g. BET SSA in K 2 qn 2,3 (36.75 m 2 /g on average) is enormously more massive than that in K 2 qn 1 (9.36 m 2 /g on average) (Figure 7b). Total pore volume in K 2 qn is in the range of 12.13-65.40 × 10 −3 cm 3 /g, with an average of 38.35 × 10 −3 cm 3 /g. Total pore volume in K 2 qn 2,3 (44.73 × 10 -−3 cm 3 /g on average) is bigger than that in K 2 qn 1 (32.30 × 10 −3 cm 3 /g on average) (Figure 7c). Micropore volume in K 2 qn is in the range of 0-8.52 × 10 −3 cm 3 /g, with an average of 1.90 × 10 −3 cm 3 /g. Micropore volume in K 2 qn 2,3 (3.85 × 10 −3 cm 3 /g on average) is dramatically more massive than that in K 2 qn 1 (0.04 × 10 −3 cm 3 /g on average). There is almost no micropore in K 2 qn 1 (Figure 7d). Mesopore volume in K 2 qn varies from 7.77 × 10 −3 cm 3 /g to 47.50 × 10 −3 cm 3 /g, with an average of 28.75 × 10 −3 cm 3 /g. Mesopore volume accounts for 74.97% of total pore volume, which indicates mesopore is the dominating pore type in K 2 qn. Mesopore volume in K 2 qn 2,3 (33.39 × 10 −3 cm 3 /g on average) is larger than that in K 2 qn 1 (24.33 × 10 −3 cm 3 /g on average) (Figure 7e). Macropore volume in K 2 qn 1 (7.93 × 10 −3 cm 3 /g on average) is a little bigger than that in K 2 qn 2,3 (7.48 × 10 −3 cm 3 /g on average), with an average of 7.71 × 10 −3 cm 3 /g in K 2 qn (Figure 7f).
These pore structure parameters in K 2 qn are much larger than those of marine shales from the Early Cambrian Niutitang Formation in the Sichuan Basin, southwest China [22].

Fractal Dimensions from Nitrogen Adsorption Isotherms
The fractal Frenkel-Halsey-Hill (FHH) model based on gas adsorption isotherms has been proven to be a useful model and is broadly used in shales [58]. The FHH model involves applying Equation (1)

Fractal Dimensions from Nitrogen Adsorption Isotherms
The fractal Frenkel-Halsey-Hill (FHH) model based on gas adsorption isotherms has been proven to be a useful model and is broadly used in shales [58]. The FHH model involves applying Equation (1) [59,60]: where P is the equilibrium pressure; P 0 is the saturation pressure of the gas; V is the volume of adsorbed gas molecules; V 0 is the volume of monolayer coverage; A is the power-law exponent, which depends on the fractal dimension (D) and the mechanisms of adsorption. D value can be calculated from the slope (A) of the straight line in the ln V versus ln(ln(P 0 /P)) FHH plot using either Equations (2) or (3): Both Equations (2) and (3) have been applied to calculate the fractal dimension D [59,61]. Equation (2) is applied in the capillary condensation regime in which the hysteresis loops start to appear, while Equation (3) is used in the case of the van der Waals force [62].
Following the FHH model, the plots of ln V versus ln(ln(P 0 /P)) in K 2 qn 1 and K 2 qn 2,3 are shown in Figure 8a,b, respectively. There are two different linear segments at the P/P 0 intervals of 0-0.5 and 0.5-1, which indicates that fractal characteristics exist in these two distinct regions [63]. The fractal dimensions in these two distinct regions represent different types of pore features. D 1 is at higher P/P 0 and represents the pore structure fractal dimension, whose behavior is controlled by capillary condensation. D 2 is at lower P/P 0 and represents the pore surface fractal dimension, whose adsorption behavior is called monolayer-multilayer adsorption controlled by van der Waals forces [60].
where is the equilibrium pressure; 0 is the saturation pressure of the gas; is the volume of adsorbed gas molecules; 0 is the volume of monolayer coverage; is the power-law exponent, which depends on the fractal dimension (D) and the mechanisms of adsorption. D value can be calculated from the slope ( ) of the straight line in the ln versus ln(ln( 0 ⁄ )) FHH plot using either Equations (2) or (3): Both Equations (2) and (3) have been applied to calculate the fractal dimension D [59,61]. Equation (2) is applied in the capillary condensation regime in which the hysteresis loops start to appear, while Equation (3) is used in the case of the van der Waals force [62].
Following the FHH model, the plots of ln versus ln(ln( 0 ⁄ )) in K2qn 1 and K2qn 2,3 are shown in Figures 8a and b, respectively. There are two different linear segments at the P/P0 intervals of 0-0.5 and 0.5-1, which indicates that fractal characteristics exist in these two distinct regions [63]. The fractal dimensions in these two distinct regions represent different types of pore features. D1 is at higher P/P0 and represents the pore structure fractal dimension, whose behavior is controlled by capillary condensation. D2 is at lower P/P0 and represents the pore surface fractal dimension, whose adsorption behavior is called monolayer-multilayer adsorption controlled by van der Waals forces [60]. In general, the fractal dimension (D) varies in the range of 2-3. It is affected by the geometrical irregularities and roughness of the surface [64]. The larger the D values are, the more complicated and irregular the surfaces are [29]. Table A2 shows fractal fitting parameters and fractal dimensions of the studied samples. Using Equation (2), D1 varies from 2.4507 to 2.7498 (2.5982 on average), which meets the conditions of the fractal theory of porous structures [65]. However, D1 calculated by Equation (3) varies from 1.3521 to 2.2494, which shows deviations from the definition of fractal dimension, so it is unreasonable and useless [65]. The same thing happens to D2. The fractal dimension varies from 2 for an entirely typical smooth solid surface to 3 for an extremely complicated surface. In this situation, only Equation (2) is used to calculate D1 and D2 in this study.
Both segments have good linear relationships, and all the coefficients of determination (R1 2 and R2 2 ) are higher than 0.9775 (Table A2), indicating that shales in K2qn are fractal [64]. Comparing with the fractal dimensions in these two formations, the values of D1 in K2qn 1 have a narrower distribution ranging from 2.4507 to 2.5600 with an average of 2.5196, while the values of D1 in K2qn 2,3 range between 2.5669 and 2.7498 with an average of 2.6809 (Figure 7g). The lower average value of D1 indicates that shales in K2qn 1 have simpler pore structures than those in K2qn 2,3 . The values of D2 in In general, the fractal dimension (D) varies in the range of 2-3. It is affected by the geometrical irregularities and roughness of the surface [64]. The larger the D values are, the more complicated and irregular the surfaces are [29]. Table A2 shows fractal fitting parameters and fractal dimensions of the studied samples. Using Equation (2), D 1 varies from 2.4507 to 2.7498 (2.5982 on average), which meets the conditions of the fractal theory of porous structures [65]. However, D 1 calculated by Equation (3) varies from 1.3521 to 2.2494, which shows deviations from the definition of fractal dimension, so it is unreasonable and useless [65]. The same thing happens to D 2 . The fractal dimension varies from 2 for an entirely typical smooth solid surface to 3 for an extremely complicated surface. In this situation, only Equation (2) Figure 7h). It suggests shales in K 2 qn 1 have smoother pore surfaces than those in K 2 qn 2,3 .

Relationships between Pore Structure Parameters
Multidimensional scaling (MDS) [66][67][68] is a method used in data sciences to visualize and compare similarities and differences of high dimensional data. Its application in geostatistics helps visually evaluate and interpret multivariate data in a lower dimension. MDS does this by projecting the multivariate distances between entities to a best-fit configuration in lower dimensions that we can see [69]. More details can be seen in Reference [69]. To investigate the relationships between pore structure parameters in K 2 qn, we used MDS in this study.
The correlation matrix and multidimensional scaling of multivariate pore structure parameters are shown in Figure 9a and b, respectively. The correlation matrix shows the Pearson correlation coefficients between these pore structure parameters (Figure 9a). The Pearson correlation coefficients vary from −1 to 1 and represent from extremely negatively correlated to extraordinarily positively correlated. The multidimensional scaling plot shows the degree of correlation between these parameters (Figure 9b). The distances of these parameters are based on three coordinates, X 1 (shown as horizontal axis), X 2 (shown as vertical axis), and X 3 (shown as color from cold to hot), respectively. The closer the two parameters are observed in the MDS plot, the more positive correlation exists between them, and vice versa.  Figure 7h). It suggests shales in K2qn 1 have smoother pore surfaces than those in K2qn 2,3 .

Relationships between Pore Structure Parameters
Multidimensional scaling (MDS) [66][67][68] is a method used in data sciences to visualize and compare similarities and differences of high dimensional data. Its application in geostatistics helps visually evaluate and interpret multivariate data in a lower dimension. MDS does this by projecting the multivariate distances between entities to a best-fit configuration in lower dimensions that we can see [69]. More details can be seen in Reference [69]. To investigate the relationships between pore structure parameters in K2qn, we used MDS in this study.
The correlation matrix and multidimensional scaling of multivariate pore structure parameters are shown in Figure 9a and b, respectively. The correlation matrix shows the Pearson correlation coefficients between these pore structure parameters (Figure 9a). The Pearson correlation coefficients vary from −1 to 1 and represent from extremely negatively correlated to extraordinarily positively correlated. The multidimensional scaling plot shows the degree of correlation between these parameters (Figure 9b). The distances of these parameters are based on three coordinates, X1 (shown as horizontal axis), X2 (shown as vertical axis), and X3 (shown as color from cold to hot), respectively. The closer the two parameters are observed in the MDS plot, the more positive correlation exists between them, and vice versa.  We analyzed the correlations among average pore diameter, BET SSA, total pore volume, micropore volume, mesopore volume, and macropore volume in this part. From the MDS plot in Figure 9b, we can get a general impression that BET SSA and micropore volume are more positively related because these two dots are close to each other, and total pore volume and mesopore volume are more positively related for the same reason.
From the MDS plot in Figure 9b, we can tell that the dot representing the average pore diameter is quite far away from other parameters, except the dot representing macropore volume. The average pore diameter has a weak positive correlation with macropore volume, with a correlation coefficient of 0.24 (Figure 9a), which is reasonable. The two dots representing BET SSA and micropore volume are quite far away from the dot representing the average pore diameter in Figure 9b, with correlation coefficients of −0.83 and −0.73 (Figure 9a), which shows they have strong negative correlations with average pore diameter. The more micropores the samples have, the smaller the average pore diameter they have, and the larger BET SSA they have, which is reasonable. Total pore volume and mesopore volume are flocking together in Figure 9b and medium negatively correlated with average pore diameter, with correlation coefficients of −0.58 and −0.59 (Figure 9a), respectively. In brief, the most important parameters influencing average pore diameter are BET SSA and micropore volume, which show strongly negative correlations with the average pore diameter.
The BET SSA is close to micropore volume in Figure 9b, and the correlation coefficient between them is 0.98 (Figure 9a), which means these two parameters are extremely positively correlated. The dot representing BET SSA is a little bit far away from the dot representing mesopore volume compared to the dot representing micropore volume in Figure 9b, and the correlation coefficient between BET SSA and mesopore volume is 0.59 (Figure 9a), which means these two parameters are medium positively correlated. The dot representing BET SSA is far away from the dot representing macropore volume in Figure 9b, with a correlation coefficient of −0.14 ( Figure 9a), which shows these two parameters barely correlate. Since micropore is smaller than mesopore, it has more BET SSA than mesopore. For the same reason, mesopore has more BET SSA than macropore. In brief, the most critical parameter influencing BET SSA is micropore volume, which is positively correlated with BET SSA.
The total pore volume and the mesopore volume are so close to each other in Figure 9b, and they have an exceedingly positive correlation (correlation coefficient 0.98) (Figure 9a), which indicates that mesopore is the primary pore type in K 2 qn. This result is consistent with the results mentioned earlier that mesopore volume accounts for 74.97% of the total pore volume.

Relationships between Fractal Dimensions and Pore Structure Parameters
From the MDS plot in Figure 9b, we can find out D 1 , D 2 , BET SSA, and micropore volume cluster together, which indicates that these four parameters are strongly positively related. Two fractal dimensions (D 1 and D 2 ) are quite close to each other, and they have a correlation coefficient of 0.89 (Figure 9a), which means D 1 and D 2 have a significantly positive correlation. This result is different from the inconsistency in D 1 and D 2 of coals [60], which might indicate lacustrine shales with more complicated pore structures usually have rougher surfaces.
Both two fractal dimensions (D 1 and D 2 ) have strongly positive correlations (correlation coefficients larger than 0.85) with BET SSA and micropore volume, which indicates more micropores and larger BET SSA lead to more complex pore structures and much rougher surfaces. Both two fractal dimensions (D 1 and D 2 ) are moderately close to the mesopore volume in Figure 9b, which means mesopore makes a moderate contribution to the complexity of the pore structure and the roughness of the pore surface. These two dots representing D 1 and D 2 are far away from the dot representing macropore volume in Figure 9b, with correlation coefficients of −0.26 and −0.23 (Figure 9a). This indicates macropore makes some contribution to simpler pore structure and smoother pore surface. These two dots representing D 1 and D 2 are exceedingly far away from the dot representing the average pore diameter in Figure 9b, with correlation coefficients of −0.96 and −0.80 (Figure 9a), which shows larger pores are useful for making pore structure simpler and pore surface smoother.

Prediction of Fractal Dimensions with Pore Structure Parameters
Multivariate linear regression is a kind of traditional regression technique of response surface modeling (RSM). Multivariate linear regression can be applied to predict a single dependent response variable using multiple independent predictor variables [22]. Stepwise regression, which is a kind of multivariate linear regression, has been used in this study to predict fractal dimensions with pore structure parameters.
Five independent variables were used in this study, including average pore diameter, BET SSA, micropore volume, mesopore volume, and macropore volume. In the criterion section, F-to-Enter is set to 4.0 (default value), and F-to-Remove is set to 3.9 (default value). P-value to reject is set to 0.05 (default value). Confidence intervals is set to 95% (default value). Generally, significance tests should be conducted on the regression equation (mainly including the test of goodness of fit and F test) and regression coefficients. However, the F test on the regression equation and the significance test on the regression coefficients are included during the regression process [22]. Therefore, only the test of goodness of fit is required. Table 1 lists the statistics from the regression models and the variance analysis results. The coefficient of determination, which is the square of multiple correlation coefficient, can be applied to judge the fitting degree of linear regression and illustrate how independent variables can explain the variations of dependent variables. The coefficient of multiple determination, which is the corrected coefficient of determination, is used to remove the effects produced by the number of explanatory variables. In Model 1, as shown in Table 1, the coefficient of multiple determination is 0.960, and the P-value is less than 0.001, which is much less than the default value (0.05). It suggests that Model 1 has passed the significance test. The regression equation of Model 1 can be expressed as: where y denotes D 1 ; x 1 denotes average pore diameter (nm); x 2 denotes BET SSA (m 2 /g). Figure 10a shows the 3D scatter plot and the response surface of Model 1. D 1 shows a negative correlation with the average pore diameter and a positive correlation with BET SSA, which is consistent with the result mentioned above.

Prediction of Fractal Dimensions with Pore Structure Parameters
Multivariate linear regression is a kind of traditional regression technique of response surface modeling (RSM). Multivariate linear regression can be applied to predict a single dependent response variable using multiple independent predictor variables [22]. Stepwise regression, which is a kind of multivariate linear regression, has been used in this study to predict fractal dimensions with pore structure parameters.
Five independent variables were used in this study, including average pore diameter, BET SSA, micropore volume, mesopore volume, and macropore volume. In the criterion section, F-to-Enter is set to 4.0 (default value), and F-to-Remove is set to 3.9 (default value). P-value to reject is set to 0.05 (default value). Confidence intervals is set to 95% (default value). Generally, significance tests should be conducted on the regression equation (mainly including the test of goodness of fit and F test) and regression coefficients. However, the F test on the regression equation and the significance test on the regression coefficients are included during the regression process [22]. Therefore, only the test of goodness of fit is required. Table 1 lists the statistics from the regression models and the variance analysis results. The coefficient of determination, which is the square of multiple correlation coefficient, can be applied to judge the fitting degree of linear regression and illustrate how independent variables can explain the variations of dependent variables. The coefficient of multiple determination, which is the corrected coefficient of determination, is used to remove the effects produced by the number of explanatory variables. In Model 1, as shown in Table 1, the coefficient of multiple determination is 0.960, and the Pvalue is less than 0.001, which is much less than the default value (0.05). It suggests that Model 1 has passed the significance test. The regression equation of Model 1 can be expressed as: = −0.0214 1 + 0.0018 2 + 2.8110 (4) where denotes D1; 1 denotes average pore diameter (nm); 2 denotes BET SSA (m 2 /g). Figure 10a shows the 3D scatter plot and the response surface of Model 1. D1 shows a negative correlation with the average pore diameter and a positive correlation with BET SSA, which is consistent with the result mentioned above.

Prediction of D 2
In Model 2, as shown in Table 1, the coefficient of multiple determination is 0.895, and the P-value is less than 0.001, which is much less than the default value (0.05). It suggests that Model 2 has passed the significance test. The regression equation of Model 2 can be expressed as: where y denotes D 2 ; x 1 denotes micropore volume (10 −3 cm 3 /g); x 2 denotes average pore diameter (nm). Figure 10b shows the 3D scatter plot and the response surface of Model 2. D 2 shows a positive correlation with micropore volume and a negative correlation with average pore diameter, which is consistent with the result mentioned above. Table A1 shows the free hydrocarbon content (S 1 ). S 1 of K 2 qn 1 is ranging from 0.06 mg HC/g to 4.30 mg HC/g, with an average of 1.95 mg HC/g. However, S 1 of K 2 qn 2,3 is in the range of 0.05-1.46 mg HC/g, with an average of 0.44 mg HC/g, which is much smaller than that of K 2 qn 1 . Therefore, K 2 qn 1 is much better than K 2 qn 2,3 in terms of free hydrocarbon content. Figure 11 shows the correlations between fractal dimensions (D 1 and D 2 ) and free hydrocarbon content (S 1 ). S 1 of lacustrine shale samples is negatively correlated with D 1 and D 2 . The coefficient of determination (r 2 ) between S 1 and D 1 is more significant than that between S 1 and D 2 , which demonstrates that D 1 can be applied to characterize the shale oil content. Shale samples with more free oil usually have a smaller D 1 , probably because free oil prefers being stored in larger pores, and a simpler pore structure makes oil flow smoother. This is quite distinct from coal and gas shale, which reveal that gas adsorption capacity is positively correlated with fractal dimensions [21]. It is because rougher surfaces (with larger fractal dimensions) can supply more adsorption sites for gas molecules. However, oil molecules, which are much bigger than gas molecules, only can be stored in larger pores (with smaller fractal dimensions).

Relationships between Free Hydrocarbon Content and Fractal Dimensions
Energies 2020, 13, x FOR PEER REVIEw 14 of 25

Prediction of D2
In Model 2, as shown in Table 1, the coefficient of multiple determination is 0.895, and the Pvalue is less than 0.001, which is much less than the default value (0.05). It suggests that Model 2 has passed the significance test. The regression equation of Model 2 can be expressed as: = 0.0186 1 − 0.0061 2 + 2.5150 (5) where denotes D2; 1 denotes micropore volume (10 −3 cm 3 /g); 2 denotes average pore diameter (nm). Figure 10b shows the 3D scatter plot and the response surface of Model 2. D2 shows a positive correlation with micropore volume and a negative correlation with average pore diameter, which is consistent with the result mentioned above. Table A1 shows the free hydrocarbon content (S1). S1 of K2qn 1 is ranging from 0.06 mg HC/g to 4.30 mg HC/g, with an average of 1.95 mg HC/g. However, S1 of K2qn 2,3 is in the range of 0.05-1.46 mg HC/g, with an average of 0.44 mg HC/g, which is much smaller than that of K2qn 1 . Therefore, K2qn 1 is much better than K2qn 2,3 in terms of free hydrocarbon content. Figure 11 shows the correlations between fractal dimensions (D1 and D2) and free hydrocarbon content (S1). S1 of lacustrine shale samples is negatively correlated with D1 and D2. The coefficient of determination (r 2 ) between S1 and D1 is more significant than that between S1 and D2, which demonstrates that D1 can be applied to characterize the shale oil content. Shale samples with more free oil usually have a smaller D1, probably because free oil prefers being stored in larger pores, and a simpler pore structure makes oil flow smoother. This is quite distinct from coal and gas shale, which reveal that gas adsorption capacity is positively correlated with fractal dimensions [21]. It is because rougher surfaces (with larger fractal dimensions) can supply more adsorption sites for gas molecules. However, oil molecules, which are much bigger than gas molecules, only can be stored in larger pores (with smaller fractal dimensions). Therefore, considering the effects of pore structure parameters on fractal dimensions, and the correlations between fractal dimensions and oil content (Figure 11), we found that lacustrine shales with larger average pore diameter, smaller BET SSA and fewer micropores will have more oil content, and these shales are better as a shale oil reservoir in terms of pore characteristics. Hence shales in K2qn 1 are better than those in K2qn 2,3 as a shale oil reservoir in the Songliao Basin. Therefore, considering the effects of pore structure parameters on fractal dimensions, and the correlations between fractal dimensions and oil content (Figure 11), we found that lacustrine shales with larger average pore diameter, smaller BET SSA and fewer micropores will have more oil content, and these shales are better as a shale oil reservoir in terms of pore characteristics. Hence shales in K 2 qn 1 are better than those in K 2 qn 2,3 as a shale oil reservoir in the Songliao Basin.

Comparison between Marine Bakken Formation Shales and Lacustrine Qingshankou Formation Shales
The Bakken Formation shale is one of the largest shale oil reservoirs worldwide located in the Williston Basin (North America) [25]. The Bakken Formation shale is a typical marine shale oil reservoir and has been studied extensively [25][26][27]. We chose marine Bakken Formation shales to compare with lacustrine Qingshankou Formation shales to see the difference of pore structure characteristics between them.
The Bakken Formation shales have organic-matter pores, interparticle pores, and intraparticle pores [27]. These kinds of pores also exist in the Qingshankou Formation shales. The isotherms from nitrogen adsorption of Middle Bakken Formation shales are closed and show distinct hysteresis loops, which indicates the existence of the mesopores in the Middle Bakken Formation [25]. However, the isotherms of Upper and Lower Bakken Formation shales are not closed and do not show visible hysteresis loops [25]. All the isotherms from nitrogen adsorption of the Qingshankou Formation shales are closed. They show apparent hysteresis loops, which are the same as Middle Bakken Formation, but different from Upper and Lower Bakken Formation. The Bakken Formation shales have both plate-shaped pores and silt-shaped pores, while the Qingshankou Formation shales have slit-shaped pores and some narrow slit-shaped pores. Figure 12 shows the statistic values (minimum value, average value, and maximum value) of pore structure parameters and fractal dimensions in the Qingshankou Formation and Bakken Formation. Pore structure parameters and fractal dimensions of the Bakken Formation are from Reference [25]. The average pore diameter of the Qingshankou Formation (11.89 nm on average) is a little bit larger than that of the Bakken Formation (10.43 nm on average) (Figure 12a). BET SSA of the Bakken Formation (4.02 m 2 /g on average) is much smaller than that of the Qingshankou Formation (22.70 m 2 /g on average) (Figure 12b). The total pore volume of the Bakken Formation, with an average of 10.30 × 10 −3 cm 3 /g, is much smaller than that of the Qingshankou Formation, with an average of 38.35 × 10 −3 cm 3 /g (Figure 12c). The micro-mesopore volume of the Qingshankou Formation (30.64 × 10 −3 cm 3 /g on average) is much larger than that of the Bakken Formation (9.27 × 10 −3 cm 3 /g on average) (Figure 12d). The pore structure fractal dimension (D 1 ) of the Qingshankou Formation (2.598 on average) is almost the same as D 1 of the Bakken Formation (2.601 on average) (Figure 12e), which demonstrates shales of these two formations have similar complexity of pore structure. However, the pore surface fractal dimension (D 2 ) of the Qingshankou Formation (2.478 on average) is much larger than D 2 of the Bakken Formation (2.135 on average) (Figure 12f), which means the Qingshankou Formation shales have much rougher pore surface than the Bakken Formation shales. The average pore diameter is strongly negatively correlated with D 1 in both Bakken Formation and Qingshankou Formation. However, the correlations between total pore volume and D 1 are opposite in these two formations. The Bakken Formation shows a negative correlation [25], while the Qingshankou Formation shows a positive correlation. This research was performed on the Qingshankou Formation shale, which is considered to be representative of lacustrine shale oil reservoirs in China. The pore characteristics of Qingshankou Formation can also represent more widely of general lacustrine organic-rich shale properties. These methods, including correlation analysis between pore structure parameters and fractal dimensions, and stepwise regression method, can be more extensively applied in other lacustrine organic-rich shales.

Conclusions
In this study, a total of 39 lacustrine shale samples from the Late Cretaceous Qingshankou Formation in the Songliao Basin were examined using the FE-SEM technique and low-pressure nitrogen adsorption. Fractal dimensions were calculated by the FHH model using nitrogen adsorption isotherms. The correlations between fractal dimensions and pore structure parameters were demonstrated using correlation matrix and multidimensional scaling. Prediction of fractal dimensions with pore structure parameters was analyzed using stepwise regression. Based on the results and discussion, the following conclusions can be drawn: This research was performed on the Qingshankou Formation shale, which is considered to be representative of lacustrine shale oil reservoirs in China. The pore characteristics of Qingshankou Formation can also represent more widely of general lacustrine organic-rich shale properties. These methods, including correlation analysis between pore structure parameters and fractal dimensions, and stepwise regression method, can be more extensively applied in other lacustrine organic-rich shales.

Conclusions
In this study, a total of 39 lacustrine shale samples from the Late Cretaceous Qingshankou Formation in the Songliao Basin were examined using the FE-SEM technique and low-pressure nitrogen adsorption. Fractal dimensions were calculated by the FHH model using nitrogen adsorption isotherms. The correlations between fractal dimensions and pore structure parameters were demonstrated using correlation matrix and multidimensional scaling. Prediction of fractal dimensions with pore structure parameters was analyzed using stepwise regression. Based on the results and discussion, the following conclusions can be drawn:

1.
Interparticle pores are highly developed and are the most significant pore type in the Qingshankou Formation. Organic-matter pores are moderately developed and show high heterogeneity. Intraparticle pores and micro-fractures can also be found in the Qingshankou Formation.

2.
The nitrogen adsorption/desorption isotherms of Qingshankou Formation shales are of Type IV, which represent mesoporous nature. According to the slope of isotherms measured by ∆V G /∆(P/P 0 ), most of the samples might have macropores or fractures. The hysteresis patterns are H3 and combination of H3-H4, which demonstrate these shales have slit-shaped pores produced by aggregates of plate-like particles and some narrow slit-shaped pores. Compared with shales in K 2 qn 2,3 , shales in K 2 qn 1 have larger average pore diameter, much smaller BET SSA, and fewer micropores.

3.
Both D 1 and D 2 in K 2 qn 1 have narrower distributions and smaller values than those in K 2 qn 2,3 , which suggests shales in K 2 qn 1 have simpler pore structures and smoother pore surfaces compared to shales in K 2 qn 2,3 .

4.
Using correlation matrix and multidimensional scaling, we can find that BET SSA, micropore volume, D 1 , and D 2 are positively correlated. Average pore diameter shows strongly negative relationships with BET SSA, micropore volume, D 1 and D 2 . Free oil content has negative relationships with D 1 and D 2 . Lacustrine shales with larger average pore diameter, smaller BET SSA, and fewer micropores will have more oil content, and these shales are better as a shale oil reservoir in terms of pore characteristics. Hence shales in K 2 qn 1 are ideal compared to those in K 2 qn 2,3 for a shale oil reservoir in the Songliao Basin. 5.
The lacustrine Qingshankou Formation shales and the marine Bakken Formation shales have similar pore types, but they have different pore shapes. Average pore diameter, BET SSA, total pore volume, and micro-mesopore volume in the Qingshankou Formation are larger than those in the Bakken Formation. These two formations have similar complexity of pore structure, but Qingshankou Formation shales have much rougher pore surface than Bakken Formation shales.