Broad ion beam-scanning electron microscopy pore microstructure and multifractal characterization of shale oil reservoir: A case sample from Dongying Sag, Bohai Bay Basin, China

Pore structure and its heterogeneity are critical factors controlling the storage capacity and transportation properties of hydrocarbons. Broad ion-beam-milling scanning-electron microscopy allows for the study of a larger planar at high resolution than other methods and can provide insight into shale microstructures. In this study, we investigate the microscopic pore structure of a shale oil reservoir sample from Paleogene Shahejie Formation in Dongying Sag, Bohai Bay Basin, based on the broad ion-beam cross-section, and discuss the heterogeneity of the major pores using multifractal theory. The representative elementary area of the sample was first inferred to be ∼100 × 100 µm2 (25 single images) for the broad ion-beam cross-section with an area of 1.054 × 0.915 mm2. Five pore types (interparticle, intraparticle clay, dissolution, inter-crystalline, and organic) were subsequently identified and analyzed in the selected typical representative elementary area. The results showed that interparticle, intraparticle clay, and dissolution pores were the major pore types and made a significant contribution to the total visible surface porosity (98.34%), whereas inter-crystalline and organic pores were not of great importance. Interparticle pores exhibited the most complex pore morphologies, the largest average pore diameter, and the simplest pore structure. Moreover, interparticle pores that were sub-parallel to the bedding plane showed the best connectivity. Intraparticle clay pores, on the other hand, had the smallest average pore diameter, the most complex pore structure, and their distribution in a two-dimensional plane was the most homogeneous. Dissolution pores were characterized by the least complex pore morphologies but more heterogeneous pore distribution. Both intraparticle clay and dissolution pores were abundant but possessed poor connectivity. We conclude that for shale oil storage and transportation in the Dongying Sag, interparticle pores play an important role in shale oil seepage, whereas intraparticle clay and dissolution pores provide the main space for the occurrence of shale oil.


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) Olatinsu et al., 2017), X-ray computed tomography (X-ray CT) (Ma et al., 2016), and scanning electron microscopy (SEM) 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 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 mm 2 at a high resolution. This method is suitable for qualitative and quantitative investigations of shale microscopic pore structure (Klaver et al., 2012(Klaver et al., , 2015(Klaver et al., , 2016Loucks et al., 2009).
In addition, the heterogeneity of shale pore structure has a controlling effect on storage and transportation capacity . 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 T 2 spectral Zhao et al., 2017). Unfortunately, few studies have applied multifractal theory to characterize the pore structure of shale oil reservoirs.
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 km 2 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 (Es 4 U ) and the lower part of the third Member (Es 3 L ) of the Paleogene Shahejie Formation, which formed in a saline, humid lacustrine environment during the syn-rift stage. The Es 4 U and Es 3 L 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 II 1 maturity kerogens (0.42-0.64 R o %) (Zhang et al., 2014).

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 Es 3 L 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%, S 1 concentration of 0.4749 mg/g, and T max of 437 C. The shale porosity and permeability were 7.40% and 0.046 Â 10 À3 mm 2 , respectively. This shale sample was considered to be representative of the shale oil reservoir with respect to organic matter abundance, maturity, mineral composition, porosity, and permeability. A subsample for BIB-SEM analysis was cut perpendicular to the bedding plane under dry conditions to prevent hydration.

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 crosssection of $1-2 mm 2 . 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 mm 2 that comprised 3360 (56 Â 60) single images. An energydispersive 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).
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 N T partitions with scale r, and that the probability of the ith partition (P i (r)) can be expressed as equation (1) (Zhao et al., 2017) where S i (r) is the surface porosity of ith partition. For pores with a multifractal property, the relationship between P i (r) and scale r can be expressed as equation (2) Figure 2. Diagram illustrating of the Feret parameters in a hypothetical pore. (Halsey et al., 1986) where a i represents singularity strength, which characterizes the density in the ith box, and N a (r) denotes the number of boxes for P i (r) with singularity strengths between a and a AE da, which is related to scale r and can be expressed by equation (3) N a r ð Þ / r Àf a ð Þ ( 3) where f(a) is usually called the multifractal spectrum. For multifractal calculation, the partition function (N(q,r)) of q with scale r is defined by equation (4) N q; r ð Þ ¼ where s(q) is the mass exponent and expressed as equation (5) s q ð Þ ¼ Àr lim According to the P i (r) and q, the generalized dimension D q can be defined as equation (6) (Halsey et al., 1986) where D q reflects the overall singularity of each box and is a strictly monotonically decreasing function of q. Subsequently, the f(a) and a can be obtained using the Legendre transformation from the q and s(q), as expressed by equations (7) and (8), respectively (Chhabra and Jesnsen, 1989) q -D q and af(a) are the basic mathematical tool for describing the heterogeneity of pore distribution.

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 . 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 mm 2 ), thereby demonstrating that the minimum REA for this sample is 25 single images ($100 Â 100 mm 2 ) (Figure 3(c)).

