Characterization of Shale Pore Structure by Multitechnique Combination and Multifractal Analysis and Its Significance

School of Human Settlement and Civil Engineering, Xi’an Jiaotong University, Xi’an 710049, China Key Laboratory of Coal Resources Exploration and Comprehensive Utilization, Ministry of Land and Resources, Xi’an 710021, China No. 8 Oil Production Plant, Changqing Oilfield Company, PetroChina, Xi’an 710021, China Hunan Geological Exploration Institute of China Metallurgical Geology Bureau, Changsha 410016, China


Introduction
Shale oil exploration and development enhance the necessity of determining the pore structures [1][2][3][4]. Previous studies have proved that oil-bearing shales generally have tiny pores and narrow throat [5][6][7]. In addition, various techniques output different pore size distributions (PSDs), which made it difficult to evaluate the pore structures precisely [8][9][10][11]. More importantly, the voids in shales could be further divided into inorganic and organic pores, and the types of pores were closely related to the sedimentary and diagenetic environment with completely different surface physical and chemical properties [12,13]. Moreover, the radius of pores varies among different regions, and determining the proportion of each region was of great signifi-cance [12,14]. Thus, the combination of various experiments and the calculation of the percentage of pores with different radius were important to the analysis of the characterization of shale pore structure [15][16][17]. After those steps, the multifractal characteristics, origin of subpore features (for example, surface roughness and heterogeneity) could be investigated from the perspective of relationships among multifractal dimensions, mineral compositions, and pore structures [16]. And these trends were of great importance to probe the pore network distribution rules and evaluate the reservoir qualities.
Since the inorganic and organic pores were found in shales, lots of equipment has been used to study the pore structure features [5,8,12,18,19]. The morphology of shale pores could be visualized via imaging equipment, like field emission scanning electron microscopy and computerized tomography [20][21][22]. Using intrusion tests (for example, pressure-controlled mercury intrusion) and nonintrusion experiments (for instance, nuclear magnetic resonance), the pore network distributions could be quantitatively gained, and the pore structure-related parameters, such as average pore radius, pore volumes, and surface areas can be calculated [23][24][25]. Besides, the pores with different radius and surfaces have unique physical and chemical properties, and previous research has illustrated that sole method leads to the limitation of space recognition [20]. In summary, from different aspects, the heterogeneity of pore structures could be studied by various techniques; however, whether different types of pores improve, deteriorate, or have no impact on pore networks remains unclear.
Multifractal analysis was prevailed due to the significance of heterogeneity of PSDs of porous media determination [26]. Fractal dimension has been considered as an important parameter to understand the complexity of pore structures and has an immense application value [27]. However, this parameter was inadequate to evaluate the irregularity of pore networks; hence, as an extension of fractals, multifractals had widely been used to study the complexity of the porous media via decomposing the multifractal network into intertwined fractal subsets [28]. Through pore radius, morphology and surface roughness affect the reservoir qualities, the degree of pore-throat heterogeneity was of crucial importance, and multifractal dimension could be seen as an effective tool [29,30]. Previous studies on multifractal analysis were commonly focused on tight sandstones and coal [31,32]. In contrast to tight sandstones that generally have relatively large pores and few bitumen, continental shales have experienced deep-water sedimentation and strong diagenesis but few episodes of intensive tectonic movements [33]. Therefore, the appearances of shales, such as mineralogy, pore structure, and hydrocarbon, were totally different from that of sandstones. Characterizing the pore-throat complexity and irregularity of the continental shales was of great significance for the understanding of the oil accumulation in the Ordos Basin.
In this study, the PSDs in a narrow range that could be probed by sole method will be extended well into the full range, to the greatest extent matches the pore network characteristics, using hybrid techniques. Intrusion technique can detect pores that are well-connected and relatively large, but the ink-bottle pores which were abundant in shales and the pores below to 30 nm may be overlooked. Adsorption methods can be used to determine more information on the pore network structures, while inversion algorithms need to be used to reveal details of radius distributions induced by adsorption volumes or relaxation times. Scanning methods were also conducted for the purpose of validation and described the pore-throat networks precisely. Further, the fractal dimensions from various tests and the multifractal information of different types of shales with different pore size ranges were analyzed. These jobs could be applied to determine the heterogeneity of the pore networks and probe the governing factors of irregularity of the pore and subpore structures.

