Pore Properties of the Lacustrine Shale in the Upper Part of the Sha-4 Member of the Paleogene Shahejie Formation in the Dongying Depression in East China

Key Laboratory of Petroleum Resources Research, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China Innovation Academy for Earth Science, Chinese Academy of Sciences, Beijing 100029, China University of Chinese Academy of Sciences, Beijing 100049, China Beijing Center for Physical and Chemical Analysis, Beijing 100089, China Roxar Software Solutions, Emerson Process Management, Lysaker 1366, Norway Shengli Oilfield Company of Sinopec, Dongying 257051, China China Institute of Disaster Prevention, Sanhe 065201, China


Introduction
The shale revolution in the USA has changed the global energy landscape with profound impact on international economies [1, 2], which has led to the assessment and explo-ration of shale oil and gas in basins worldwide. It has been conservatively estimated that China has shale oil resources of 40:3 × 10 9 tons, and lacustrine shale holds the vast majority of the shale oil resources in China (40:2 × 10 9 tons) [3]. However, the studies of lacustrine shale oil reservoirs are very limited, compared with many studies of marine shale gas reservoirs [4][5][6][7][8][9][10]. As a result, the pore properties of lacustrine shale and especially their controlling factors are still unclear. The pore properties help determine the petroleum storage and transport properties of shale reservoirs. Therefore, the studies of the pore properties of lacustrine shale and their controlling factors are of great importance.
Compared with conventional reservoirs (e.g., sandstones), shale reservoirs have complex pore structures. Pores in shale are divided into micropores (with diameter < 2 nm), mesopores (with diameter of 2-50 nm), and macropores (with diameter > 50 nm) [11]. Among the techniques for quantitative analysis of shale pore systems [4,[12][13][14], low-pressure gas adsorption (LPGA) is important in the measurement of pore properties over a wide range of pore sizes. Low-pressure N 2 adsorption is frequently used to analyze mesopore characteristics in shale reservoirs [15]. One of reasons is that mesopores are usually much more than micropores and macropores in shale reservoirs [16]. The others are that pore volume (V), pore size distribution (PSD), specific surface area (S), and fractal dimensions can be calculated from a N 2 adsorption isotherm [17][18][19]. However, N 2 adsorption is less accurate in characterizing micropores due to the limitation arising from N 2 molecule and pore throat sizes. A CO 2 molecule has a higher thermal energy at the experimental temperature (273.15 K) and, as a result, can penetrate smaller pores more easily than a N 2 molecule (at 77.35 K) [20]. Low-pressure CO 2 adsorption can be used to analyze micropores. Therefore, both N 2 and CO 2 techniques need to be used to study pore properties in shale.
For shales in the oil window, residual oil in pores dramatically affects the results of pore properties obtained from the N 2 and CO 2 techniques. Removing the residual oil in shale samples via the Soxhlet (solvent) extraction can solve this problem in the analysis of LPGA [21][22][23]. However, published data on shale pore properties in the oil window are rarely obtained from solvent-extracted samples [24][25][26][27][28][29]. Therefore, there is an urgent need to restudy pore properties of shale in the oil window and their controlling factors, based on the solvent extraction.
The geological factors controlling pore properties are usually studied using the bivariate cross-plots and univariate regression analysis [23,25,30,31]. But in most cases, the correlation coefficients of pore properties with geological factors are very low (0.33-0.86) [23,25,30,31]. The reason is that the pore properties are controlled or affected by multiple factors in most cases. Therefore, multivariate regression analysis is necessary. The traditional multivariate regression analysis can develop the multiple regression models but may not work well with the independent variables that are correlated to each other [32,33]. In this case, the traditional method is not applicable to studying the relationships of pore properties with geological factors because the geological factors as independent variables are often correlated to each other. Partial least square regression (PLSR) is a multivariate analytical method that can process the correlated independent variables. PLSR has been used to analyze correlations between pore properties and geological factors, based on seven shale samples [12,27,28]. But the residual oil in samples was not extracted before testing pore properties [12,27,28]. The potential geological factors did not include burial depth and Ro [12,27,28], and the importance of the geological factor (descriptor) for a dependent variable (one of pore properties; response) was not evaluated with the parameter of variable importance in projection (VIP) [12,27,28,33,34]. As a result, the established correlations between the pore properties and the geological factors of the Soxhlet-extracted samples from shale oil reservoirs may be unreliable, and the geological factors controlling pore properties of shale oil reservoirs are still unclear. Therefore, PLSR analysis needs to be used with VIP evaluation when the univariate correlation relationships are poor between pore properties of the Soxhlet-extracted samples and geological factors.
The Dongying Depression is located in the southern part of the Bohai Bay Basin in East China [35]. The lacustrine shale in the upper part of the Sha-4 Member (Es 4 U) of the Paleogene Shahejie Formation displays a great potential for shale oil exploration and development in this depression [22], with total organic carbon (TOC) content of 0.5-11.2 wt% [36], Ro values of 0.5-1.3% [37], thickness up to 400 m [38], and the area of about 3000 km 2 for the shale thicker than 50 m. In this work, the pore properties in fifteen Soxhlet-extracted samples of the Es 4 U shale were analyzed using N 2 and CO 2 adsorption techniques. These samples were also analyzed for geological factors including mineral composition, Ro, and IOC. The relationships between pore properties and geological factors show poor correlations in the bivariate cross-plots and the results of univariate regression analyses. The partial least square regression was then used. The high correlation coefficients (0.72-0.94) show reliability of the PLSR analyses. The importance of the geological factors (descriptors) for the pore properties (responses) was evaluated with VIP scores. The results illustrate that clay minerals and carbonates with high VIP scores are two key factors that affect pore properties of the lacustrine shale. The reason is that many pores in shale samples are developed between clay minerals but some are filled with carbonates. As the shale with higher volumes of micro-and mesopores contains more clay minerals and fewer carbonates, this shale is unfavorable for production via hydraulic fracturing. The analysis of the shale pore evolution helped us find the shale with the modest micromesopore volume, brittle mineral, and clay mineral contents which is conducive to the development of the Es 4 U shale oil.

Geological Setting
The Dongying Depression of the Bohai Bay Basin in eastern China is a lacustrine basin with a faulted margin in the north and a ramp margin in the south (see Figure 1 in Liu et al. [39] and Figure 1). The Paleogene and Neogene successions in the depression consist of clastics with subordinate lacustrine carbonates and evaporites, and they mainly comprise four formations from bottom to top: the Paleogene Shahejie (Es) and Dongying (Ed) Formations and the Neogene Guantao (Ng) and Minghuazhen (Nm) Formations (see the strata column ( Figure 2) in Zhang et al. [35]). The Shahejie Formation is further divided into the Sha-4 (Es 4 ), -3 (Es 3 ), -2 (Es 2 ), 2 Geofluids and -1 (Es 1 ) Members from bottom to top. The shale in the upper part of the Sha-4 Member (Es 4 U) is an important source rock in the Dongying Depression [22]. The thickness of Es 4 U shale is up to 400 m [38]. The area with thickness greater than 50 m (154 m on average) occupies about 3000 km 2 and is mainly located in the central-north part of the Dongying Depression ( Figure 1). The thick Es 4 U shale is good source rock, mainly containing type I kerogen and followed by type II kerogen, and has TOC content of 0.5-11.2%, Ro values of 0.5-1.3%, and hydrocarbon generation potential (free and cracked hydrocarbon, i.e., S 1 + S 2 ) values of 0.62-76.51 mg/g [29,[36][37][38]40]. The average oil saturation of Es 4 U shale within the oil window (samples with depths greater than 3000 m in Zhang et al. [41]) is 45.8%. The paleosalinity was high during the Es 4 U deposition [42], facilitating the development of carbonates (see the strata column in Zhang et al. [35]).

Samples and Methods
3.1. Samples. To reveal the pore properties of the main bodies of the shale, fifteen core samples of the Es 4 U shale were selected from nine wells (Figures 1 and S1), based on the sample locations in Es 4 U, and observation of rock textures and carbonate content in specimens (Table S1). These samples were divided into two types according to the rock textures: laminated shale (13 samples) and massive mudstone (2 samples; see Table S1). The laminated shale is dark, gray, and brown in color and includes horizontal and low-angle parallel laminations. They are mainly composed of interlaminated and interbedded mudstone and carbonates. The carbonates are mainly calcite (Figure 2(a)). In addition, several shale samples contain silty sands (Figure 2(b)). Massive mudstone is dark to gray. Its structure is shown in Figure 2(c).

Soxhlet
Extraction. Approximately 20-40 g of each sample was crushed to the powder finer than 60 mesh [23,44] and processed with Soxhlet extractors and 300 ml [44] of a dichloromethane/methanol mixture (DCM/MeOH, 93 : 7 vol/vol) for 72 h [23,44]. The water bath temperature of the extraction was 48°C. The complete extraction of residual oil was guaranteed by the observation that solvent in the siphon and thimble-holder becomes colorless [45]. The Soxhletextracted samples were placed in beakers and dried under a fume hood. Each dried sample was split into four aliquots. Three of them were grounded into the powder finer than 100 mesh.

3.3.
Microscopy. The particles of the original (Soxhlet-unextracted) and Soxhlet-extracted shale samples were embedded in Araldite 502 resin and polished into the fluorescent thin sections. Whether there is residual oil in the shale samples was checked by fluorescence observation of thin section through an Olympus BX-51 optical microscope equipped with a mercury lamp (USH102D), an excitation filter (UMNU2, 360-370 nm), and a longpass emission filter (LP400, >400 nm) [46]. Thin sections were also polished with argon ions by using a Leica EM TIC 3X with an acceleration voltage of 5 kV and a    4 Geofluids current of 2 mA for 2 h to create an advanced smooth surface for electron microscope observation [47]. After coating with gold at a thickness of 2 nm on the surface, the thin sections were observed by an scanning electron microscope (SEM, Zeiss Crossbeam 540) with field emission gun [47]. The acceleration voltage was 1.2 kV.
3.4. X-Ray Diffraction, Insoluble Organic Carbon Content, and Vitrinite Reflectance Methods. One ground aliquot was tested for mineral composition by using a Rigaku D/max-2500PC X-ray diffractometer. The method and analytical conditions are documented in Zhang et al. [48]. Another ground aliquot of each sample was tested for insoluble organic carbon (IOC) [49]. IOC is the total organic carbon content of the samples after removal of soluble organic matter by Soxhlet extraction with DCM/MeOH. The Soxhletextracted samples finer than 100 mesh were analyzed for IOC with a LECO CS230 carbon/sulfur analyzer [44]. The kerogen in the other ground aliquot of each sample, separated from the shale sample and polished into a thin section, was tested using a J&M TIDAS MSP 400 microspectrophotometer for vitrinite reflectance (Ro) [50].
3.5. Low-Pressure Gas Adsorption. Low-pressure N 2 and CO 2 adsorption experiments were conducted in turn using a Micromeritics ASAP 2020 surface area analyzer. One of the Soxhlet-extracted aliquots was sieved to obtain particles between 60 and 80 mesh. Four grams of each sieved sample was degassed at 383.15 K in vacuum (pressure < 0:01 torr) for 20 h to remove residual volatile material before the measurements. The N 2 and CO 2 adsorption isotherms were measured at 77.35 K and 273.15 K, respectively. The tested ranges of relative equilibrium adsorption pressure (the ratio of adsorption equilibrium pressure to saturated vapor pressure, denoted by P/P 0 ) in N 2 and CO 2 adsorption experiments were set as 0.005-0.980 and 0.0005-0.0300, respectively.
Experiments were regarded to reach equilibration when the ratio of the pressure variation to the average pressure per time interval (60 s) was less than 0.01%. Using the density functional theory (DFT) [17,51], the pore volumes (V) and pore size distributions (PSD) were analyzed from the N 2 and CO 2 adsorption data. Based on the international standard ISO 9277:2010 [52], the specific surface areas (S, ratio of the total surface of a sample to its mass) were calculated with the Brunauer-Emmett-Teller (BET) model [18] from the N 2 adsorption data under the relative equilibrium pressures ranging from 0.05 to 0.30. The average pore diameter was customarily calculated under the assumption that all pores are cylindrically shaped, equivalent to quadrupling the quotient of the N 2 adsorption pore volume and specific surface area [53]. The specific surface areas tested by CO 2 adsorption were not reported in this paper, due to large uncertainties in the area occupied by an adsorbed CO 2 molecule in the pores of molecular dimensions and the CO 2 -measured monolayer capacity (the amount of gas molecules of a complete, closepacked monolayer covered on the solid surface) [24,54].
3.6. Fractal Analysis. Fractal geometry has been applied to the irregularity of porous media [55]. In this study, the fractal FHH (Frenkel-Halsey-Hill) model [19] was utilized to analyze the fractal dimensions from the N 2 adsorption data. Based on the FHH model, there are two methods for calculation of fractal dimensions, using the van der Waals force regime or capillary condensation regime [56]. The method based on the capillary condensation regime is used in this paper, as it is suitable for fractal analysis of porous medium [57]. It can be described as follows:

Geofluids
where K is the slope of the line of ln V vs. ln ½ln ðP 0 /PÞ, V is the adsorbed nitrogen volume (cm 3 /g) at the equilibrium pressure P (Pa), P 0 (Pa) is the saturated vapor pressure of nitrogen at 77.35 K; C is a constant (dimensionless), and D (dimensionless) is the fractal dimension. The details are described in Yao et al. [57]. If a ln V vs. ln ½ln ðP 0 /PÞ line shows two segments with different slopes, two kinds of fractal dimensions (D 1 and D 2 ) need to be calculated also using Equations (1) and (2). The fractal dimension D 1 , calculated with the adsorption data at P/P 0 < 0:45, reflects the complexity of pore surfaces, and D 2 , calculated with the adsorption data at P/P 0 > 0:45, represents the complexity of pore space structures [58].
3.7. Partial Least Square Regression. Partial least square regression (PLSR) was employed for establishing the relationships of pore property parameters with geological factors in this study. The pore properties, including micropore volume (V mic ), mesopore volume (V mes ), micromesopore volume (V mic-mes , sum of V mic and V mes ), specific surface area (S), average pore diameter (APD), and fractal dimensions (D 1 and D 2 ), with the appropriate forms (e.g., ln V mes and ln S; see Section 4.3 for details) used as dependent variables (responses) in the PLSR analyses, and the geological factors (mineralogy, IOC, burial depth, and lnRo) as independent variables (descriptors). Each variable (descriptor or response) was firstly standardized through subtracting its mean and dividing by its standard deviation. Using the PLSR analysis, orthogonal (uncorrelated to each other) components and linear combinations of descriptors are first set up: where F i is the ith component (i = 1, 2, ⋯, m), x j is the jth standardized descriptor (j = 1, 2, ⋯, n), and a j is a coefficient. Then, the regression linear equation of a standardized response versus one or several orthogonal components is established: where y is the standardized response, b j is a coefficient, and C 1 is a residual. The linear equation of a response versus descriptors is obtained from Equations (3) and (4): where c j is a coefficient and C 2 is a residual. The components in PLSR are constructed through maximizing the information of the descriptors to efficiently predict the standardized responses [59]. The number of the components was determined by cross-validation (CV) [60]. The predicted value of the fitting model of an unstandardized response is equal to the mean plus the result of multiplying the predicted value of a standardized response and standard deviation (i.e., inverse standardization). The means and standard deviations used for the inverse standardization are equal to those that were used for the standardization of the responses. In PLSR analyses, the variable importance in projection (VIP) of a descriptor to a response is used to quantitatively describe the importance of a descriptor for a response and is computed from the weights of descriptors in components and correlation coefficients of descriptors with components [33,61]. Larger VIP indicates more importance of a descriptor for a response. The descriptors with VIP values greater than one (the average of square VIP values) show the above-average importance for the prediction of the response and are considered to be the relatively important descriptors [33,34].  Table S1. The mineral composition of massive mudstone samples is not evidently different from that of laminated shale samples. The samples from the top, middle, and bottom of the Es 4 U shale do not show evident difference in mineral composition either. These samples are not statistically analyzed separately according to their texture types and locations in Es 4 U. Clay minerals range from 6 wt% (sample DYS14) to 53 wt% (sample DYS12) with an average of 30.1 wt%. Quartz content varies from 3 wt% (sample DYS14) to 48 wt% (sample DYS13), and the average is 23.6 wt%. Feldspar content is between 4 wt% (sample DYS14) and 37 wt% (sample DYS3), and its average is 18.0 wt%. Calcite content varies widely from 3 wt% (sample DYS7) to 85 wt% (sample DYS14) with an average of 15.7 wt% and dolomite from 0 wt% (sample DYS2) to 54 wt% (sample DYS4) (10.4 wt% on average). The abundant carbonates arose from the high paleosalinity during the Es 4 U deposition [42].

Results and Discussion
Carbonates are negatively related to clay minerals and quartz (Figures S2(a) and S2(b)). Comparison of brittle mineral (including carbonates, quartz, and feldspar) content of shales in different areas is helpful for shale evaluation. In this study, we compared the brittle mineral content of Es 4 U shale in the Dongying Depression with those of Chang-7 shale in the Ordos Basin, from which commercial flows of oil have already been produced using horizontal drilling and hydraulic fracturing techniques. The content of brittle minerals in the Es 4 U shale is 68.3 wt% on average (from 47 wt% to 94 wt%), higher than that of the lacustrine Chang-7 (49.1 wt%) [23]. Relatively high brittle mineral content in the Es 4 U shale is favorable for hydraulic fracturing.
The other shales were studied with Soxhlet-unextracted samples (see Section 4.2.2 for details), such as the Bakken shale in the Williston Basin and Kong-2 shale in the Cangdong Depression [25][26][27][28]. The results about the pores in these studies cannot be compared with the pores in the Soxhlet-extracted samples in our study. We did not compare the brittle mineral content in the Soxhlet-unextracted 6 Geofluids samples with those in the Soxhlet-extracted samples in our study, as this comparison is not helpful to study the pore properties of the Es 4 U shale.
The IOC content of the Es 4 U shale ranges from 1.14 wt% (sample DYS14) to 5.52 wt% (sample DYS8), with an average of 2.30 wt% (Table S1). The IOC content of the Es 4 U shale is lower than that of the Chang-7 shale (average IOC = 6:93wt% [23]) but higher than the TOC values of the lacustrine Qianfoya shale (on average 1.21 wt%) from which commercial flows of oil have also already been produced using horizontal drilling and hydraulic fracturing techniques [62]. The relatively abundant organic matter in the Es 4 U shale in the Dongying Depression is a result of high paleoproductivity, warm humid climate, and deep-water environment during deposition [38].
The measured Ro values are listed in Table S1. Sample DYS13 was not tested due to insufficient vitrinite. Its Ro value was calculated by using the regression equation of the measured Ro values versus the burial depth ( Figure S3). These measured and calculated Ro values vary from 0.60% (sample DYS14) to 1.29% (sample DYS7) (Table S1), with an average of 0.83%. This is very similar to Chang-7 shale whose Ro values vary from 0.64% to 1.34%, with an average of 0.88% [23]. Most Es 4 U shale samples are in the oil window (Ro = 0:6 -1:0%) and some in the gas condensatewet gas window (Ro = 1:0 -1:3%) according to U.S. EIA [63] and Dembicki [64]. . Therefore, most of the residual oil in our samples was removed. Then, all of fifteen Soxhlet-extracted samples were analyzed using low pressure CO 2 and N 2 adsorption. 4.2.1. CO 2 and N 2 Adsorption Isotherms. The total adsorbed CO 2 volumes in the CO 2 adsorption experiments range from 0.15 cm 3 /g to 2.07 cm 3 /g at standard temperature and pressure (STP). The CO 2 adsorption isotherms of all samples are convex in shape when plotted against the relative pressure (P/P 0 ) ( Figure S4) and belong to type I, which displays the features of micropores, according to the IUPAC (International Union of Pure and Applied Chemistry) classification [11]. In the N 2 adsorption experiments, the total adsorbed N 2 volumes vary from 2.12 cm 3 /g to 29.61 cm 3 /g at STP. The N 2 adsorption isotherms of all shale samples belong to type IV with hysteresis loops (Figure 3), showing the features of mesopores [11]. These hysteresis loops display rapid increases of desorption branches at relative pressures of about 0.45. These features of the CO 2 and N 2 adsorption isotherms indicate that the samples from the Es 4 U all contain micro-and mesopores.
To further check the effect of the residual oil removal, three original samples (DYS3, DYS10, and DYS11) were also analyzed. Both total adsorbed CO 2 and N 2 volumes were evidently increased after the Soxhlet extraction (Figures 3 and  S4). This indicates that the pore space occupied by residual oil is indeed released.
The N 2 hysteresis loops are categorized into types H1, H2, H3, and H4 according to the IUPAC classification [11]. For the Soxhlet-extracted samples from the Es 4 U shale, all hysteresis loops are not similar to that of type H1, but those loops of five samples (33.3%, 5/15) are similar to type H2 (Table S2) whose desorption branch never plots parallel to the adsorption branch, with an inflexion at a relative pressure of 0.5, as shown in Figure 3(a). The pores in samples with the hysteresis loops of type H2 are mainly inkbottle-shaped (narrow necks and wide bodies) [11] and are usually intergranular, intercrystalline, and organic pores [24,65]. Figure 2(g) shows the intercrystalline pores of pyrite, an example of an inkbottle-shaped pore. Type H3 loop is characterized by parallel adsorption-desorption branches above medium relative pressure (~0.5) and steeply increasing isotherms near the saturated vapor pressure (i.e., P/P 0 = 1) (Figure 3(b)). The pores in samples with the hysteresis loops of the type H3 loop are mainly slit-shaped [11]. The hysteresis loops of eight samples (53.3%, 8/15) are kin to type H3 (Table S2). The type H4 loop is similar to type H3 at low-medium relative pressure, but the isotherms increase slowly at high relative pressure (see Figure 3(c)). The pores in samples with the type H4 loop are mainly narrow slit-shaped [11]. The hysteresis loops of two samples (13.3%, 2/15) are similar to that of type H4 (Table S2). The (narrow) slit-shaped pores are usually the pores between plate-like particles (e.g., clay sheets) and fractures [24,65]. Figures 2(h) and 2(i) show the fractures and pores between clay minerals, respectively. Fractures usually are wider than narrow slit-shaped pores between plate-like particles [40]. The samples mainly containing slit-shaped pores are significantly more than the narrow slit-shaped pore-dominated samples, probably because there develop massive fractures in the Es 4 U shales. The hydrocarbon flow capability of different types of pores is discussed together with the fractal characteristics in Section 4.2.3.

Pore Parameters and
Pore-Size Distribution. The pore parameters discussed here include micropore volume, mesopore volume, specific surface area, and average pore diameter. Mesopore volume, specific surface area, and average pore diameter were calculated using the methods introduced in Section 3.5 and data from N 2 and CO 2 adsorption experiments. The results are listed in Table S2.
The micro-and mesopore volumes and specific surface area for the three pairs of original and Soxhlet-extracted samples are listed in Table S2. Figure 4 shows the distribution curves of micropores (0.3-2 nm in diameter) and mesopores (2-50 nm in diameter) in the original and Soxhlet-extracted samples in Es 4 U. The amount of the pores with the diameter less than 10 nm increase significantly. Therefore, the average pore diameter decreases after the Soxhlet extraction (Table S2).

Geofluids
For the Soxhlet-extracted shale samples in Es 4 U, the micropore volumes (Table S2) vary from 0.0006 cm 3 /g (sample DYS14) to 0.0071 cm 3 /g (sample DYS6) (0.0038 cm 3 /g on average). The mesopore volumes range from 0.0028 cm 3 /g (sample DYS14) to 0.0381 cm 3 /g (sample DYS15) with an average of 0.0181 cm 3 /g. The mesopore volume is much higher than the micropore volume. The micromesopore volume (sum of micropore and mesopore volumes) varies from 0.003 cm 3 /g to 0.045 cm 3 /g with an average of 0.022 cm 3 /g. The specific surface areas, tested by N 2 adsorption, range from 1.04 m 2 /g (sample DYS14) to 30.55 m 2 /g (sample DYS15) with an average of 12.56 m 2 /g. The average pore diameters of the samples, obtained from the N 2 adsorption data, are between 4.97 nm (sample DYS1) and 12.46 nm (sample DYS13) with an average of 7.25 nm. The average pore diameter and specific surface areas were not calculated using the CO 2 adsorption data, due to their large uncertainties (see subsection 3.5).
To understand the pores in the Es 4 U shale, the pore parameters of the Es 4 U shale need to be compared with those of other shales. But in most of previous studies, the shale samples were analyzed without solvent extraction [24][25][26][27][28]. As a result, the pore data from these previous studies whose samples were not extracted cannot be compared with data in this study, due to the pore parameter values before and after Soxhlet extraction change (e.g., samples DYS3, DYS10, and DYS11 in Table S2). However, the samples from the lacustrine Chang-7 shale in the Erdos Basin were analyzed using low-pressure N 2 and CO 2 adsorption methods, after the solvent extraction [23,44]. The Chang-7 shale has the average micropore volume of 0.0035 cm 3 /g, mesopore volume of 0.0113 cm 3 /g, specific surface area of 4.99 m 2 /g, and average pore diameter of 17.54 nm [23,44]. Compared with the Soxhlet-extracted samples from the Chang-7 shale, the Es 4 U shale has significantly higher micro-and mesopore volumes (0.0038 cm 3 /g and 0.0181 cm 3 /g on average, respectively), but the average pore diameter is lower. The average oil saturation of mature Es 4 U shale (45.8%) [41] is close to that in the Chang-7 shale (44.63%) [66]. The Es 4 U shale is much thicker (see Section 3.1) and has more brittle minerals (see Section 4.1) than the Chang-7 shale (with thickness mostly lower than 100 m [67]). The roof and floor of the Chang-7 shale are tight sandstone [68], whereas the roof and floor of the Es 4 U shale are mudstone, shale, and gypsum [38,69]. All the characteristics of shale reservoirs except the average pore diameter display that the Dongying Es 4 U shale is favorable for the development of shale oil and may achieve higher oil flow than the commercial oil flow   [19] applied to the N 2 adsorption data (see Section 3.6 for details). The ln V vs. ln ½ln ðP 0 /PÞ lines of all samples show two segments divided at the relative pressure of 0.45 (as displayed in Figure 5(a)). The correlation coefficients (R) of these segments with different slopes are all higher than 0.97 (Table S2), indicating that the fractal characteristics are evident. From these two segments, two kinds of fractal dimensions (D 1 and D 2 ) were calculated for each sample. The values of the fractal dimension D 1 , calculated with the adsorption data at P/P 0 < 0:45 (region 1 in Figure 5(a)), range from 2.25 (sample DYS4) to 2.50 (sample DYS1) (2.42 on average). In region 1, the monolayer adsorption is dominant, and thus, the D 1 reflects the complexity of pore surfaces [58]. The D 2 values, calculated with the adsorption data at P/P 0 > 0:45 (region 2 in Figure 5(a)), are from 2.54 (sample DYS14) to 2.80 (sample DYS1) with an average of 2.71. In region 2, nitrogen molecules fill the pore space with multilayer adsorption, and thus, the D 2 reflects the complexity of pore space structures [58]. Higher D 1 and D 2 indicate that the samples have rougher pore surfaces and more complicated pore space structures, respectively [24]. Therefore, pores with low D 1 and D 2 values are favorable for oil flow. Figure 5(b) illustrates the average values of the fractal dimensions of the pores in the samples with the hysteresis loops of types H2, H3, and H4 that mainly contain inkbottle-, slit-, and narrow slit-shaped pores, respectively [11]. Slit-shaped pores with better openness have less curvature internal surface and simpler network of pores than inkbottle-shaped and narrow slit-shaped pores and thus have a higher hydrocarbon flow capacity [70]. In the Es 4 U shale, the slit-shaped pore-dominated samples are the most (53.3%) and are favorable for hydrocarbon flow. The shale samples mainly containing slit-shaped pores have lower D 1 and D 2 values than those mainly containing inkbottle-and narrow slit-shaped pores ( Figure 5(b)). One reason is that inkbottle-shaped pores have more complex pore surfaces and space structures than slit-shaped pores. The other reason is that although slit-shaped pores and narrow slit-shaped pores are similar in shape, the aggregation of small pores has more complicated pore surfaces and space structures than that of the relatively large pores. Therefore, inkbottleand narrow slit-shaped pores (types H2 and H4) have higher values of fractal dimensions than slit-shaped pores and display certain similarities in fractal dimensions. These illustrate that the classification of the hysteresis loop types corresponds to the values of fractal dimensions and thus is reliable.  Figures 6(a) and S5(a) show that the average pore diameter is negatively correlated with both mesopore volume and specific surface area. Therefore, the quantity increase of relatively small pores in the shale is an important factor leading to the increases of mesopore volume and specific surface area in Figure 6(b). Furthermore, the samples with the hysteresis loops of types H2, H3, and H4 show the similar correlation relationships between the mesopore volume and the specific surface area (Figure 6(b)), suggesting that the influence of the geometrical shapes of the pores on the mesopore volume and specific surface area is relatively weak and masked by the strong influence of the pore quantity. Figure 6(c) shows that the average pore diameter has a significantly negative correlation with fractal dimension D 2 (R = 0:98). The reason is that the aggregation of small pores has more complicated pore space structures than that of the relatively large pores. Moreover, small pores have more curvature internal surface. Therefore, the average pore diameter is also negatively related to the D 1 (Figure 6(d)), which is also supported by the positive correlation between D 1 and D 2 ( Figure S5(b)). Because of the correlations between the average pore diameter, specific surface area, and mesopore volume (Figures 6(a), 6(b), and S5(a)), D 1 and D 2 are also positively correlated to the specific surface area and mesopore volume (Figures S5(c)-S5(f)). As mentioned above, the samples with the relatively high mesopore volume and specific surface area usually have higher percentage of relatively small pores and thus have higher D 1 and D 2 .

Correlations between Pore Parameters.
In summary, the increase of small pores in the shale leads to the increases of the pore volumes, specific surface area, and fractal dimensions (D 1 and D 2 ) and the decrease of the average pore diameter. But the correlations of D 1 with average pore diameter (R = 0:52, Figure 6(d)) and D 2 (R = 0:52, Figure S5(b)) are not very evident. This may suggest that the D 1 of the Es 4 U shale is also affected by other factors, such as mineral composition, burial depth, and Ro (see below).

Geological Factors Controlling Shale Pores.
We analyzed the relationships of geological factors with pore properties using cross-plots and univariate regression analyses. The pore properties are pore volume, diameter, specific surface area, and fractal dimensions. The geological factors (i.e., variables) include insoluble organic carbon (IOC), burial depth, vitrinite reflectance (Ro), and content of clay minerals, quartz, feldspar, carbonates (sum of calcite, dolomite, and siderite), calcite, and dolomite (Table S1). There are correlations between the pore parameters and geological factors, as shown in Figures 6(e)-6(h) and S5(g)-S5(l), but most of the correlation coefficients are low. Similar weak correlations have been found in many other shales [23,25,31].
In fact, pore formation is usually controlled or affected by multiple geological factors. Therefore, PLSR analyses (see Section 3.7 for details) need to be conducted to reveal multivariate correlations. Before the PLSR analyses, some pore parameters and geological factors were transformed into a logarithmic form in this work. Univariate regression analyses show that the logarithms of the pore volumes are correlated with the average pore diameter (Figure 6(a)) and fractal dimensions D 1 (Figure S5(d)) and D 2 ( Figure S5(f)), as well as the clay minerals and carbonates (Figures 6(e), 6(f), S5(g), and S5(h)). Therefore, the pore volumes were transformed into the logarithmic form. The mesopore volume has remarkable positive correlation with the specific surface area (Figure 6(b)). The logarithm of the specific surface area was taken as a dependent variable. Considering that shale porosity (total pore volume × apparent density) has an exponential relationship with the burial depth under equilibrium compaction [71] and that the log of Ro increases linearly with burial depth (Figure S3), we transformed Ro into the logarithmic form.    Table S2. 10 Geofluids

Geofluids
(independent variables) include burial depth, mineral contents, IOC, and lnRo. The results of the PLSR analysis, including standardized coefficients (coefficients in regression equations obtained from standardized data) and VIP (variable importance in projection, see Section 3.7) values of the descriptors, are listed in Table S3. The predicted values of the PLSR fitting model of ln V mic were computed through the inverse standardization (see Section 3.7) of the regression values obtained from the standardized coefficients of the multivariate regression equation in Table S3 and the standardized data. These standardized data were calculated from the data in Tables S1 and S2. The results are listed in Table S4. The relationship between the experimental and predicted values of the V mic have a much higher correlation coefficient (R = 0:91, Figure 7(a)) than those of univariate regression analyses (R = 0:69 -0:70, Figure 6(e) and S5(g)).
The standardized coefficients of the multivariate regression equation and VIP values (Figure 8(a) and Table S3) of the PLSR fitting model of ln V mic indicate that the micropore volume is mainly correlated positively to clay minerals, quartz, and IOC but negatively to the variables of carbonates and calcite, with the VIP values higher than 1 (see Section 3.7 for details). Pyrite, feldspar, dolomite, burial depth, and lnRo are subordinate variables in the multivariate regression equation with the VIP values lower than 1 (also see Section 3.7 for details). Micropores can be developed between clay minerals [72] and are affected by diagenesis although clay minerals are mainly terrigenous. From the early middle diagenetic stage, smectite illitization can increase intercrystalline pores [30,73]. Besides, some pores along cleavage planes of clay minerals can be produced by the mineral distortion due to compaction [72]. Carbonates are unfavorable to pore development when carbonates fill pores [74]. At the middle diagenetic stage, the acid fluid produced from the organic-rich shale in the Es 4 U shale can dissolve some carbonates, but subsequently, the cementation of carbonates occurred when the solution became saturated again due to shale usually is a relatively closed system [30,40]. Quartz grains probably play a supporting role in compaction or are negatively related to carbonates ( Figure S2(b)). The positive correlation of micropore volume with IOC supports that some micropores can be produced from organic matters (mainly kerogen) with thermal maturation [31]. The eight Soxhlet-extracted samples from the Chang-7 lacustrine shale in the Ordos Basin also show a positive relation of the micropore volume with the IOC [23]. This correlation with IOC illustrates that organic pores are an important part of micropores. But, it is worth noting that for the Es 4 U shale in the Dongying Depression, the VIP value of IOC (1.14) is evidently lower than those of the carbonates (1.50) and clay minerals (1. 46), which indicate that the inorganic factors are more important for controlling micropores than IOC in the Es 4 U shale.

Geological Factors Controlling Mesopore Volume and Specific Surface
Area. Table S3 also shows the standardized coefficients and VIP values of the descriptors obtained from the PLSR analysis of the log of the mesopore volume (ln V mes ). The predicted values of the PLSR fitting model of ln V mes (Table S4) were calculated through the inverse standardization (see Section 3.7) of the regression values obtained from the standardized data and standardized coefficients of the multivariate regression equation in Table S3. The relationship between the experimental and predicted values of V mes (Figure 7(b)) has a much higher correlation coefficient (R = 0:94) than those of univariate regression analyses (R = 0:72 -0:85, Figures 6(f) and S5(h)).
The VIP values and standardized coefficients of the  Figure 6: Correlation relationships between pore property parameters and between geological factors and pore property parameters for the Soxhlet-extracted samples from the Es 4 U shale: (a-d) between pore parameters (mesopore volume, specific surface area, average pore diameter, and fractal dimensions D 1 and D 2 ); (e-h) between pore parameters (micro-and mesopore volumes, average pore diameter, and fractal dimension D 1 ) and mineral content (carbonates and clay minerals). The samples with the hysteresis loops of types H2, H3, and H4 mainly contain inkbottle-, slit-, and narrow slit-shaped pores, respectively. The data are listed in Tables S1 and S2. multivariate regression equation (Figure 8(a) and Table S3) of the PLSR fitting model of ln V mes indicate that the mesopore volume is mainly correlated positively to clay minerals and negatively to carbonates and calcite, with the VIP values higher than 1 (see Section 3.7 for details). This correlation with clay minerals and carbonates is similar to the model of ln V mic , but IOC is only a subordinate factor for mesopore volume (Figure 8(a)). The uncorrelated relationship with IOC coincides with the fact that the mesopores of the Es 4 U samples with the hysteresis loops of types H3 and H4 are dominant (Figure 7(b)), which are mainly inorganic pores (see Section 4.2.1). The difference in the relationships of micropore and mesopore volumes with IOC suggests that organic pores are very subordinate in mesopores but are important in micropores. As mentioned above, carbonates are unfavorable to pore development when they fill pores. Pores between clay minerals can be increased by diagenesis, such as smectite illitization and distortion of clay minerals due to compaction [30,72,73]. Therefore, the mesopore volume is mainly correlated to carbonates and clay minerals. Quartz, pyrite, burial depth, feldspar, dolomite, and lnRo are all subordinate variables (factors) in the multivariate regression equation with the VIP values lower than 1 (see Section 3.7 for details), except for IOC.
In the Sichuan Basin and Williston Basin, the mesopore volumes of the Longmaxi (thirteen samples) and Bakken (eight samples) marine shales are all related to TOC content or Ro [25,31]. These correlations are different from that for the Es 4 U lacustrine shale in the Dongying Depression. One of the reasons probably is that terrigenous minerals in the lacustrine shale are more abundant than those in the marine   Figure 7: Correlation relationships between experimental and predicted pore property parameters (V mic , V mes , APD, and D 1 ) in the PLSR fitting models. V mic : micropore volume; V mes : mesopore volume; APD: average pore diameter; D 1 : fractal dimension D 1 . The samples with the hysteresis loops of types H2, H3, and H4 mainly contain inkbottle-, slit-and narrow slit-shaped pores, respectively. The data are listed in Tables S2 and S4. 13 Geofluids shales so that the mineral composition of the lacustrine shale becomes more important for controlling pores than organic matter.
The log-transformed specific surface area (ln S) and micromesopore volume (ln V mic-mes ) were also analyzed with PLSR. The specific surface area and micromesopore volume all have significantly positive correlations with the mesopore volume ( Figure 6(b) and S5(m)), due to the mesopore volume being related to the specific surface area (see Section 4.2.4 for details) and being much higher than the micropore volume (Table S2). The standardized coefficients of the multivariate regression equations in the PLSR fitting models of ln V mic-mes and ln S are close to those in the ln V mes model (Table S3). The clay minerals, carbonates, and calcite also have VIP values higher than 1 (Figures 8(b) and 8(c) and Table S3) and are the relatively important factors for characterizing ln V mic-mes and ln S (see Section 3.7 for details). These features are the same as those in the PLSR fitting model of ln V mes .
In summary, the PLSR fitting models have higher correlation coefficients than those of univariate regression for analyzing the correlations of pore volumes with the multiple geological variables (factors), supporting that the pore vol-umes are controlled or affected by multiple geological factors. The mesopore and micromesopore volumes as well as the specific surface area of the Dongying Es 4 U shale are mainly controlled by the mineral composition that are determined by sedimentation and diagenesis. The micropore volume is mainly correlated positively to clay minerals, quartz, and IOC but negatively to the variables of carbonates and calcite. The difference in the relationships of micropore and mesopore volumes with IOC suggests that organic pores are very subordinate in mesopores but are important in micropores. The shale with high clay mineral and low carbonate content has more micromesopores for hydrocarbon storage but may be unfavorable for hydraulic fracturing. In contrast, although the pore volumes in the shale with abundant carbonates are low, the high content of carbonates facilitates fracture forming. Therefore, pore volumes and brittle minerals in the Es 4 U shale need to be synthetically taken into account in the "sweet spot" selection and determination.     Table S3. 14 Geofluids standardized coefficients and VIP values of geological factors in the PLSR fitting model of the APD are listed in Table S3. The relationship between the experimental and predicted values of APD, computed through the inverse standardization of the regression values of the standardized data, have a higher correlation coefficient (R = 0:76, Figure 7(c)) than those of univariate regression analyses (R = 0:42 -0:53, Figures 6(g) and S5(i)-S5(k)). APD is mainly correlated negatively to the clay minerals, feldspar, burial depth, and lnRo and positively to the carbonates and calcite, based on the VIP values higher than 1 (Figure 8(c) and Table S3). The reason for the negative correlations of APD with clay minerals and feldspar is probably that mesopores between crystal palettes of clays and cleavages in feldspar grains are usually small. The average pore diameter decreases with the increasing burial depth, due to the compaction. As the burial depth is correlated to lnRo ( Figure S3), APD is also correlated negatively to lnRo.
The standardized coefficients and VIP values higher than 1 in the PLSR fitting model of D 2 (Figure 8(d) and Table S3) indicate that clay minerals, burial depth, carbonates, calcite, and feldspar are relatively important factors. lnRo has a VIP value close to 1. If lnRo is also considered as a relatively important variable, the main geological factors controlling D 2 are similar to those controlling APD. But the effects of these geological factors on the response (D 2 or APD) are the opposite, as APD has a negative correlation with D 2 (Figure 6(c)).  Table S3) have relatively important influences on D 1 (see Section 3.7 for details). As mentioned above, the developed clay minerals lead to the increase of surface area and surface roughness of pores. But the dolomite in the Es 4 U shale is mainly self-structured crystals with the flat surface [75]. Therefore, D 1 is mainly correlated positively to clay minerals and negatively to dolomite.
In summary, average pore diameter and fractal dimensions are also controlled by multiple geological factors and their PLSR fitting models have higher correlation coefficients than those of univariate regression. The average pore diameter and fractal dimension D 2 of the Es 4 U shale are mainly controlled by multiple geological factors including clay minerals, feldspar, burial depth, lnRo, carbonates, and calcite, while the fractal dimension D 1 is mainly correlated positively to clay minerals and negatively to dolomite. The most important factor is the content of clay minerals that affect both pore size and fractal dimensions.
4.4. Pore Evolution. As inorganic pores are dominant in the Es 4 U solvent-extracted samples, absolute pore volumes were directly used to analyze the pore evolution, instead of TOCor IOC-normalized pore volumes [76]. The cross-plots of micropore and mesopore volumes, specific surface area, average pore diameter, and fractal dimensions with Ro and burial depth are shown in Figure 9. The burial depth was calculated with the regression equation in Figure S3, so that the uniform cross-plots in Figure 9 can be drawn by using both Ro and burial depth.
Over the range of 0.6-1.3% Ro, Figures 9(a)-9(f) all display scatter features. However, the envelope curves of micropore and mesopore volumes as well as specific surface area and average pore diameter show a bimodal distribution with two peaks at 0.7% and 0.9% Ro (Figures 9(a)-9(d)). The bimodal distribution of envelope curves of fractal dimensions D 1 and D 2 are not remarkable. The mesopore volume, specific surface area, and average pore diameter of the samples with the hysteresis loop of types H3 and H4 are bimodal, but those of the samples with the hysteresis loop of type H2 are substantially unchanged with the increasing burial depth and Ro (Figures 9(b)-9(d)). Most organic pores are inkbottle-shaped pores that mainly occur in the samples with the hysteresis loop of type H2. Therefore, the mesopore evolution is not mainly controlled by organic matter transformation. These phenomena are different from the pore evolution model reflected by the marine New Albany shale [77] showing the unimodal trend over the similar maturity range (Ro = 0:55 -1:15%). The unimodal trend was interpreted as the result of organic matter transformation [77]. The pore evolution model from the marine New Albany shale [77] cannot be used to interpret the pore evolution of the lacustrine Es 4 U shale in the Dongying Depression. The main reasons probably are that the Es 4 U shale is lacustrine and that organic pores in the Es 4 U shale are not main contributors to the volume values of micro-and mesopores. In Figures 9(a)-9(c) the micro-and mesopore volumes as well as the specific surface area of the sample with the Ro value near 0.6% are very low, as are clay mineral and quartz content ( Figures S6(a) and S6(b)). But carbonates comprise more than 80% ( Figure S6(c)). The high carbonate content may lead to low pore volumes and specific surface areas. The average pore diameter is high, probably due to the fractures (i.e., slit-shaped pores) and less compaction.
When the thermal maturity reaches 0.7% Ro, the envelope curves of all micro-and mesopore volumes as well as specific surface area display peaks (Figures 9(a)-9(c)). At this stage, kerogen was generating hydrocarbons, and the clay minerals are dewatering, resulting in high pressures which prompted the formation of fractures in the Es 4 U shale [73,78]. Simultaneously, organic acids were produced from the organic-rich shale which can dissolve carbonates and feldspar [78]. These processes probably result in the increases of pore volumes and specific surface area (Figures 9(a)-9(c)) and a high average pore diameter (Figure 9(d)). There are only the samples with the hysteresis loop of type H3 near the first peaks of the mesopore volume, specific surface area, and average pore diameter envelopes (Figures 9(b)-9(d)). Fractures may be developed in this type of samples [40] and are favorable for the development of shale oil. One of the reasons for the high average pore diameter and mesopore volume probably is that the shale is not completely compacted, and so high pore volumes can be retained. Another reason is that fewer carbonates and more clay minerals and 15 Geofluids quartz in these samples ( Figure S6) are favorable for pore development. This is supported by the correlations that indicate micro-and mesopore volumes are all correlated positively to clay minerals and quartz and negatively to carbonates, as mentioned above.
Over the range of 0.85-0.90% Ro, the envelope curves for micro-and mesopore volumes, specific surface area, and average pore diameter, as well as carbonates, increase with increasing Ro (Figures 9(a)-9(d) and S6(c)). The high microand mesopore volumes as well as the specific surface area of   Figure 9: Cross-plots of main pore property parameters with Ro and burial depth for the Es 4 U shale. The samples with the hysteresis loops of types H2, H3, and H4 mainly contain inkbottle-, slit-, and narrow slit-shaped pores, respectively. The data are listed in Tables S1 and S2. 16 Geofluids the sample at the second peak (Ro = 0:9%) of these envelopes (Figures 9(a)-9(c)) may arise from the fracture formation after the carbonate cementation. This is supported by the high carbonate content and type H3 loop. The authigenic carbonates in the Es 4 U shale of the Dongying Depression were mainly formed during the period of 28.1-4.6 Ma [30,78]. The Bohai Bay Basin, including the Dongying Depression, was subjected to the tectonic movement since the Neogene (24.6 Ma to present), which may have resulted in the formation of fractures in the Es 4 U shale with the relatively high content of carbonates [79][80][81]. The sample (No. DYS6) at the second peaks of the micro-and mesopore volumes as well as APD envelope curves has modest micromesopore volume (0.036 cm 3 /g) and brittle mineral (71 wt%) and clay mineral content (29 wt%) (Tables S1 and S2) and may reflect a good shale reservoir for development. Over the range from 0.95% to 1.3% Ro (3904-4335 m), the envelope curves for the micro-and mesopore volumes and specific surface as well as the average pore diameter of the Soxhlet-extracted samples from the lacustrine Es 4 U shale all display the decreasing trends with the increasing Ro and burial depth (Figures 9(a)-9(d)). One reason probably is the strong compaction. With the increase of the depth, the carbonates and quartz in the shale become fewer but soft clay minerals increase ( Figure S6). This variation in mineral composition makes the rocks easily compacted under a high pressure. In the process of compaction, pores between clay minerals must be reduced greatly, and some fractures in carbonates may also be reduced or closed.

