Abstract
Keywords
Introduction
Unconventional shale oil is regarded as a future worldwide energy source and is widely distributed in the major petroliferous basins in China (Zou et al., 2013). Unsuccessful production of shale oil in China has led to a great demand for fundamental investigations of the pore structures of shale oil reservoirs. Pore structures control the storage capacity and transportation properties of hydrocarbons (Cai et al., 2013). Shale oil reservoirs are complex and heterogeneous porous mediums, which are characterized by pore sizes that range widely from the nanoscale to the microscale, as well as rich organic matter and clay minerals, which hinder the understanding of pore structures when using the conventional reservoir characterization method (Chen et al., 2019; Jarvie et al., 2007; Zhang et al., 2017a).
Specialized techniques have been developed for revealing pore networks in shales, for example, mercury intrusion capillary pressure (MICP) (Hu et al., 2017; Zhang et al., 2017b), gas adsorption (GA) (Shao et al., 2017; Wang and Ju, 2015; Zargari et al., 2015), low-field nuclear magnetic resonance (NMR) (Li et al., 2017a, 2017b; Olatinsu et al., 2017), X-ray computed tomography (X-ray CT) (Ma et al., 2016), and scanning electron microscopy (SEM) (Wang et al., 2016; Xu et al., 2018; Yang et al., 2016). The quantitative pore structure parameters (i.e. porosity, permeability, specific surface area, and pore size distribution (PSD)) can be characterized by the methods of MICP, GA, and NMR. However, these quantitative techniques cannot provide direct information about real pore types, morphologies, or the porosity associated with mineralogy and microstructures (Klaver et al., 2012). However, SEM images at high resolution can provide insight into shale microstructures, thus providing information regarding the actual pore type, size, morphology, and distribution, in addition to revealing the heterogeneity of a microstructure. The pore network in shales can be directly studied based on the combination of argon ion-beam milling and SEM (Wang et al., 2016; Yang et al., 2016). Focused ion-beam milling (FIB) in combination with SEM has provided an alternative tool for observing three-dimensional pore networks (Dewers et al., 2012; Zhou et al., 2016). However, these two methods are severely affected by the contradiction between resolution and representation (Klaver et al., 2012). As an enhancement of these methods, broad-ion-beam milling (BIB) combined with SEM allows for the study of a larger planar up to 1–2 mm2 at a high resolution. This method is suitable for qualitative and quantitative investigations of shale microscopic pore structure (Klaver et al., 2012, 2015, 2016; Loucks et al., 2009).
In addition, the heterogeneity of shale pore structure has a controlling effect on storage and transportation capacity (Wang et al., 2016). Fractal theory is an effective method for characterizing pore-structure heterogeneity to reveal the complex pore-networks in shale. The fractal dimension has been widely calculated using nitrogen adsorption data to characterize shale pore-structure (Hu et al., 2016; Shao et al., 2017). However, although the fractal dimension is very suitable for homogeneous rocks, it cannot reveal all of the features of the fractals with heterogeneity and singularity because it characterizes the average properties (Gould and Vadakkan, 2011). Multifractal theory decomposes the self-similar measures into intertwined fractal sets and divides the complex fractals into several regions, which are characterized by the singularity strength and generalized fractal dimensions (Hu et al., 2009). Multifractal analyses can provide more information about the pore-structure properties than the single fractal dimension and has been used to characterize the pore structures of gas shales and tight oil reservoirs based on SEM images or
In the present study, the microscopic pore structure, including pore type, size, morphology, orientation, and PSD, was analyzed in a representative area of a BIB cross-section of the lamellar organic-rich calcareous shale, which is regarded as the most favorable reservoir for shale oil in the Dongying Sag, Bohai Bay basin, China. In addition, the heterogeneity of the major pores and their effects on the occurrence and seepage of shale oil are discussed, which is relevant for providing insight into shale oil reservoir characteristics.
Geological settings
The Bohai Bay basin includes a variety of Mesozoic–Cenozoic fault basins, which were formed and preserved during the process of geological evolution (Figure 1). The Dongying Sag is located in the southeastern area of the Bohai Bay basin and is one of the most petroliferous depressions in China. The sag has an area of ∼5800 km2 and can be divided into four sub-sags: Boxing, Lijin, Minfeng, and Niuzhuang (Figure 1). According to the geological history and sedimentary sequences, the evolution of the Dongying sag can also be divided into syn-rift and post-rift stages (Xie et al., 2006). The main source rocks are dark shales (including mudstone and shale) in the upper part of the fourth Member (Es4U) and the lower part of the third Member (Es3L) of the Paleogene Shahejie Formation, which formed in a saline, humid lacustrine environment during the syn-rift stage. The Es4U and Es3L are the primary target strata for shale oil exploration in the Dongying Sag and contain a high total organic carbon (TOC) content as well as types I and II1 maturity kerogens (0.42–0.64

Sample location (well Y556) and geological setting of the study area (after Zhang et al., 2017a, 2018a, 2019).
Methodology
Sample
Lamellar organic-rich calcareous shale is regarded as the most favorable reservoir for shale oil (Ning et al., 2017). Thus, in this study, the lamellar organic-rich calcareous shale was selected to discuss the effect of microscopic pore structures on the storage and seepage of shale oil. The lamellar organic-rich calcareous shale sample originated from the Es3L member at depth of 2516.16 m in Well Y556. This sample mainly consisted of clay minerals (32.8%), calcite (25.6%), and quartz (19.8%), but contained small amounts of dolomite (6.9%), siderite (5.3%), and pyrite (4.6%). The organic matter characteristics were revealed through the TOC content of 2.62%,
Methods
BIB–SEM imaging
Prior to BIB–SEM imaging, the subsample was pre-polished using silicon carbide sandpaper and was then ion-milled by argon BIB polishing to obtain a planar cross-section of ∼1–2 mm2. The BIB polished cross-section was carbon-coated and imaged using an FEI Helios NanoLab 650 DualBeam SEM with a back-scatter detector (BSE). The BIB cross-section was mapped at a pixel of 10 nm using a 10–20% overlap to create a large area mosaic of ∼1.054 × 0.915 mm2 that comprised 3360 (56 × 60) single images. An energy-dispersive spectroscopy detector was used for semi-quantitative identification of minerals from the BIB–BSE mosaic using a gray scale.
Image processing and determination of REA
The BIB–BSE image gray scale is typically associated to the densities of the materials in shale. The minimum density results in the smallest scale values for pores and is usually represented by black. Pores were segmented and analyzed automatically using ImageJ software, which initially transformed BSE images into an eight-bit bitmap using the tape function. Subsequently, BSE images were converted to pore binary images based on the threshold calculated by the Multi-Otsu thresholding algorithm (Zhang et al., 2017c). The representative elementary area (REA) was determined using the box-counting method on the pore binary images. From these binary images, the surface porosity was calculated within increasing box sizes ranging from 1 × 1 to 60 × 56 in single image. According to fractal theory, the liner relationship of surface porosity versus increasing box size on a log–log plot indicated that the REA was reached.
Moreover, based on the mineral identification in the BSE images, each pore area in a REA pore binary image was manually filled with a designated color using Adobe Photoshop 6.4 software to distinguish the different pore types. Each pore type image was then extracted by the “color threshold” tool in ImageJ. The ImageJ software analyzed these pore type images and provided quantitative pore size (pore area, pore perimeter, Feret diameter, and MinFeret diameter), direction (FeretAngle), and morphological (circularity, convexity, and elongation) parameters for each single segmented pore. The Feret diameter refers to the longest distance between two points in a pore region, while the MinFeret diameter is the minimum pore width, as illustrated in Figure 2. The FeretAngle is the counterclockwise angle between the Feret diameter and horizontal direction (Figure 2). The pore morphological parameters were calculated following the method described in Klaver et al. (2015).

Diagram illustrating of the Feret parameters in a hypothetical pore.
Multifractal method
In the present study, the multifractal was defined based on the generalized dimension. The process for multifractal calculation is as follows.
It was assumed that the SEM pore binary images (fractals) can be segmented in
For multifractal calculation, the partition function (
According to the
Results and discussion
REA determination
The accuracy of image analysis depends on the resolution and scale of an image (Liu and Ostadhassan, 2017). According to the box-counting method, the properties do not vary dramatically when increasing the box size beyond a certain scale (Li et al., 2019). This scale is interpreted as the REA. In this study, point counting was performed by considering the surface porosity (Figure 3(a)). As illustrated in Figure 3(b), the surface porosity varies dramatically with a box size ranging from 1 (1 × 1) to 16 (4 × 4) single images, whereas it shows a gradual decline and does not appear to reach stability when the box size ranges from 25 (5 × 5) to 3300 (60 × 56) single images. Thus, fractal theory was used to determine the REA. For the BIB cross-section, the surface porosity showed a liner relationship with box size on a log–log plot when increasing the box size above 25 (5 × 5) single images (∼100 × 100 µm2), thereby demonstrating that the minimum REA for this sample is 25 single images (∼100 × 100 µm2) (Figure 3(c)).

Principle of REA determination.
Pore characterization of REA
Mineralogy of REA
A typical REA with size of approximately 100 × 100 µm2 (i.e. 25 single images) located in the lower left-hand corner of the BIB cross-section was selected to characterize the pore structure of the shale oil reservoir (Figures 3(a) and 4(a)). The REA corresponded well to the mineral composition as described in Section Sample. The mineral fabric consisted mainly of clay mineral matrix, calcite, quartz, and framboidal pyrite (Figure 4). The calcite and quartz grains were supported by the clay mineral matrix and were generally sub-parallel to the bedding plane. Calcite grains were characterized by plenty of pores as a result of dissolution during diagenesis (Figure 4(d)). Moreover, clay aggregates contained numerous nanometer-sized pores between clay mineral crystals in the clay matrix (Figure 4(c)). The organic matter in the REA was mixed with fine-grained clay minerals and sub-parallel to the bedding plane (Figure 4(e)), which was considered to be solid bitumen (Fishman et al., 2012; Klaver et al., 2015).

SEM analysis of the REA (resolution of 10 nm): (a) interparticle pore, (b) intraparticle clay pore, (c) inter-crystalline pore, dissolution pore; (d) organic pore, (f) interparticle pore, (g) intraparticle clay pore, (h) dissolution pore, (i) inter-crystalline pore, and (j) organic pore.
Pore types and morphologies
According to Loucks et al. (2012), five pore types were encountered: interparticle pores, intraparticle clay pores, dissolution pores, inter-crystalline pores, and organic pores. Interparticle pores predominantly occurred between calcite/quartz grains and the clay mineral matrix and were primarily slit-shaped pores (Figure 4(a), (b), and (f)). The pore area of interparticle pores ranged from 500 to 3,296,700 nm2 (mean 199,642 nm2), and the perimeter varied from 99 to 26,678 nm (mean 3475 nm) (Table 1). The mean Feret and MinFeret diameters were 1164 and 312 nm, respectively. The number of interparticle pores accounted for 3.22% (375/11,638) of the total number of pores, but contributed 17.24% to the total visible surface porosity, thus indicating that interparticle pores were mainly composed of larger pores.
Pore extraction results for five pore types.
Note: Superscripts a, b and c indicate the minimum, maximum and average values, respectively.
Intraparticle clay pores exhibited a complex pore morphology, such as slit-shape to irregular jagged (Figure 4(a), (c), and (g)). The pore size of intraparticle pores was much smaller than the size of interparticle pores, with mean pore area, perimeter, Feret, and MinFeret diameters of 17,659 nm2, 527 nm, 183 nm, and 89 nm, respectively (Table 1). Similarly, intraparticle clay pores had the highest pore number content (49.61%; 5774/11,638), but contributed 30.13% to the total visible surface porosity.
Dissolution pores were mainly developed within calcite grains as a result of diagenetic dissolution and were equidimensional, polygonal, and jagged in shape (Figure 4(a), (d), and (h)). These pores contributed significantly to the total visible surface porosity with an area percentage 50.97% (Table 1). The mean pore size parameters (pore area, perimeter, Feret, and MinFeret diameters) of the dissolution pores were 50,132 nm2, 921 nm, 302 nm, and 159 nm, respectively (Table 1).
Inter-crystalline pores were typical evident in framboidal pyrite aggregates and were mainly polygonal or triangular in shape (Figure 4(a), (d), and (i)). Inter-crystalline pores were rare and did not significantly contribute to the total visible surface porosity. The quantity and area percentage were 0.51 and 1.60%, respectively (Table 1). A total of 12 (0.10%, in number) organic pores were detected from the organic matter and contributed 0.06% to the total visible surface porosity (Figure 4(a), (e), and (j); Table 1), thereby indicating that organic pores are not of great importance to the Dongying Sag shale oil reservoir (Klaver et al., 2012; Li et al., 2015).
Three quantitative morphological parameters (circularity, convexity, and elongation) were applied to further characterize the morphology of four pore types. Organic pores were excluded because the low number of organic pores could not provide representative distributions. Circularity provides an indicator of roundness: the higher the circularity, the closer the pore to a circle. Convexity is a measure of the smoothness of a pore surface: the closer to 1, the smoother the pore surface. Elongation reflects the difference between the pore length and width: when close to 1, elongation indicates a relatively long pore compared to its width (Klaver et al., 2015).
Figure 5 compiles the distribution of the three morphological parameters for all pore types excluding organic pores, because the amount of organic matter pores is too small. Figure 5(a) illustrates that the interparticle pores exhibited circularity with an obvious left-skewed distribution with a frequency peak of between 0.1 and 0.2, whereas dissolution and inter-crystalline pores have comparable distributions of circularity with a weak right-skewed distribution with a peak of between 0.5 and 0.6. The circularity distribution of intraparticle clay pores is characterized by two peaks at 0.3–0.5 and 0.9–0.1. For convexity, Figure 5(b) shows that these four pore types have comparable distributions and present distinct right skewness with the peaks ranging from 0.6–0.7 to 0.7–0.8. Figure 5(c) illustrates that the elongations distribution of interparticle pores are much larger with a peak at 0.8–0.9, while the elongation distributions of intraparticle clay, dissolution, and inter-crystalline pores are similar and show a weak left skewness with a peak at 0.3–0.4. In addition, pore orientation (FeretAngle) distributions of Feret diameters for four pore types were obtained and plotted in Figure 5(d). These four FeretAngle distributions show an advantage angle of 120°–150°, thus indicating a preferential orientation of the longest pore axis sub-parallel to the bedding plane. In general, the pore morphologies of interparticle pores were the most complex and were characterized by the minimal circularity (mean 0.2534) and convexity (mean 0.6307) and the maximum elongation (mean 0.6889). In contrast, the dissolution pore morphologies were the least complicated and were associated with the maximum circularity (mean 0.5487), convexity (mean 0.7630), and the minimal elongation (mean 0.4270).

Frequency plots of pore morphologies and orientations of interparticle (red), intraparticle clay (blue), dissolution (yellow), and inter-crystalline pores (green). (a) Circularity, (b) convexity, (c) elongation, and (d) FeretAngle.
Pore size distributions
The PSDs of five pore types and total pores were calculated using the continuous PSD method proposed by Münch and Holzer (2008), which is also detailed in our previous work (Zhang et al., 2017b, 2018b). In addition, because of the very low number of pores in framboidal pyrite and organic matter (typically <100 pores), the micropore (<100 nm (Zhang et al., 2017b)) contents and average pore diameters (by equation (9)) were only determined for the three major pore types (interparticle, intraparticle, and dissolution pores) and total pores, as exhibited in Figure 6

Pore size distributions of different pore types.
Multifractal characteristics
Based on the SEM pore binary images, the generalized dimension and multifractal spectra of three major pores types were obtained using MATLAB and multifractal geometric principles. In the present study, the range of

(a) Generalized dimensional and (b) multifractal spectra of interparticle (red), intraparticle clay (blue), and dissolution pores (yellow).
Three important parameters (
Multifractal parameters and box-counting dimensions of interparticle, intraparticle clay, and dissolution pores.
The multifractal spectra of three major pore types were calculated and plotted in Figure 7(b). As is observed from this figure, the variation of the widths and the crest of the spectra can be observed, which presents a strong multifractal characteristic. Intraparticle clay pores have the largest
The asymmetry of the singularity spectrum (A) was calculated as (

Homogeneity of interparticle (red), intraparticle clay (blue), and dissolution pores (yellow).
Overall, the fractal characteristics in combination with the SEM image information showed that interparticle pores were mainly slit-shaped pores sub-parallel to the bedding and that they had the largest pore size, the least complex pore structure, and the best connectivity; therefore, they are considered to be crucial to the seepage of shale oil. In contrast, intraparticle clay pores had the smallest pore size, the most complex pore structure, and a lower connectivity; thus, they are interpreted as mainly providing the shale oil occurrence space and as playing a non-important role in the seepage of shale oil in the Dongying Sag. Dissolution pores were found to be the most abundant, but possessed the poorest connectivity with a heterogeneous pore distribution; thus, they are interpreted as being the main occurrence space as opposed to seepage channels for shale oil. Hence, interparticle pores are understood to play a crucial role in shale oil seepage in the Dongying Sag, while intraparticle clay and dissolution pores provide the main shale oil occurrence space.
Conclusions
The REA for the mm-scale BIB cross-section was determined to be approximately 100 × 100 µm2 (25 single images) for the shale oil reservoir sample from Dongying Sag in the Bohai Bay basin. A typical REA from the lower left-hand corner of the BIB cross-section was selected to characterize the pore structure of the shale sample. In combination with multifractal theory, the microstructure and heterogeneity of the shale sample were revealed.
Five pores types were encountered in the typical REA: interparticle, intraparticle clay, dissolution, inter-crystalline, and organic pores. Interparticle, intraparticle clay, and dissolution pores were the major pore types and contributed significantly to the total visible surface porosity with a content of 98.34%, whereas inter-crystalline and organic pores were found to be of lesser importance to the Dongying Sag shale oil reservoir. Of the three major pore types, interparticle pores were characterized by the most complex pore morphologies, the largest average pore diameter (3475 nm), but the least complex pore structure as they had the smallest
Interparticle pores were mainly slit-shaped pores sub-parallel to the bedding plane and had the best connectivity. Intraparticle clay pores showed a lower connectivity and the most complex pore structure. Of the three pore types, dissolution pores were found to be the most abundant but had the poorest connectivity. Hence, we conclude that interparticle pores are crucial to the shale oil seepage in the Dongying Sag, whereas intraparticle clay pores and dissolution pores mainly provide the shale oil occurrence space rather than seepage channels.