Methodology
2.1. Geologic Settings. The Shanbei Slope is one of the secondary structural units in Ordos Basin (Fig. S1a) [34]. The Yanchang group was deposited in lacustrine with fluctuated lake level, and it could be divided into ten members from bottom to top (Fig. S1b) [35]. The Chang 7 Member (C7), which is located in the low part of the group, mainly deposited lake shale that has thickness more than 80 m intercalated by sand body and mudstone (Figure 1(b)) [36]. The C7 shale was considered as the source rock in previous unconventional petroleum development, possessing high total organic carbon (TOC) and vitrinite reflectance (R o ) values with low thermal maturity, and type II1 was the predominant organic matter [36]. Previous exploitation job revealed that the wells carried out in C7 shales obtained commercial oil flow, suggesting that this member has a bright prospect for shale oil development [36].

Specimens and Experiments.
The total five specimens were collected from five wells with the depth range of 1433.43-1521.21 m. A series of tests were conducted, including field emission-scanning electron microscopy (FE-SEM), X-ray diffraction (XRD), geochemical analysis, Nitrogen (N 2 ) adsorption, pressure-controlled mercury intrusion (PCMI), and nuclear magnetic resonance (NMR).
The specimens were first grinded and mechanically polished before the tests in order to create a flat surface, and then the surfaces were argon ion milled and coated with carbon. After those procedures, the micromorphological and structural features were detected by Zeiss GeminiSEM 500 apparatus (FE-SEM).
For XRD tests, the specimens were first crushed and then mixed with ethanol. Subsequently, each specimen was smear mounted on glass slides for X-ray diffraction analysis with the help of Bruker D8 ADVANCE apparatus at 40 kV and 40 mA. Finally, the mineral compositions and contents were acquired by deciphering the spectra.
For the geochemical analysis, maceral types, total organic carbon (TOC), and vitrinite reflectance (R o ) were tested. Before the test, the specimens were powdered and treated with 10% hydrochloric acid to remove carbonate, and then placed into a liquiTOC II Analyzer to measure the TOC contents. Organic petrography was carried out on polished specimens made with a cold-setting epoxy-resin mixture. Zeiss Axioimager II microscope system equipped with ultraviolet (UV) light source and the Diskus-Fossil system were adopted to analyze the maceral types and R o , and the measurements were within the framework and restricted by the yttriumaluminum-garnet standard reference under oil immersion. N 2 adsorption was a significant method for the analysis of the pore structures, and the determination of PSDs ranged from 0.30 to 300 nm by Micromeritics ASAP 2020 Plus HD88 analyzer. Before the tests, volatile matter and moisture were removed by the way of degassing, and then all specimens were exposed to N 2 at 77 K. The surface area and pore volume were calculated by the five-point Brunauer, Emmett, and Teller (BET) method and the Barrett, Johner, and Halenda (BJH) method.

Geofluids
For PCMI tests, the specimens were first treated in high temperatures at 110°C over 24 h to dry the moisture. Then, all specimens were evacuated to 0.001 psi, and mercury was intruded into the blocks from 0.03 to 200 MPa. Micromeritics Autopore 9400 porosimetry was used to determine the equivalent pore radius with the help of Washburn equation (1921) [37] and the findings from Purcell [38].
Transverse relaxation time (T 2 ) from NMR apparatus was made by Niumag Analytical Instrument Corporation (23 MHz). The specimens were saturated with brine and then were putted into the apparatus, and the echo time interval, signal superposition time, and echo numbers were set as 0.07 ms, 64, and 1000, respectively. After the detection in saturated status, the core plugs were centrifuged and then detected the spectra in the same experimental parameters.

Multifractal Parameters.
In order to execute multifractal calculation in shales, the box-counting method was applied to the data to determine the multifractal parameters in our research [39,40]. The boxes with equal length l were conducted for the purpose of data record, and the boxes were labeled by index I where NðlÞ represented the total number of boxes [41]. For different tests and outputted data, the box length had different meanings, for example, the variations of T 2 relaxation time were taken as the length l in NMR tests, while the relative pressure could be seen as the size of the box for N 2 adsorption [42]. Here, we take the variations of T 2 relaxation time as an instance; therefore, the probability mass function for the i th box could be defined as [41,43] follows: where N i ðlÞ represents the volume of the saturated brine for the ith box and N T corresponded to the total volume of brine that was saturated in the pores and throats. When the interval of size was very small, P i ðlÞ could be written as follows: where α i was the singularity exponent. For the porous media which had multifractal properties, N i ðlÞ increased when i decreased and observed a power law: where NαðlÞ represented for the number of the boxes between differential of α in multifractals; f ðαÞ stood for the multifractal singularity spectra, and those parameters could be calculated by the following equations: where q was the exponent expressing the multifractal dimensions, varied from -10 to 10 in increments of 1. This parameter served as a scanning tool to scrutinize the denser or sparser regions of measure μ i ðq, lÞ, and it followed a power law function of l that was defined as follows: where τ q represented the mass scaling function, which could be written as follows:

