# 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

Authors: Ziyi Wang; Liuping Zhang; Yuan Gao; Wenxiu Yang; Yongshi Wang; Huimin Liu; Zhaoyang Li; Shanshan Zhou
Journal: Geofluids
 (2021)

---

## Abstract

Lacustrine shales hold a huge potential oil resource in China. Pore properties (pore volume, diameter, specific surface area, and fractal dimensions) and their relationships with geological factors (mineralogy, insoluble organic carbon, burial depth, and vitrinite reflectance) are critical for evaluating shale oil resource. However, the factors controlling pores for lacustrine shale oil remain unclear, as the relationships between pore properties of Soxhlet-extracted samples and geological factors have not been studied using multivariate analytical methods. In this paper, the samples from the lacustrine shale in the upper part of the Sha-4 Member of the Paleogene Shahejie Formation in the Dongying Depression were tested with a set of experiments including Soxhlet (solvent) extraction, X-ray diffraction mineral analysis, insoluble organic carbon, vitrinite reflectance, and low-pressure CO2 and N2 adsorption experiments. The micromesopore volume varies from 0.003 cm3/g to 0.045 cm3/g. The relationships of pore properties with geological factors were studied with partial least square regression analysis (PLSR analysis, a powerful multivariate regression analysis). The results of the PLSR analyses indicate that clay minerals and carbonates are two key factors affecting the pore properties of the lacustrine shale. Compared with marine shales, more clay minerals in the lacustrine shale make them become more important for controlling pores than organic matter. The PLSR results also illustrate that the shale with higher pore volume contains more clay minerals and fewer carbonates and thus is unfavorable for hydraulic fracturing. Therefore, the shale with high micromesopore volume may be unfavorable for shale oil production. The shale with the modest micromesopore volume (~0.036 cm3/g), relatively high content of brittle minerals (~71 wt%), and low clay mineral content (~29 wt%) is conducive to both oil storage and hydraulic fracturing for the development of the Es4U shale oil in the Dongying Depression in East China.


---

## Body

## 1. 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 exploration of shale oil and gas in basins worldwide. It has been conservatively estimated that China has shale oil resources of 40.3×109tons, and lacustrine shale holds the vast majority of the shale oil resources in China (40.2×109tons) [3]. However, the studies of lacustrine shale oil reservoirs are very limited, compared with many studies of marine shale gas reservoirs [4–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 (withdiameter<2nm), mesopores (with diameter of 2–50 nm), and macropores (with diameter>50nm) [11]. Among the techniques for quantitative analysis of shale pore systems [4, 12–14], low-pressure gas adsorption (LPGA) is important in the measurement of pore properties over a wide range of pore sizes. Low-pressure N2 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 N2 adsorption isotherm [17–19]. However, N2 adsorption is less accurate in characterizing micropores due to the limitation arising from N2 molecule and pore throat sizes. A CO2 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 N2 molecule (at 77.35 K) [20]. Low-pressure CO2 adsorption can be used to analyze micropores. Therefore, both N2 and CO2 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 N2 and CO2 techniques. Removing the residual oil in shale samples via the Soxhlet (solvent) extraction can solve this problem in the analysis of LPGA [21–23]. However, published data on shale pore properties in the oil window are rarely obtained from solvent-extracted samples [24–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 (Es4U) 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 km2 for the shale thicker than 50 m. In this work, the pore properties in fifteen Soxhlet-extracted samples of the Es4U shale were analyzed using N2 and CO2 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 Es4U shale oil.


## 2. 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 (Es4), -3 (Es3), -2 (Es2), and -1 (Es1) Members from bottom to top. The shale in the upper part of the Sha-4 Member (Es4U) is an important source rock in the Dongying Depression [22]. The thickness of Es4U shale is up to 400 m [38]. The area with thickness greater than 50 m (154 m on average) occupies about 3000 km2 and is mainly located in the central-north part of the Dongying Depression (Figure 1). The thick Es4U 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., S1+S2) values of 0.62–76.51 mg/g [29, 36–38, 40]. The average oil saturation of Es4U 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 Es4U deposition [42], facilitating the development of carbonates (see the strata column in Zhang et al. [35]).Figure 1
Outline map of the Dongying Depression, showing the isopachs (m) of the Es4U shale (data from Zhang et al. [43]), the Ro (%) contours of the bottom boundary of the Es4U shale (data from the database of Sinopec Shengli Oilfield Company) and the locations of sampling wells.


## 3. Samples and Methods

### 3.1. Samples

To reveal the pore properties of the main bodies of the shale, fifteen core samples of the Es4U shale were selected from nine wells (Figures 1 and S1), based on the sample locations in Es4U, 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).Figure 2
Thin-section photomicrographs of the Es4U shale samples in the Dongying Depression: (a) monopolarization photo of the interlaminated mudstone and calcite (dyed red by Alizarin red) in the shale sample DYS2 of the FS1 well at 3688.0 m (burialdepth=3679.4m); (b) monopolarization photo of the silty sands in the laminated shale sample DYS12 of the T746 well at 3545.2 m (burialdepth=3538.7m); (c) monopolarization photo of the massive structure in the mudstone sample DYS11 of the T720 well at 3677.3 m (burialdepth=3666.9m); (d) monopolarization photo of the sample DYS9 of the L988 well at 3979.3 m (burialdepth=3966.1m); (e) ultraviolet fluorescence photo of the original sample DYS9 showing the residual oil in the shale; (f) ultraviolet fluorescence photo of the Soxhlet-extracted sample showing most residual oil in the sample DYS9 was removed; (g) SEM image showing inkbottle-shaped pores between pyrites in the sample DYS10 of the L988 well at 4131.9 m (burialdepth=4117.1m); (h) SEM image showing slit-shaped pores (fractures) in the sample DYS3 of the FS1 well at 3819.6 m (burialdepth=3810.9m); (i) SEM image showing narrow slit-shaped pore between clay minerals in the sample DYS11.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