Pore characterization of REA
Mineralogy of REA. A typical REA with size of approximately 100 Â 100 mm 2 (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).
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 nm 2 (mean 199,642 nm 2 ), 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. 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 nm 2 , 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 nm 2 , 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 leftskewed 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).
Pore size distributions. The PSDs of five pore types and total pores were calculated using the continuous PSD method proposed by Mu¨nch and Holzer (2008), which is also detailed in our previous work (Zhang et al., , 2018b. In addition, because of the very low number of pores in framboidal pyrite and organic matter (typically <100 pores), the micropore (<100 nm ) 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 D a ¼ where m is the total number of the pores, p i is the proportion, and d i is the diameter of each pore. The results showed that the pores of this shale oil reservoir sample were mainly composed of nanopores with the micropore content and average diameter of 27.80% and 203 nm, respectively. Of the three major pore types, the intraparticle clay pores had the maximum micropore content (38.95%), followed by dissolution (23.94%) and interparticle pores (23.52%), which corresponds to the average diameters (from low to high) of 184, 209, and 221 nm.
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 q was set from -20 to 20 with steps of 0.5. The generalized dimensions (D q ) versus variable q of the three major pore types are shown in Figure 7(a). This shows that all pore types exhibit decreasing D q values with increasing q. However, D q shows a pronounced decrease as q increases when q < 0, while a smaller difference is observed for D q when q > 0, which is similar to the  generalized dimension spectral of tight oil reservoirs based on the T 2 spectral for water saturated conditions (Zhao et al., 2017). Three important parameters (D 0 , D 1 , and D 2 ) that are commonly used for multifractal analysis were determined and are listed in Table 2. D 0 (capacity dimension) indicates the pore structure complexity by providing the average values of the analyzed structure distribution, and D 1 and D 2 are the information dimension and correlation dimension, respectively. The same relationship (i.e. D 0 > D 1 > D 2 ) was observed for all three pore types, thereby revealing the multifractal nature of the pore distributions for these three pore types. Intraparticle clay pores had the highest D 0 value (1.971), whereas interparticle pores had the smallest (D 0 ¼ 1.832), thus indicating that intraparticle and interparticle pores had the most and least complex pore structures, respectively. The same distributions of D 1 and D 2 were also found for the three pore types ( Table 2).
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 f(a) max value corresponding to its largest D 0 value. In addition, the values of a max and a min can be obtained from Figure 7(b), which reflect the fluctuation of the minimum and maximum probability sunsets, respectively (Li et al., 2012). The singularity strength range (Da) refers to the difference between the a max and a min , which is a quantitative measurement of the degree of multifractality: the larger the Da, the stronger the multifractality (Hu et al., 2009). The Da values ranged from 3.148 to 4.086, which are much larger than the Da values of Bakken shales in the Williston Basin in Montana, North Dakota (USA), which is one of the largest unconventional shale oil plays in the world (Liu and Ostadhassan, 2017). This indicates the presence of a strong degree of multifractality.
The asymmetry of the singularity spectrum (A) was calculated as (a 0 -a min )/(a max -a 0 ) (Posadas et al., 2003) for each of the three pore types (Table 2) and were all < 1, which indicates a low exponent and slight fluctuation. Moreover, the magnitude of the difference in values between a 0 and D 0 provides a measure of the heterogeneity of pore distributions (Li et al., 2012). If all of the data points fall on the 45 line, the pore distribution is considered to be homogeneous and can be characterized by D 0 alone. However, Figure 8(a) shows that the data points of all three pore types all deviate from the 45 line, thus indicating that the pore distributions of these three pore types are heterogeneous and can be more appropriately described by the multifractal spectra rather than a mono-fractal model. The larger the distance between a point and the 45 line is, the more heterogeneous the pore distribution is. This demonstrates that the interparticle and dissolution pores are more heterogeneous in a two-dimensional plane. The box-counting dimensions (D b ) were also analyzed based on the binary images for the three main pore types. Figure 8(b) shows that the interparticle and dissolution pores have the larger box-counting dimension values, which is consistent with the more heterogeneous distributions of interparticle and dissolution pores.
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 mm 2 (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 D 0 value (1.832). Conversely, intraparticle clay pores had the smallest average pore diameter (572 nm) and the most complex pore structure (D 0 ¼ 1.971). The least complex pore morphologies were the main features of the dissolution pores. Moreover, larger singularity strength range (Da) and box-counting dimensions (D b ) indicated that interparticle and dissolution pore distributions in a two-dimensional plane were more heterogeneous.
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.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding
The author(s) disclosed the following financial support for the research, authorship, and publication of this article: National Natural Science Foundation (Grants 41972123, 41602131, 41672130, 41772131).