Geofluids
Then, the generalized multifractal dimensions could be expressed as follows: 3. Results

Mineral Compositions and
Petrography. Clay mineral, quartz, feldspar, and carbonate were major mineral constituents of the C7 shales in the Ordos Basin ( Figure 1). Illite and I/S mixed layer were the dominant clay minerals that formed the total clay compositions ( Figure 1). Based on the three-end diagram of the compositions of different minerals, it was demonstrated that the C7 shales were mainly classified into two groups: brittle mineral high shales (BMHS) and clay mineral high shales (CMHS). Quartz, feldspar, and pyrite were the main brittle minerals in C7 shales, while others were plastic minerals. Apart from sample 3, the plastic minerals play dominant role in C7 shales, indicating that the shales in our research area were easy to transform. Table S1 showed the detailed mineral compositions for those samples. According to the SEM results, the pores in the C7 shales were mostly at the nanoscale ( Figure 2). Pore types existing in the C7 shales include a mix of primary interparticle pores, micropores within clay minerals, and organic matter pores that were created by kerogen to petroleum conversion [44]. The interparticle pores, which were the predominant pore types for BMHS, were located between particle and crystals, while for the CMHS, micropores and organic matter pores were more important. The coexistence of different types of pores resulted in the complex relationships among physical properties, pore structures, and fluid movability.

Geochemistry Parameters.
The maceral types of all the selected shale samples showed distinct variation in the different groups. Although exinite plays a significant role in the maceral and all samples belong to type II1 according to previous research, there were still few differences among different samples, and kerogen index (KI) could be used for the purpose of maceral type determination [45]: where Sap is sapropelinite, Exi is exinite, Vit is vitrinite, and Ine is inertinite. Table S2 showed that samples with high brittle mineral content have smaller KI values than CMHS. For the TOC values, CMHS have the greatest TOC content which ranged from 2.40% to 5.91%, with an average of 4.71% and a TOC content of greater than 5% accounting for over 60%, while the TOC of BMHS ranging from 3.67% to 4.12%, with an average of 3.90%. There was no statistical significant difference between two kinds of shales (~0.90%).
3.3. Nitrogen Adsorption. The isotherm curves for N 2 adsorption and desorption were obtained, and all the samples were belong to types III and IV according to IUPAC (Figure 3) [46]. The capillary hysteresis loops were distinct and similar to those of H3 types, suggesting that slit-shaped pores were the main pore types in C7 shale, corresponding to the SEM images [47][48][49].
3.4. PCMI. PCMI tests were the commonly used method for characterizing pore size distributions, and the pore volume measured by nonwetting phase (mercury) was equal to the volume normalized by the pore volume; therefore, the information of pore networks could be determined by figuring out the relationships between mercury saturation and intrusion pressure [50][51][52]. The views of threshold pressure, median pressure, maximum intrusion saturation, and other parameters Figure 4 showed that the BMHS have relatively uniform distributed capillary curve distributions, while that of CMHS vary a lot with strong heterogeneous pore volume distributions, and different samples have distinct boundaries for the latter group.
3.5. NMR. The NMR T 2 curves were expressed by the pore size distributions where each signal amplitude was represented by its volume [53,54]. In this research, the T 2 spectrum before and after centrifugation was investigated, and the curves were presented in Figure 5. All samples have a distinct left peak and a minor right peak; the former one represented matrix pores while the latter one corresponded to microcracks and larger pores. The amplitudes of the samples almost keep unchanged, revealing that the movable fluid percentages of the C7 shales were very low, no matter what shale types they were in. The T 2 relaxation time of water or brine in the porous media could be described as follows: where T 2 is the relaxation time, F s refers to the pore shape factor, ρ 2 stands for relaxivity, and r represents pore radius. F s and ρ 2 are constant, and then, we could define [55][56][57] Thus, Eq. (11) could be written as follows: Therefore, the NMR T 2 curves could be converted into the pore size distributions when C and n are determined, and the methods are listed below ( Figure 6)  (1) First, we calculated the cumulative amplitude percentage of the curves derived from NMR, LTA, and PCMI (2) Since sole method could only represent the pore size distributions in a short range, the distributions derived from those two hybrid methods have some differences. As shown in Figure 6, any pore-throat radius rðiÞ has a corresponding cumulative pore size percentage CFðiÞ.
(3) In the end, with the help of spline interpolation methods, the cumulative T 2 curve amplitude could be corresponding to a single CFðiÞ, and the NMR spectrum is converted to pore size distributions 4.1.2. Pore Size Distributions Derived from Two Hybrid Methods. The pore size distributions derived from NMR-LTA and NMR-PCMI data are shown in Figure 7 and Fig.  S2. Comparing the pore size distribution curves of two differ-ent hybrid methods, the pore networks of C7 shale could be further investigated. The similar variation trends between NMR and LTA can be seen in Figure 7(a), indicating the feasibility of revealing the pore size distributions using the NMR-LTA hybrid methods. However, as for the NMR-PCMI hybrid methods, the pore peaks of those tests have some differences, suggesting that the pore networks that detected intrusion and nonintrusion methods varied a lot (Figure 7(b)). The detection methods are considered as the reason for these phenomena. The NMR and LTA belong to nonintrusion methods; adsorption is the way of fluid entry of the pore space; however, as for the PCMI, mercury enters the pores due to high injection pressure. Therefore, the pore size distributions detected by those two hybrid methods were totally different. Besides, the parameters that derived from hybrid methods have distinct features. Figure 7 and Fig. S2 show that the BMHS have narrow brine-nitrogen gap area while these gaps in the CMHS are relatively wide. Low scope in large pore of brine adsorption as well as high small and  Table S5.