### 3.2. 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 Soxhlet-extracted 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 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 Soxhlet-extracted 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 N2 and CO2 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.01torr) for 20 h to remove residual volatile material before the measurements. The N2 and CO2 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/P0) in N2 and CO2 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 N2 and CO2 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 N2 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 N2 adsorption pore volume and specific surface area [53].The specific surface areas tested by CO2 adsorption were not reported in this paper, due to large uncertainties in the area occupied by an adsorbed CO2 molecule in the pores of molecular dimensions and the CO2-measured monolayer capacity (the amount of gas molecules of a complete, close-packed 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 N2 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:
(1)lnV=K·lnlnP0/P+C,(2)D=K+3,where K is the slope of the line of lnV vs. lnlnP0/P, V is the adsorbed nitrogen volume (cm3/g) at the equilibrium pressure P (Pa), P0 (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 lnV vs. lnlnP0/P line shows two segments with different slopes, two kinds of fractal dimensions (D1 and D2) need to be calculated also using Equations (1) and (2). The fractal dimension D1, calculated with the adsorption data at P/P0<0.45, reflects the complexity of pore surfaces, and D2, calculated with the adsorption data at P/P0>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 (Vmic), mesopore volume (Vmes), micromesopore volume (Vmic‐mes, sum of Vmic and Vmes), specific surface area (S), average pore diameter (APD), and fractal dimensions (D1 and D2), with the appropriate forms (e.g., lnVmes and lnS; 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:
(3)Fi=∑j=1najxj,where Fi is the ith component (i=1,2,⋯,m), xj is the jth standardized descriptor (j=1,2,⋯,n), and aj is a coefficient. Then, the regression linear equation of a standardized response versus one or several orthogonal components is established:
(4)y=∑i=1mbiFi+C1,where y is the standardized response, bj is a coefficient, and C1 is a residual. The linear equation of a response versus descriptors is obtained from Equations (3) and (4):
(5)y=∑j=1ncjxj+C2,where cj is a coefficient and C2 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].


## 3.1. Samples

To reveal the pore properties of the main bodies of the shale, fifteen core samples of the Es4U shale were selected from nine wells (Figures 1 and S1), based on the sample locations in Es4U, 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).Figure 2
Thin-section photomicrographs of the Es4U shale samples in the Dongying Depression: (a) monopolarization photo of the interlaminated mudstone and calcite (dyed red by Alizarin red) in the shale sample DYS2 of the FS1 well at 3688.0 m (burialdepth=3679.4m); (b) monopolarization photo of the silty sands in the laminated shale sample DYS12 of the T746 well at 3545.2 m (burialdepth=3538.7m); (c) monopolarization photo of the massive structure in the mudstone sample DYS11 of the T720 well at 3677.3 m (burialdepth=3666.9m); (d) monopolarization photo of the sample DYS9 of the L988 well at 3979.3 m (burialdepth=3966.1m); (e) ultraviolet fluorescence photo of the original sample DYS9 showing the residual oil in the shale; (f) ultraviolet fluorescence photo of the Soxhlet-extracted sample showing most residual oil in the sample DYS9 was removed; (g) SEM image showing inkbottle-shaped pores between pyrites in the sample DYS10 of the L988 well at 4131.9 m (burialdepth=4117.1m); (h) SEM image showing slit-shaped pores (fractures) in the sample DYS3 of the FS1 well at 3819.6 m (burialdepth=3810.9m); (i) SEM image showing narrow slit-shaped pore between clay minerals in the sample DYS11.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