Conclusions
In this study, low-pressure gas (N 2 and CO 2 ) adsorption methods were applied to shale samples from the upper part of the Sha-4 Member of the Paleogene Shahejie Formation (Es 4 U) in the Dongying Depression, after the Soxhlet extraction of residual oil. On the basis of the results and discussion, the following conclusions have been reached: (1) The mesopore volume is much higher than the micropore volume in the Soxhlet-extracted samples. The Soxhlet-extracted samples of the Es 4 U shale have an average micropore volume of 0.0038 cm 3 /g and an average mesopore volume of 0.0182 cm 3 /g, higher than those of the Soxhlet-extracted samples of Chang-7 shale in the Ordos Basin. In the Es 4 U shale, the samples with type H3 hysteresis loop make up the most (53.3%) in the studied samples. Furthermore, our samples from the Es 4 U shale in the Dongying have an average brittle mineral (quartz, feldspar, and carbonates) content of 68.3 wt%, higher than that reported for the Chang-7 shale. Therefore, relatively high oil flow from the Es 4 U shale may potentially be achieved through horizontal drilling and hydraulic fracturing. However, the relationships of micro-and mesopore volumes with brittle minerals in the Es 4 U shale particularly need to be taken into account in the "sweet spot" selection and determination (2) PLSR (partial least square regression) analysis with VIP (variable importance in projection) evaluation is a powerful tool for the analyses of the main factors controlling shale pores. The PLSR results show that micro-and mesopore volumes, specific surface area, and the fractal dimension are mainly correlated positively to clay minerals and negatively to carbonate content. The increase in abundance of relatively small mesopores, related to clay minerals, is the main reason for high mesopore volume, with the high specific surface area and the fractal dimension reflecting the complexity of pore space structures. As a result, the samples with the relatively high mesopore volume for hydrocarbon storage have a high percentage of relatively small mesopores, complicated pore space structures, high content of clay minerals, and less carbonates. The shale represented by these samples may be unfavorable for shale oil development due to low oil flow, even if hydraulic fracturing is conducted. Therefore, the shale with the modest pore volumes and brittle and clay minerals may be conducive to the development of the Es 4 U shale oil (3) Organic pores in the Es 4 U shale mainly exist as micropores but are not a main contributor to mesopore volume. As the mesopore volume is much higher than the micropore volume, organic pores are not important for the storage space in the Soxhlet-extracted samples from the Es 4 U shale. Probably because of this, the IOC (insoluble organic carbon) is not a main factor controlling micromesopores in the lacustrine shale in the Dongying Depression (4) Over the maturity range of 0.6-1.3% Ro, the envelope curves of the pore volumes, specific surface area, and average pore diameter with the maturity show bimodal distributions with two peaks at 0.7% and 0.9% Ro. The samples at the first peak (0.7% Ro) have high pore volumes, specific surface area, and diameters with more clay minerals and less carbonates. The sample DYS6 from well L672 at the second peak (Ro = 0:9%) has relatively high micromesopore volume (0.036 cm 3 /g), specific surface area (18.85 m 2 /g), and average pore diameter (6.45 nm) with modest brittle mineral (71 wt%) and clay mineral content (29 wt%) and may reflect the shale favorable for shale oil development

Data Availability
The experimental data used to support the findings of this study are included within the manuscript and the supplementary materials.