Multifractal Characteristics from Different Methods
with Same Shale Type. The multifractal features from LTA, PCMI, and NMR methods are obtained, and the generalized dimension D q spectra are depicted in Figure 8. The figure shows that D q spectra demonstrate a monotonic decrease with the growth of q. The parameters that related with generalized dimension are listed in Tables S3, S4, and S5. Here, only some important multifractal dimension-related parameters are presented for the purpose of pore network analysis.
D 0 , which indicates the nonempty boxes with density probability, is the box dimension, and when we segmented the LTA spectra, we made each box has no zero data, meaning that D 0 will be equal to 1 [33,58,59]. D 1 , considered as information dimension, which represents the local scaling property of PSDs [33,59,60]. D −10 and D 10 do not have certain definition yet; however, they can characterize the heterogeneity of the PSDs over the entire pore size range. Apart from those sole parameters, some hybrid parameters are also

Geofluids
included to indicate the pore network features, and those parameters were usually been used for the propose of reservoir feature determination. H, which is ðD 2 + 1Þ/2, is defined as the Hurst exponent, which is generally been used to represent the degree of the positive autocorrelation; with the increase of H, the autocorrelation in distribution of porosity increases [33,61]. D 0 -D 1 can be used as an indicator to reflect the concentration degree of PSDs; high D 0 -D 1 values are corresponding to more clustered style of PSD, that is to say, less homogeneous. D −10 -D 10 is another important parameter to characterize the heterogeneity of the PSDs, with the increase of D −10 -D 10 , the heterogeneity of the objects increases. When we compare the characteristics of the generalized dimension curves derived from different methods, we found that they all follow self-similarity. However, the NMRderived curves are flat when compared to other two method-derived curves (Figure 8). This trend indicates that

Geofluids
NMR method generally detects the voids with good connectivity, and tiny pores (such as intercrystal pores) are hard to be indicated due to the limitations of water molecular size and adsorption ability. Besides, the D 0 -D 1 and D −10 -D 10 from NMR data are significantly low when compared to LTA and PCMI methods, also suggesting that narrow range of pore spaces could by detected by NMR, and it means that for the propose of PSD determination, hybrid methods are necessary due to the limited detection range for different apparatus (Figure 9(a)).

Multifractal Characteristics from Different Shale Types
with Same Method. The multifractal parameters from different methods are calculated and are displayed in Figure 9(b). The multifractal dimensions derived from those methods show that the shale which is abundant in clay minerals has more complex pore network (high multifractal dimension). These phenomena may result from the occurrence of the clay mineral-related pores with wide range of pore radius (Figure 2(e)). Besides, the gap of those two types of shales for LTA method is relatively large, indicating that the pore  Figure 7: Comparison between pore size distributions converted from NMR-LTA and NMR-PCMI measurements (sample 1). 8 Geofluids network detection method which is based on the adsorption ability is more sensitive to the multifractal dimension variation.
The correlations between multifractal parameters (Hurst exponent) from LTA, PCMI, and NMR could reflect the inconsistence of the detecting range from the three  Figure 10). The figure shows that there are not any distinct correlations between H from PCMI and LTA (NMR), while there are relatively distinct correlations between H from LTA and NMR. These trends can be interpreted as that intrusion and nonintrusion methods are totally different mechanisms on detecting the pore size distributions, and the appearance of microfractures during the process of mercury injection may cause the inhomogeneous of pore size distributions. Hence, a combination of the two techniques to detect the pore size distributions may not properly characterize its heterogeneity, as we can see from Figures 10(a) and 10(c).