## 3.2. 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 Soxhlet-extracted 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 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 Soxhlet-extracted 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 N2 and CO2 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.01torr) for 20 h to remove residual volatile material before the measurements. The N2 and CO2 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/P0) in N2 and CO2 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 N2 and CO2 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 N2 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 N2 adsorption pore volume and specific surface area [53].The specific surface areas tested by CO2 adsorption were not reported in this paper, due to large uncertainties in the area occupied by an adsorbed CO2 molecule in the pores of molecular dimensions and the CO2-measured monolayer capacity (the amount of gas molecules of a complete, close-packed 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 N2 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:
(1)lnV=K·lnlnP0/P+C,(2)D=K+3,where K is the slope of the line of lnV vs. lnlnP0/P, V is the adsorbed nitrogen volume (cm3/g) at the equilibrium pressure P (Pa), P0 (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 lnV vs. lnlnP0/P line shows two segments with different slopes, two kinds of fractal dimensions (D1 and D2) need to be calculated also using Equations (1) and (2). The fractal dimension D1, calculated with the adsorption data at P/P0<0.45, reflects the complexity of pore surfaces, and D2, calculated with the adsorption data at P/P0>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 (Vmic), mesopore volume (Vmes), micromesopore volume (Vmic‐mes, sum of Vmic and Vmes), specific surface area (S), average pore diameter (APD), and fractal dimensions (D1 and D2), with the appropriate forms (e.g., lnVmes and lnS; 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:
(3)Fi=∑j=1najxj,where Fi is the ith component (i=1,2,⋯,m), xj is the jth standardized descriptor (j=1,2,⋯,n), and aj is a coefficient. Then, the regression linear equation of a standardized response versus one or several orthogonal components is established:
(4)y=∑i=1mbiFi+C1,where y is the standardized response, bj is a coefficient, and C1 is a residual. The linear equation of a response versus descriptors is obtained from Equations (3) and (4):
(5)y=∑j=1ncjxj+C2,where cj is a coefficient and C2 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].


## 4. Results and Discussion

### 4.1. Mineral Composition, Insoluble Organic Carbon, and Thermal Maturity

The results of the XRD analysis for the mineral composition are listed in TableS1. 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 Es4U 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 Es4U. 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 Es4U deposition [42].Carbonates are negatively related to clay minerals and quartz (FiguresS2(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 Es4U 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 Es4U 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 Es4U shale is favorable for hydraulic fracturing.The other shales were studied with Soxhlet-unextracted samples (see Section4.2.2 for details), such as the Bakken shale in the Williston Basin and Kong-2 shale in the Cangdong Depression [25–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 samples with those in the Soxhlet-extracted samples in our study, as this comparison is not helpful to study the pore properties of the Es4U shale.The IOC content of the Es4U 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 Es4U shale is lower than that of the Chang-7 shale (averageIOC=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 Es4U 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 TableS1. 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 Es4U shale samples are in the oil window (Ro=0.6–1.0%) and some in the gas condensate-wet gas window (Ro=1.0–1.3%) according to U.S. EIA [63] and Dembicki [64].


### 4.2. Pore Properties

Before low-pressure gas adsorption experiments for shale pores, the effect of residual oil removal was checked with the thin sections of all original (Soxhlet-unextracted) and Soxhlet-extracted shale samples. Figures2(d)–2(f) are the thin-section photomicrographs of the sample DYS9. Ultraviolet fluorescence photos indicate there is residual oil with yellow-green fluorescence in the unextracted sample (Figure 2(e)), but the residual oil is hardly observed in the Soxhlet-extracted sample (Figure 2(f)). Therefore, most of the residual oil in our samples was removed. Then, all of fifteen Soxhlet-extracted samples were analyzed using low pressure CO2 and N2 adsorption.


#### 4.2.1. CO2 and N2 Adsorption Isotherms

The total adsorbed CO2 volumes in the CO2 adsorption experiments range from 0.15 cm3/g to 2.07 cm3/g at standard temperature and pressure (STP). The CO2 adsorption isotherms of all samples are convex in shape when plotted against the relative pressure (P/P0) (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 N2 adsorption experiments, the total adsorbed N2 volumes vary from 2.12 cm3/g to 29.61 cm3/g at STP. The N2 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 CO2 and N2 adsorption isotherms indicate that the samples from the Es4U all contain micro- and mesopores.Figure 3
N2 adsorption-desorption isotherms of original and Soxhlet-extracted samples from the Dongying Es4U shale with typical types of hysteresis loops: (a) sample DYS10, (b) sample DYS3, and (c) sample DYS11.
(a)(b)(c)To further check the effect of the residual oil removal, three original samples (DYS3, DYS10, and DYS11) were also analyzed. Both total adsorbed CO2 and N2 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 N2 hysteresis loops are categorized into types H1, H2, H3, and H4 according to the IUPAC classification [11]. For the Soxhlet-extracted samples from the Es4U 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/P0=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 Es4U shales. The hydrocarbon flow capability of different types of pores is discussed together with the fractal characteristics in Section 4.2.3.


#### 4.2.2. 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 Section3.5 and data from N2 and CO2 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 TableS2. 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 Es4U. 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).Figure 4
Pore size distribution of the original and Soxhlet-extracted samples from the Dongying Es4U shale: (a) sample DYS10, (b) sample DYS3, and (c) sample DYS11. V: volume; D: diameter; dV/dD: the derivative of volume.
(a)(b)(c)For the Soxhlet-extracted shale samples in Es4U, the micropore volumes (Table S2) vary from 0.0006 cm3/g (sample DYS14) to 0.0071 cm3/g (sample DYS6) (0.0038 cm3/g on average). The mesopore volumes range from 0.0028 cm3/g (sample DYS14) to 0.0381 cm3/g (sample DYS15) with an average of 0.0181 cm3/g. The mesopore volume is much higher than the micropore volume. The micromesopore volume (sum of micropore and mesopore volumes) varies from 0.003 cm3/g to 0.045 cm3/g with an average of 0.022 cm3/g. The specific surface areas, tested by N2 adsorption, range from 1.04 m2/g (sample DYS14) to 30.55 m2/g (sample DYS15) with an average of 12.56 m2/g. The average pore diameters of the samples, obtained from the N2 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 CO2 adsorption data, due to their large uncertainties (see subsection 3.5).To understand the pores in the Es4U shale, the pore parameters of the Es4U 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–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 N2 and CO2 adsorption methods, after the solvent extraction [23, 44]. The Chang-7 shale has the average micropore volume of 0.0035 cm3/g, mesopore volume of 0.0113 cm3/g, specific surface area of 4.99 m2/g, and average pore diameter of 17.54 nm [23, 44]. Compared with the Soxhlet-extracted samples from the Chang-7 shale, the Es4U shale has significantly higher micro- and mesopore volumes (0.0038 cm3/g and 0.0181 cm3/g on average, respectively), but the average pore diameter is lower. The average oil saturation of mature Es4U shale (45.8%) [41] is close to that in the Chang-7 shale (44.63%) [66]. The Es4U 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 Es4U shale are mudstone, shale, and gypsum [38, 69]. All the characteristics of shale reservoirs except the average pore diameter display that the Dongying Es4U shale is favorable for the development of shale oil and may achieve higher oil flow than the commercial oil flow from the Chang-7 shale in the Ordos Basin. But the influence of the average pore diameter on the oil flow needs to be studied further.


#### 4.2.3. Fractal Dimensions

The fractal dimension values of the Soxhlet-extracted samples from the Es4U shale were calculated using the FHH method [19] applied to the N2 adsorption data (see Section 3.6 for details). The lnV vs. lnlnP0/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 (D1 and D2) were calculated for each sample. The values of the fractal dimension D1, calculated with the adsorption data at P/P0<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 D1 reflects the complexity of pore surfaces [58]. The D2 values, calculated with the adsorption data at P/P0>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 D2 reflects the complexity of pore space structures [58]. Higher D1 and D2 indicate that the samples have rougher pore surfaces and more complicated pore space structures, respectively [24]. Therefore, pores with low D1 and D2 values are favorable for oil flow.Figure 5
Plot of lnV vs. ln[ln(P0/P)] from the N2 adsorption isotherm of the Soxhlet-extracted sample DYS1 (a) and average fractal dimensions of the pores in the Soxhlet-extracted samples with the hysteresis loops of types H2, H3, and H4 that mainly contain inkbottle-, slit-, and narrow slit-shaped pores, respectively (b). The data are listed in Table S2.
(a)(b)Figure5(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 Es4U 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 D1 and D2 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, inkbottle- and 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.


#### 4.2.4. Correlations between Pore Parameters

Figures6(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
Correlation relationships between pore property parameters and between geological factors and pore property parameters for the Soxhlet-extracted samples from the Es4U shale: (a–d) between pore parameters (mesopore volume, specific surface area, average pore diameter, and fractal dimensions D1 and D2); (e–h) between pore parameters (micro- and mesopore volumes, average pore diameter, and fractal dimension D1) 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.
(a)(b)(c)(d)(e)(f)(g)(h)Figure6(c) shows that the average pore diameter has a significantly negative correlation with fractal dimension D2 (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 D1 (Figure 6(d)), which is also supported by the positive correlation between D1 and D2 (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)), D1 and D2 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 D1 and D2.In summary, the increase of small pores in the shale leads to the increases of the pore volumes, specific surface area, and fractal dimensions (D1 and D2) and the decrease of the average pore diameter. But the correlations of D1 with average pore diameter (R=0.52, Figure 6(d)) and D2 (R=0.52, Figure S5(b)) are not very evident. This may suggest that the D1 of the Es4U shale is also affected by other factors, such as mineral composition, burial depth, and Ro (see below).


### 4.3. 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 (TableS1). 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 Section3.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 D1 (Figure S5(d)) and D2 (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.


#### 4.3.1. Geological Factors Controlling Micropore Volume

The log of the micropore volume (lnVmic) was set as a response (dependent variable) in the PLSR analysis. The descriptors (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 lnVmic 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 Vmic 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)).Figure 7
Correlation relationships between experimental and predicted pore property parameters (Vmic, Vmes, APD, and D1) in the PLSR fitting models. Vmic: micropore volume; Vmes: mesopore volume; APD: average pore diameter; D1: fractal dimension D1. 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.
(a)(b)(c)(d)The standardized coefficients of the multivariate regression equation and VIP values (Figure8(a) and Table S3) of the PLSR fitting model of lnVmic 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 Es4U 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 Es4U 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 Es4U shale.Figure 8
Variable importance in projection (VIP) of geological factors for the PLSR fitting models of pore properties.Vmic: micropore volume; Vmes: mesopore volume; Vmic‐mes: micromesopore volume, sum of Vmic and Vmes; S: specific surface area; APD: average pore diameter; D1 and D2: fractal dimensions. The data are listed in Table S3.
(a)(b)(c)(d)


#### 4.3.2. Geological Factors Controlling Mesopore Volume and Specific Surface Area

TableS3 also shows the standardized coefficients and VIP values of the descriptors obtained from the PLSR analysis of the log of the mesopore volume (lnVmes). The predicted values of the PLSR fitting model of lnVmes (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 Vmes (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 multivariate regression equation (Figure 8(a) and Table S3) of the PLSR fitting model of lnVmes 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 lnVmic, 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 Es4U 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 Es4U 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 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 (lnS) and micromesopore volume (lnVmic‐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 lnVmic‐mes and lnS are close to those in the lnVmes 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 lnVmic‐mes and lnS (see Section 3.7 for details). These features are the same as those in the PLSR fitting model of lnVmes.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 volumes are controlled or affected by multiple geological factors. The mesopore and micromesopore volumes as well as the specific surface area of the Dongying Es4U 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 Es4U shale need to be synthetically taken into account in the “sweet spot” selection and determination.


#### 4.3.3. Geological Factors Controlling Geometrical Parameters of Mesopores

The geometrical parameters of pores including average pore diameter (APD) and fractal dimensionsD1 and D2 were set as the responses for the PLSR analyses. The 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 ofD2 (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 D2 are similar to those controlling APD. But the effects of these geological factors on the response (D2 or APD) are the opposite, as APD has a negative correlation with D2 (Figure 6(c)).Figure8(d) and Table S3 also show the standardized coefficients and VIP values of the PLSR analysis for D1. The relationship between the experimental and predicted values of the PLSR fitting model of D1 has a higher correlation coefficient (R=0.72, Figure 7(d)) than those of univariate regression analyses (R=0.61, Figures 6(h) and S5(l)). The dolomite and clay minerals with VIP values higher than 1 (Figure 8(d) and Table S3) have relatively important influences on D1 (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 Es4U shale is mainly self-structured crystals with the flat surface [75]. Therefore, D1 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 dimensionD2 of the Es4U shale are mainly controlled by multiple geological factors including clay minerals, feldspar, burial depth, lnRo, carbonates, and calcite, while the fractal dimension D1 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 Es4U solvent-extracted samples, absolute pore volumes were directly used to analyze the pore evolution, instead of TOC- or 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.Figure 9
Cross-plots of main pore property parameters with Ro and burial depth for the Es4U 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.
(a)(b)(c)(d)(e)(f)Over the range of 0.6–1.3% Ro, Figures9(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 D1 and D2 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 Es4U shale in the Dongying Depression. The main reasons probably are that the Es4U shale is lacustrine and that organic pores in the Es4U 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 (Figures9(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 Es4U 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 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 (Figures9(a)–9(d) and S6(c)). The high micro- and mesopore volumes as well as the specific surface area of 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 Es4U 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 Es4U shale with the relatively high content of carbonates [79–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 cm3/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 Es4U 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.


## 4.1. Mineral Composition, Insoluble Organic Carbon, and Thermal Maturity

The results of the XRD analysis for the mineral composition are listed in TableS1. 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 Es4U 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 Es4U. 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 Es4U deposition [42].Carbonates are negatively related to clay minerals and quartz (FiguresS2(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 Es4U 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 Es4U 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 Es4U shale is favorable for hydraulic fracturing.The other shales were studied with Soxhlet-unextracted samples (see Section4.2.2 for details), such as the Bakken shale in the Williston Basin and Kong-2 shale in the Cangdong Depression [25–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 samples with those in the Soxhlet-extracted samples in our study, as this comparison is not helpful to study the pore properties of the Es4U shale.The IOC content of the Es4U 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 Es4U shale is lower than that of the Chang-7 shale (averageIOC=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 Es4U 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 TableS1. 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 Es4U shale samples are in the oil window (Ro=0.6–1.0%) and some in the gas condensate-wet gas window (Ro=1.0–1.3%) according to U.S. EIA [63] and Dembicki [64].


## 4.2. Pore Properties

Before low-pressure gas adsorption experiments for shale pores, the effect of residual oil removal was checked with the thin sections of all original (Soxhlet-unextracted) and Soxhlet-extracted shale samples. Figures2(d)–2(f) are the thin-section photomicrographs of the sample DYS9. Ultraviolet fluorescence photos indicate there is residual oil with yellow-green fluorescence in the unextracted sample (Figure 2(e)), but the residual oil is hardly observed in the Soxhlet-extracted sample (Figure 2(f)). Therefore, most of the residual oil in our samples was removed. Then, all of fifteen Soxhlet-extracted samples were analyzed using low pressure CO2 and N2 adsorption.


### 4.2.1. CO2 and N2 Adsorption Isotherms

The total adsorbed CO2 volumes in the CO2 adsorption experiments range from 0.15 cm3/g to 2.07 cm3/g at standard temperature and pressure (STP). The CO2 adsorption isotherms of all samples are convex in shape when plotted against the relative pressure (P/P0) (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 N2 adsorption experiments, the total adsorbed N2 volumes vary from 2.12 cm3/g to 29.61 cm3/g at STP. The N2 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 CO2 and N2 adsorption isotherms indicate that the samples from the Es4U all contain micro- and mesopores.Figure 3
N2 adsorption-desorption isotherms of original and Soxhlet-extracted samples from the Dongying Es4U shale with typical types of hysteresis loops: (a) sample DYS10, (b) sample DYS3, and (c) sample DYS11.
(a)(b)(c)To further check the effect of the residual oil removal, three original samples (DYS3, DYS10, and DYS11) were also analyzed. Both total adsorbed CO2 and N2 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 N2 hysteresis loops are categorized into types H1, H2, H3, and H4 according to the IUPAC classification [11]. For the Soxhlet-extracted samples from the Es4U 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/P0=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 Es4U shales. The hydrocarbon flow capability of different types of pores is discussed together with the fractal characteristics in Section 4.2.3.


### 4.2.2. 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 Section3.5 and data from N2 and CO2 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 TableS2. 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 Es4U. 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).Figure 4
Pore size distribution of the original and Soxhlet-extracted samples from the Dongying Es4U shale: (a) sample DYS10, (b) sample DYS3, and (c) sample DYS11. V: volume; D: diameter; dV/dD: the derivative of volume.
(a)(b)(c)For the Soxhlet-extracted shale samples in Es4U, the micropore volumes (Table S2) vary from 0.0006 cm3/g (sample DYS14) to 0.0071 cm3/g (sample DYS6) (0.0038 cm3/g on average). The mesopore volumes range from 0.0028 cm3/g (sample DYS14) to 0.0381 cm3/g (sample DYS15) with an average of 0.0181 cm3/g. The mesopore volume is much higher than the micropore volume. The micromesopore volume (sum of micropore and mesopore volumes) varies from 0.003 cm3/g to 0.045 cm3/g with an average of 0.022 cm3/g. The specific surface areas, tested by N2 adsorption, range from 1.04 m2/g (sample DYS14) to 30.55 m2/g (sample DYS15) with an average of 12.56 m2/g. The average pore diameters of the samples, obtained from the N2 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 CO2 adsorption data, due to their large uncertainties (see subsection 3.5).To understand the pores in the Es4U shale, the pore parameters of the Es4U 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–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 N2 and CO2 adsorption methods, after the solvent extraction [23, 44]. The Chang-7 shale has the average micropore volume of 0.0035 cm3/g, mesopore volume of 0.0113 cm3/g, specific surface area of 4.99 m2/g, and average pore diameter of 17.54 nm [23, 44]. Compared with the Soxhlet-extracted samples from the Chang-7 shale, the Es4U shale has significantly higher micro- and mesopore volumes (0.0038 cm3/g and 0.0181 cm3/g on average, respectively), but the average pore diameter is lower. The average oil saturation of mature Es4U shale (45.8%) [41] is close to that in the Chang-7 shale (44.63%) [66]. The Es4U 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 Es4U shale are mudstone, shale, and gypsum [38, 69]. All the characteristics of shale reservoirs except the average pore diameter display that the Dongying Es4U shale is favorable for the development of shale oil and may achieve higher oil flow than the commercial oil flow from the Chang-7 shale in the Ordos Basin. But the influence of the average pore diameter on the oil flow needs to be studied further.


### 4.2.3. Fractal Dimensions

The fractal dimension values of the Soxhlet-extracted samples from the Es4U shale were calculated using the FHH method [19] applied to the N2 adsorption data (see Section 3.6 for details). The lnV vs. lnlnP0/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 (D1 and D2) were calculated for each sample. The values of the fractal dimension D1, calculated with the adsorption data at P/P0<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 D1 reflects the complexity of pore surfaces [58]. The D2 values, calculated with the adsorption data at P/P0>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 D2 reflects the complexity of pore space structures [58]. Higher D1 and D2 indicate that the samples have rougher pore surfaces and more complicated pore space structures, respectively [24]. Therefore, pores with low D1 and D2 values are favorable for oil flow.Figure 5
Plot of lnV vs. ln[ln(P0/P)] from the N2 adsorption isotherm of the Soxhlet-extracted sample DYS1 (a) and average fractal dimensions of the pores in the Soxhlet-extracted samples with the hysteresis loops of types H2, H3, and H4 that mainly contain inkbottle-, slit-, and narrow slit-shaped pores, respectively (b). The data are listed in Table S2.
(a)(b)Figure5(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 Es4U 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 D1 and D2 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, inkbottle- and 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.


### 4.2.4. Correlations between Pore Parameters

Figures6(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
Correlation relationships between pore property parameters and between geological factors and pore property parameters for the Soxhlet-extracted samples from the Es4U shale: (a–d) between pore parameters (mesopore volume, specific surface area, average pore diameter, and fractal dimensions D1 and D2); (e–h) between pore parameters (micro- and mesopore volumes, average pore diameter, and fractal dimension D1) 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.
(a)(b)(c)(d)(e)(f)(g)(h)Figure6(c) shows that the average pore diameter has a significantly negative correlation with fractal dimension D2 (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 D1 (Figure 6(d)), which is also supported by the positive correlation between D1 and D2 (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)), D1 and D2 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 D1 and D2.In summary, the increase of small pores in the shale leads to the increases of the pore volumes, specific surface area, and fractal dimensions (D1 and D2) and the decrease of the average pore diameter. But the correlations of D1 with average pore diameter (R=0.52, Figure 6(d)) and D2 (R=0.52, Figure S5(b)) are not very evident. This may suggest that the D1 of the Es4U shale is also affected by other factors, such as mineral composition, burial depth, and Ro (see below).


## 4.2.1. CO2 and N2 Adsorption Isotherms

The total adsorbed CO2 volumes in the CO2 adsorption experiments range from 0.15 cm3/g to 2.07 cm3/g at standard temperature and pressure (STP). The CO2 adsorption isotherms of all samples are convex in shape when plotted against the relative pressure (P/P0) (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 N2 adsorption experiments, the total adsorbed N2 volumes vary from 2.12 cm3/g to 29.61 cm3/g at STP. The N2 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 CO2 and N2 adsorption isotherms indicate that the samples from the Es4U all contain micro- and mesopores.Figure 3
N2 adsorption-desorption isotherms of original and Soxhlet-extracted samples from the Dongying Es4U shale with typical types of hysteresis loops: (a) sample DYS10, (b) sample DYS3, and (c) sample DYS11.
(a)(b)(c)To further check the effect of the residual oil removal, three original samples (DYS3, DYS10, and DYS11) were also analyzed. Both total adsorbed CO2 and N2 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 N2 hysteresis loops are categorized into types H1, H2, H3, and H4 according to the IUPAC classification [11]. For the Soxhlet-extracted samples from the Es4U 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/P0=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 Es4U shales. The hydrocarbon flow capability of different types of pores is discussed together with the fractal characteristics in Section 4.2.3.


## 4.2.2. 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 Section3.5 and data from N2 and CO2 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 TableS2. 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 Es4U. 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).Figure 4
Pore size distribution of the original and Soxhlet-extracted samples from the Dongying Es4U shale: (a) sample DYS10, (b) sample DYS3, and (c) sample DYS11. V: volume; D: diameter; dV/dD: the derivative of volume.
(a)(b)(c)For the Soxhlet-extracted shale samples in Es4U, the micropore volumes (Table S2) vary from 0.0006 cm3/g (sample DYS14) to 0.0071 cm3/g (sample DYS6) (0.0038 cm3/g on average). The mesopore volumes range from 0.0028 cm3/g (sample DYS14) to 0.0381 cm3/g (sample DYS15) with an average of 0.0181 cm3/g. The mesopore volume is much higher than the micropore volume. The micromesopore volume (sum of micropore and mesopore volumes) varies from 0.003 cm3/g to 0.045 cm3/g with an average of 0.022 cm3/g. The specific surface areas, tested by N2 adsorption, range from 1.04 m2/g (sample DYS14) to 30.55 m2/g (sample DYS15) with an average of 12.56 m2/g. The average pore diameters of the samples, obtained from the N2 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 CO2 adsorption data, due to their large uncertainties (see subsection 3.5).To understand the pores in the Es4U shale, the pore parameters of the Es4U 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–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 N2 and CO2 adsorption methods, after the solvent extraction [23, 44]. The Chang-7 shale has the average micropore volume of 0.0035 cm3/g, mesopore volume of 0.0113 cm3/g, specific surface area of 4.99 m2/g, and average pore diameter of 17.54 nm [23, 44]. Compared with the Soxhlet-extracted samples from the Chang-7 shale, the Es4U shale has significantly higher micro- and mesopore volumes (0.0038 cm3/g and 0.0181 cm3/g on average, respectively), but the average pore diameter is lower. The average oil saturation of mature Es4U shale (45.8%) [41] is close to that in the Chang-7 shale (44.63%) [66]. The Es4U 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 Es4U shale are mudstone, shale, and gypsum [38, 69]. All the characteristics of shale reservoirs except the average pore diameter display that the Dongying Es4U shale is favorable for the development of shale oil and may achieve higher oil flow than the commercial oil flow from the Chang-7 shale in the Ordos Basin. But the influence of the average pore diameter on the oil flow needs to be studied further.


## 4.2.3. Fractal Dimensions

The fractal dimension values of the Soxhlet-extracted samples from the Es4U shale were calculated using the FHH method [19] applied to the N2 adsorption data (see Section 3.6 for details). The lnV vs. lnlnP0/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 (D1 and D2) were calculated for each sample. The values of the fractal dimension D1, calculated with the adsorption data at P/P0<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 D1 reflects the complexity of pore surfaces [58]. The D2 values, calculated with the adsorption data at P/P0>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 D2 reflects the complexity of pore space structures [58]. Higher D1 and D2 indicate that the samples have rougher pore surfaces and more complicated pore space structures, respectively [24]. Therefore, pores with low D1 and D2 values are favorable for oil flow.Figure 5
Plot of lnV vs. ln[ln(P0/P)] from the N2 adsorption isotherm of the Soxhlet-extracted sample DYS1 (a) and average fractal dimensions of the pores in the Soxhlet-extracted samples with the hysteresis loops of types H2, H3, and H4 that mainly contain inkbottle-, slit-, and narrow slit-shaped pores, respectively (b). The data are listed in Table S2.
(a)(b)Figure5(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 Es4U 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 D1 and D2 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, inkbottle- and 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.


## 4.2.4. Correlations between Pore Parameters

Figures6(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
Correlation relationships between pore property parameters and between geological factors and pore property parameters for the Soxhlet-extracted samples from the Es4U shale: (a–d) between pore parameters (mesopore volume, specific surface area, average pore diameter, and fractal dimensions D1 and D2); (e–h) between pore parameters (micro- and mesopore volumes, average pore diameter, and fractal dimension D1) 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.
(a)(b)(c)(d)(e)(f)(g)(h)Figure6(c) shows that the average pore diameter has a significantly negative correlation with fractal dimension D2 (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 D1 (Figure 6(d)), which is also supported by the positive correlation between D1 and D2 (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)), D1 and D2 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 D1 and D2.In summary, the increase of small pores in the shale leads to the increases of the pore volumes, specific surface area, and fractal dimensions (D1 and D2) and the decrease of the average pore diameter. But the correlations of D1 with average pore diameter (R=0.52, Figure 6(d)) and D2 (R=0.52, Figure S5(b)) are not very evident. This may suggest that the D1 of the Es4U shale is also affected by other factors, such as mineral composition, burial depth, and Ro (see below).


## 4.3. 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 (TableS1). 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 Section3.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 D1 (Figure S5(d)) and D2 (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.


### 4.3.1. Geological Factors Controlling Micropore Volume

The log of the micropore volume (lnVmic) was set as a response (dependent variable) in the PLSR analysis. The descriptors (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 lnVmic 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 Vmic 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)).Figure 7
Correlation relationships between experimental and predicted pore property parameters (Vmic, Vmes, APD, and D1) in the PLSR fitting models. Vmic: micropore volume; Vmes: mesopore volume; APD: average pore diameter; D1: fractal dimension D1. 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.
(a)(b)(c)(d)The standardized coefficients of the multivariate regression equation and VIP values (Figure8(a) and Table S3) of the PLSR fitting model of lnVmic 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 Es4U 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 Es4U 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 Es4U shale.Figure 8
Variable importance in projection (VIP) of geological factors for the PLSR fitting models of pore properties.Vmic: micropore volume; Vmes: mesopore volume; Vmic‐mes: micromesopore volume, sum of Vmic and Vmes; S: specific surface area; APD: average pore diameter; D1 and D2: fractal dimensions. The data are listed in Table S3.
(a)(b)(c)(d)


### 4.3.2. Geological Factors Controlling Mesopore Volume and Specific Surface Area

TableS3 also shows the standardized coefficients and VIP values of the descriptors obtained from the PLSR analysis of the log of the mesopore volume (lnVmes). The predicted values of the PLSR fitting model of lnVmes (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 Vmes (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 multivariate regression equation (Figure 8(a) and Table S3) of the PLSR fitting model of lnVmes 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 lnVmic, 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 Es4U 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 Es4U 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 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 (lnS) and micromesopore volume (lnVmic‐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 lnVmic‐mes and lnS are close to those in the lnVmes 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 lnVmic‐mes and lnS (see Section 3.7 for details). These features are the same as those in the PLSR fitting model of lnVmes.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 volumes are controlled or affected by multiple geological factors. The mesopore and micromesopore volumes as well as the specific surface area of the Dongying Es4U 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 Es4U shale need to be synthetically taken into account in the “sweet spot” selection and determination.


### 4.3.3. Geological Factors Controlling Geometrical Parameters of Mesopores

The geometrical parameters of pores including average pore diameter (APD) and fractal dimensionsD1 and D2 were set as the responses for the PLSR analyses. The 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 ofD2 (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 D2 are similar to those controlling APD. But the effects of these geological factors on the response (D2 or APD) are the opposite, as APD has a negative correlation with D2 (Figure 6(c)).Figure8(d) and Table S3 also show the standardized coefficients and VIP values of the PLSR analysis for D1. The relationship between the experimental and predicted values of the PLSR fitting model of D1 has a higher correlation coefficient (R=0.72, Figure 7(d)) than those of univariate regression analyses (R=0.61, Figures 6(h) and S5(l)). The dolomite and clay minerals with VIP values higher than 1 (Figure 8(d) and Table S3) have relatively important influences on D1 (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 Es4U shale is mainly self-structured crystals with the flat surface [75]. Therefore, D1 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 dimensionD2 of the Es4U shale are mainly controlled by multiple geological factors including clay minerals, feldspar, burial depth, lnRo, carbonates, and calcite, while the fractal dimension D1 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.3.1. Geological Factors Controlling Micropore Volume

The log of the micropore volume (lnVmic) was set as a response (dependent variable) in the PLSR analysis. The descriptors (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 lnVmic 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 Vmic 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)).Figure 7
Correlation relationships between experimental and predicted pore property parameters (Vmic, Vmes, APD, and D1) in the PLSR fitting models. Vmic: micropore volume; Vmes: mesopore volume; APD: average pore diameter; D1: fractal dimension D1. 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.
(a)(b)(c)(d)The standardized coefficients of the multivariate regression equation and VIP values (Figure8(a) and Table S3) of the PLSR fitting model of lnVmic 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 Es4U 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 Es4U 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 Es4U shale.Figure 8
Variable importance in projection (VIP) of geological factors for the PLSR fitting models of pore properties.Vmic: micropore volume; Vmes: mesopore volume; Vmic‐mes: micromesopore volume, sum of Vmic and Vmes; S: specific surface area; APD: average pore diameter; D1 and D2: fractal dimensions. The data are listed in Table S3.
(a)(b)(c)(d)


## 4.3.2. Geological Factors Controlling Mesopore Volume and Specific Surface Area

TableS3 also shows the standardized coefficients and VIP values of the descriptors obtained from the PLSR analysis of the log of the mesopore volume (lnVmes). The predicted values of the PLSR fitting model of lnVmes (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 Vmes (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 multivariate regression equation (Figure 8(a) and Table S3) of the PLSR fitting model of lnVmes 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 lnVmic, 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 Es4U 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 Es4U 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 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 (lnS) and micromesopore volume (lnVmic‐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 lnVmic‐mes and lnS are close to those in the lnVmes 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 lnVmic‐mes and lnS (see Section 3.7 for details). These features are the same as those in the PLSR fitting model of lnVmes.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 volumes are controlled or affected by multiple geological factors. The mesopore and micromesopore volumes as well as the specific surface area of the Dongying Es4U 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 Es4U shale need to be synthetically taken into account in the “sweet spot” selection and determination.


## 4.3.3. Geological Factors Controlling Geometrical Parameters of Mesopores

The geometrical parameters of pores including average pore diameter (APD) and fractal dimensionsD1 and D2 were set as the responses for the PLSR analyses. The 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 ofD2 (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 D2 are similar to those controlling APD. But the effects of these geological factors on the response (D2 or APD) are the opposite, as APD has a negative correlation with D2 (Figure 6(c)).Figure8(d) and Table S3 also show the standardized coefficients and VIP values of the PLSR analysis for D1. The relationship between the experimental and predicted values of the PLSR fitting model of D1 has a higher correlation coefficient (R=0.72, Figure 7(d)) than those of univariate regression analyses (R=0.61, Figures 6(h) and S5(l)). The dolomite and clay minerals with VIP values higher than 1 (Figure 8(d) and Table S3) have relatively important influences on D1 (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 Es4U shale is mainly self-structured crystals with the flat surface [75]. Therefore, D1 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 dimensionD2 of the Es4U shale are mainly controlled by multiple geological factors including clay minerals, feldspar, burial depth, lnRo, carbonates, and calcite, while the fractal dimension D1 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 Es4U solvent-extracted samples, absolute pore volumes were directly used to analyze the pore evolution, instead of TOC- or 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.Figure 9
Cross-plots of main pore property parameters with Ro and burial depth for the Es4U 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.
(a)(b)(c)(d)(e)(f)Over the range of 0.6–1.3% Ro, Figures9(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 D1 and D2 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 Es4U shale in the Dongying Depression. The main reasons probably are that the Es4U shale is lacustrine and that organic pores in the Es4U 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 (Figures9(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 Es4U 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 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 (Figures9(a)–9(d) and S6(c)). The high micro- and mesopore volumes as well as the specific surface area of 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 Es4U 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 Es4U shale with the relatively high content of carbonates [79–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 cm3/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 Es4U 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.


## 5. Conclusions

In this study, low-pressure gas (N2 and CO2) adsorption methods were applied to shale samples from the upper part of the Sha-4 Member of the Paleogene Shahejie Formation (Es4U) 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 Es4U shale have an average micropore volume of 0.0038 cm3/g and an average mesopore volume of 0.0182 cm3/g, higher than those of the Soxhlet-extracted samples of Chang-7 shale in the Ordos Basin. In the Es4U shale, the samples with type H3 hysteresis loop make up the most (53.3%) in the studied samples. Furthermore, our samples from the Es4U 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 Es4U shale may potentially be achieved through horizontal drilling and hydraulic fracturing. However, the relationships of micro- and mesopore volumes with brittle minerals in the Es4U 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 Es4U shale oil(3)
Organic pores in the Es4U 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 Es4U 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 cm3/g), specific surface area (18.85 m2/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



---

*Source: 6616843-2021-03-16.xml*