The Governing Factors of Pore Structure Heterogeneity.
That has been proven that the pore structure heterogeneity has various impact factors, including geochemical features, mineralogy, and pore network characteristics. However, how those factors affect the multifractal dimension and the heterogeneous differences between BMHS and CMHS still be a myth. In this work, the impacts of geological parameters on multifractal dimension of different types of shales are investigated.
4.3.1. The Impact of Mineralogy. The C7 shale belongs to the mixed deposits with various kinds of minerals, including quartz, clay, and carbonate (Table S1). The relationships between different mineral contents and the multifractal dimension are acquired to study what roles various kinds of minerals play in the pore network inhomogeneity. The results of the impact from the mineral compositions on pore network heterogeneity are displaced in Figure 11. Not all minerals have influence on the multifractal dimensions, and the data derived from PCMI generally show opposite trends when compared with other two methods. Those findings also prove that the intrusion and nonintrusion methods have totally different detection mechanisms, and high injection pressure, which derived from PCMI, could lead to the microcrack, and made the pore size distributions from PCMI became unreliable. The results also demonstrate that the development of carbonate and dolomite would lead to high pore network heterogeneity; while with the increase of detrital grains (like buddingtonite) and chlorite, the pore structure becomes homogeneous. However, the intercrystalline pores which are contributed by illite and I/S mixed layer made the void space becomes less heterogeneous (Figure 12(f)).

4.3.2.
The Impact of Geochemical Parameters. Due to the abundant organic matters in shale reservoirs, the geochemical features, such as organic matter types, total organic carbon, and organic maturity, generally play a dominant role in the pore structure evolution [59]. The cross-plots between the multifractal dimensions and geochemical parameters are demonstrated in Figure 11. For the D 0 -D 1 derived from   12 Geofluids

Geofluids
PCMI, which represent the pore structure heterogeneity in a wide range, show that R o and TOC have distinct negative with it. These phenomena reveal that with the increase of maturity, the shale becomes more homogeneous. That is due to the development of plastic component, which would transform but not crack when mercury enters the pores, and lead to less microcrack, which would increase the heterogeneity because of the wide distributed tiny matrix pores. However, when it comes to the heterogeneity in a wide range, NMR became an important tool (Figures 11(d) and 11(e)). This is because that NMR method could detect the full range pores. Kerogen index is the most significant parameter here, with the increase of this parameter, the porous media become more heterogeneity (Figures 13(c) and 13(f)). In summary, the type of kerogen does more impact on the shale pore network heterogeneity.

4.3.
3. The Impact of Pore Structure Parameters. Figure 13 demonstrates the correlation among D 0 -D 1 , D −10 -D 10 , and different kinds of pores. As to D 0 -D 1 , an increasing averaged small pore radius will enhance homogeneity of pore networks, while an increasing averaged middle and large pore radius will improve heterogeneity of pore networks. However, for the D −10 -D 10 , a decreasing averaged middle and large pore radius will decrease heterogeneity of pore structures. Those opposite trends demonstrate that the heterogeneity of pore networks in different pore size ranges may totally different. Besides, only LTA in Figures 13(a)-13(c) and PCMI in Figures 13(d)-13(f) have distinct relationships, reflecting that tiny and large pores have significant role in the variations of the heterogeneity of narrow and wide ranges, respectively.

Conclusions
Through a series of tests and theory, characterization of shale pore structure of the C7 shales in Ordos Basin is investigated. The following conclusions could be drawn: (1) The combination of LTA and NMR could properly characterize almost full pore size distributions in C7 shales, and pores which radius between 1 to 100 nm predominate the pore system (2) Multifractal theory provides an effective solution to determine the pore structure heterogeneity, and the pore network detection method which is based on the adsorption ability is more sensitive to the multifractal dimension variation. The Hurst exponent shows good consistency only in LTA/NMR group, indicating that intrusion and nonintrusion methods are totally different mechanisms on detecting the pore size distributions, and the appearance of microfractures during the process of mercury injection may cause the inhomogeneous of pore size distributions (3) Buddingtonite and carbonate are the fundamental governing mineral factors on heterogeneity of pore structures, while chlorite and illite have obvious impact on the pore networks in narrow and wide ranges, respectively. In geochemical parameters, the type of kerogen does more impact on the shale pore network heterogeneity. A decreasing averaged small pore radius enhances the heterogeneity of pore structures in narrow ranges, while a decreasing averaged small pore radius has the same trends in wide ranges

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