Effect of matching relation of multi-scale, randomly distributed pores on geometric distribution of induced cracks in hydraulic fracturing

Reservoir rock contains many multi-scale, unevenly distributed pores, and the pore structures of shale in different reservoirs and geological environments vary greatly. Because the seepage velocity and pressure field are related to the pore spatial variations, the inhomogeneity of the seepage is superimposed on the anisotropy of the rock’s physical properties, which will affect the distribution of the induced cracks. A method for calculating the pore size in the bonded particle model, based on Delaunay triangulation, is proposed. A modeling approach capable of simulating the multi-scale pore distribution of actual rock is presented based on the proposed method. To understand how microcracks connect micropores in the process of fracturing, several bonded particle model samples with different pore structures were established, and numerical experiments were conducted based on the coupling calculation of the discrete seepage algorithm and discrete element method. The focus of this study was on the interactions between the distribution characteristics of multi-scale pores, the specific physical properties of the fracturing fluid, and the distribution differences of the induced cracks caused by the special seepage characteristics when using different fracturing fluids. The numerical results showed that the advantages of supercritical CO2 fracturing are maximized in deep reservoirs (high in-situ stress) and that a suitable in-situ stress condition is required (i.e. a stress ratio close to 1).


Introduction
The fracturing technique is important for oil production from low-permeability reservoirs. The effect of fracturing stimulation is directly related to the distribution and geometric structure of the induced fracture network. There are many factors that affect the fracture patterns, such as the rock lithology, the mechanical properties of shale rock, the porosity of reservoir, horizontal stress ratio, geological structure, tectonic movement, fracturing fluid viscosity, and fracturing pressurization rate. Many scholars have been devoted to determining the factors that affect the induced fracture network, and have conducted considerable research on this subject. Norman et al. (1963) investigated the effect of an existing fracture or joint plane on the extension of a hydraulically induced fracture through a series of hydraulic fracturing experiments. Mukherjee et al. (2000) studied how the hydraulic fracture geometry evolves when a stress field is disturbed by pore pressure depletion. Haddad et al. modeled hydraulic fracturing through intact reservoir rock using a cohesive zone model and inspected the sensitivity of the fracture geometry and characteristics to the pumping rate, fluid viscosity, leak-off coefficient, and mechanical properties of the reservoir rock (Haddad and Sepehrnoori, 2015). Yuzhang et al. (2015) proposed a parameter called the fracture network growth index to study the multi-factor coupling effect of the geological and engineering factors on the hydraulic fracture network. Over the years, scholars have sought to understand the factors that affect the characteristics of hydraulic fractures and determine what the optimal fracture network is to produce the most effective and economic output.
In most of the previous studies, the reservoir rocks were regarded as isotropic, homogeneous media. In reality, low-permeability reservoirs usually have composite or inhomogeneous structures. The inhomogeneity of a reservoir is very complex and manifests in many aspects, such as nonuniformly distributed natural microcracks, relatively hard gravel randomly distributed in the matrix, anisotropy of the mechanical properties and strength of the rock, and nonuniform pore-size distributions (Chen et al., 2014;Wang et al., 2016b;Yang et al., 2004).
Many types of pores are distributed non-uniformly in shale, such as organic pores, inorganic pores, and microfractures, and the various kinds of pores and fractures usually have different sizes (Chen et al., 2014;Wang et al., 2016b;Yang et al., 2016). The influence of the distributive characteristics of multi-scale pores on the induced fractures is very complicated. First, fluid flow and transport in inhomogeneous porous material are strongly influenced by the interconnectivity of the pores. Second, anisotropic complex pore structures have a significant impact on the mechanical properties and strength of the rock mass. In the process of fracturing stimulation, the coupling between the different seepage characteristics of the fracturing fluid and the isotropic mechanical properties of the rock can lead to very complex behaviors of the fracture propagation.
Homogeneous models cannot reflect the effects of these inhomogeneous structures, and they result in inconsistencies due to the neglect of the underlying physical processes governing the fracture growth and fluid flow in the vastly different pore structures. One reason for ignoring the heterogeneity in the past is that the required number of computations for inhomogeneous models is very large. However, with the advancement of computer technology and the improvement of computing power, scholars have begun to focus on the crucial effect of heterogeneity on hydraulic fracturing. Wang et al. (2016c) investigated how reservoirs respond to variable injection rates in different hydraulic fracturing stages, and they considered the macroscopic characteristics and microscopic heterogeneous characteristics simultaneously by introducing heterogeneity of the failure strength and elastic modulus into their numerical discrete-fracture model. Yang et al. conducted triaxial hydrofracturing tests of heterogeneous rock with irregular, randomly distributed gravel to study the heterogeneity effects of unconventional reservoir rock on the 3D initiation, growth, and distribution of hydrofracturing cracks. Corresponding numerical simulations were also performed. The results showed that material heterogeneity greatly influenced the hydrofracturing behavior Liu et al., 2016). Wang developed a fully coupled non-planar hydraulic-fracture-propagation model that could account for heterogeneous properties (Hanyi, 2015).
However, most of the previous studies considered the heterogeneity of randomly specified strengths or mechanical properties of different lithologies. There have been few studies of the influence of inhomogeneities of the pore structure on the induced fracture. Shale is a fine-grained sedimentary rock, the pore structures of which vary greatly in different reservoirs and geological environments. Some pore structures exhibit unimodal logarithmic normal size distributions, while others exhibit bimodal logarithmic normal size distributions, and the spatial and quantitative distributions of the pore structures are inhomogeneous (Cao et al., 2016;Ju and Wu, 2016;Pan et al., 2017;Wang et al., 2016a;Zhang et al., 2017;Yu et al., 2016). Because the seepage velocity and pressure field are related to the pore spatial variations, the inhomogeneity of the seepage is superimposed on the anisotropy of the physical properties of the rock, and this will affect the distribution of the induced cracks. However, to the best of our knowledge, the dependence of the distribution characteristics of multi-scale pores on the induced cracks is not yet clear.
To understand the mechanism of hydraulic fracturing in reservoir rock with an inhomogeneous multi-scale pore structure, and quantitatively describe the influence of the distribution characteristics of the multi-scale pores on the fracturing performance, a method for calculating the pore size in the bonded particle model (BPM) based on Delaunay triangulation is proposed. A modeling approach capable of simulating the multi-scale pore structure of actual rock and a characterization method describing a matching relation of multi-scale, randomly distributed pores in shale are presented based on this method. Several numerical samples with different quantitative and spatial distributions of the pore structure are established. The discrete seepage algorithm for the BPM was used to simulate the fracturing process. To quantitatively study the distribution characteristics of the induced cracks, the distribution uniformity index (DUI) was used to quantify the complexity of the fracture geometry . This study focused on the distribution of the induced cracks caused by the special seepage characteristics due to the matching relation for the multi-scale pores when using different fracturing fluids. Furthermore, the peculiar characteristics of supercritical CO 2 fracturing due to the interactions between specific physical properties of the fracturing fluid and the inhomogeneous pore structure were explored. The effects of this mechanism under different in-situ stress conditions were also studied.

Bonded particle method (BPM)
The BPM is an efficient tool for simulating the movement and interactions of spherical particles using the discrete element method (DEM) (Potyondy and Cundall, 2004). Brittle solids (such as rock) are modeled as closely packed assemblies of bonded particles in the BPM. The resulting assembly is capable of "fracturing" when the bonds break. It has an innate advantage of simulating the breaking of rock, and it allows the formation and coalescence processes of microcracks to be investigated.
To study the propagation of induced cracks in the process of hydraulic fracturing, several 3D numerical samples with the same shapes and different pore structures were established based on a BPM approach utilizing PFC3D (one sample is shown in Figure 1). These samples were composed of different size particles that were connected using parallel bonds. Particles with different spatial distributions produce different unevenly distributed pores. Comprehensive details regarding the generation of an inhomogeneous pore structure will be described in the following several sections.
Three principal stresses (S v is the vertical principal stress, S h is the minimum horizontal principal stress, and S H is the maximum horizontal principal stress) were applied to the sides of the sample as stress boundary conditions. The boreholes were located in the centers of the samples, and the other dimensions are shown in Figure 1(a). In the current study, a quantitative description of the influence of the distributive characteristics of multi-scale pores on the fracturing performance was investigated. The effect of different in-situ stress conditions on this relationship is also discussed.

Pore-size calculation method in BPM
Because the BPM is a particle-based discrete-element-calculation method, it has excellent simulation abilities for microcrack evolution. However, it does not have built-in tools that can directly simulate the pore structure of rock. Thus, it is difficult to fit the pore distribution of the numerical sample to that of actual rock. We propose a method for calculating the pore size in the BPM. The quantitative distribution and spatial distribution of the pores in the numerical sample were fit to the actual rock based on this method.
Our calculation method for the pore size in the BPM was as follows. First, we created a 3D Delaunay triangulation from the points of every particle center and divided the space occupied by the model into many tetrahedra (a simple schematic diagram is shown in Figure 2(a)). Next, we analyzed each tetrahedron of this triangulation, the volume of solid particles was subtracted from the volume of the tetrahedron, and the remaining void volume represented a pore (as shown in Figure 2(b)).
The shapes of the void volumes in the BPM were very irregular. We used the maximum value of the distance from the Voronoi vertex to the particle center minus the corresponding radius of the particle as the pore size ( Figure 3).
When the radii of the surrounding spheres were similar, the Voronoi vertex was almost at the center of the void, as shown in Figure 3. The radii of the particles were subtracted from the distance between the Voronoi vertex and the four centers of surrounding particles (there are four surrounding particles in 3D; Figure 3 shows 2D examples for clarity). We used the maximum length of the resulting values to characterize the size of the pore. When the radii of the surrounding spheres varied greatly, the Voronoi vertex deviated from the center of the void and approached the largest particle. Thus, the values of the pore sizes could be slightly different because of the sizes of the surrounding particles when the void volume was the same. The above method was always used to calculate the pore size. This did not affect the statistics of the pore size distribution, because the pore-size distribution in our numerical sample was multi-scale, and we counted the distribution of pores at different orders of magnitude.
Finally, the sizes of all the pores were calculated for a specific numerical sample based on the above method. Statistical data of the number of pores at different scales and the characteristics of spatial distribution were obtained. Simulated fracturing experiments were carried out to study the influence of the quantitative distribution and the spatial distribution of the multi-scale pore structure on the induced cracks.

Fluid-solid coupling numerical method
The following is an introduction to the simulation algorithm of fracturing fluid seepage in the pores and the fluid-solid coupling in the BPM. First, a BPM is established, i.e. a compacted particle assembly is established, to represent rock. Next, a 3D Delaunay triangulation is applied, as shown in Figure 4. The void surrounded by four adjacent particles is defined as a "fluid domain" (FD), and the connection between two FDs is defined as a "fluid pipe" (FP), the properties of which are closely related to one side of the tetrahedron.
Each connection (i.e. FP) between two adjacent FDs has a small control interface between three adjacent particles. It controls the flow of fracturing fluid between adjacent FDs, and its shape is complex (as shown in Figure 5).
When the fracturing fluid flows from a FD to its adjacent FD through a control interface, the flow passage can be equivalent to a cylindrical pipe with length (L) and aperture (a). Capillary effects are introduced to the seepage calculations. When FD 2 is full of fracturing  fluid and is about to flow into FD 1 where there is no fracturing fluid, if the pressure difference between two adjacent FDs is less than the capillary force, that is, if Then, there is no flow between the two FDs. In the formula (1), c denotes the interfacial tension, h is contact angle, and a is the equivalent aperture of the flow passage related to the control interface. If the pressure difference between two adjacent FDs is greater than the capillary force, that is, if Then, fracturing fluid will flow from FD2 to FD 1.If FD 1 and FD 2 are filled with fracturing fluid, the capillary effect is not considered. When the fracturing fluid flows between two adjacent FDs, the flow rate (volume of fracturing fluid per unit time) is calculated by the following formula where k is the micropermeability of the FP, ðP 2 À P 1 Þ is the pressure difference between two adjacent FDs, l is the viscosity of the fracturing fluid, and L is the length of the FP, that is, the distance between the centers of two adjacent FDs.
In the simulation of fracturing fluid seepage, each FD receives fluid from the surrounding FPs associated with it. In one time step, assuming that the inflow is positive, the increase in the fluid pressure in the corresponding FD is given by the following equation where K f is the volume modulus of the fracturing fluid and V d is the volume of the FD. The second term in the equation represents the volume change of the FD due to the stress on the BPM. It is assumed that there is a uniformly distributed pressure on the surfaces of the particles surrounded by a polygonal path connecting the contact points constituting the FD ( Figure 6). The force vector applied by the fracturing fluid on the corresponding particles is as follows where n i is the unit normal vector connecting the FD and particle centers and s is the projection area determined by three contact points. The forces from all the FDs in contact with the particle are applied to this particle as volumetric forces, and these volume forces are transmitted to the mechanical DEM calculation process. In this way, the coupling between the deformable fluid and the particles in the BPM is established, as shown in Figure 7.
The following method for deriving the microscopic properties of solids in the DEM can be used to set the relevant microscopic parameters in the seepage calculations, such as the micropermeability. A set of permeability coefficients is assumed, and simulated flow tests are conducted on a model consisting of many FDs and FPs. The permeability coefficients are continuously adjusted until the total permeability of the model matches the actual permeability of the rock, and the relationship between the micropermeability coefficient and macroscopic physical and mechanical parameters, such as the stress and strain, can also be established. This algorithm, which consists of FPs and FD networks, can reproduce quite complex mechanisms, such as Darcy flow and fluid-solid coupling.
The above algorithm accounts for three coupling forms between the fracturing fluid and solid particles: 1. The change in the contact force or the fracture of intergranular bonds will cause a change in the pore size (and the flow resistance caused by it). 2. The displacement of particles will cause a volume change of the FD, which will lead to a change in pressure. 3. The fracturing fluid in the FD exerts a pressure on the particles associated with it.

Generation of numerical sample with multi-scale pores
Shale is formed by sedimentation over a long period of geological time. Its pore structures are varied because of the complex tectonic movement and diagenetic environment. The pore structure not only affects the anisotropy of the mechanical properties of the rock but also influences the flow through the rock matrix. Thus, it can influence the induced cracks in hydraulic fracturing. However, the technical bottleneck of observing the configurations of real pore structures has obstructed understanding of this effect.
In general, there are two methods for generating a pore structure in numerical models. The first is direct modeling. A numerical model is established using computerized tomography (CT) scanning images of the actual rock. The pore structure of this kind model is very close to the real structure, but it is very computationally demanding and requires many CT images to describe the detailed pore geometry (Øren et al., 2007;Zhou and Xiao, 2017;Zhou et al., 2016a).
The second method uses a stochastic pore network model. The pore structure of this kind of model is generated by a random function. It has the advantages of fast generation and high calculation efficiency. However, the detailed pore geometry of the numerical model is not identical to that of the actual rock, and only the macroscopic statistical parameters are consistent (Ju et al., 2014;Li and Zheng, 2015;Zhou et al., 2016b).
In this study, we used a combination of statistical and stochastic methods to create numerical samples of shale rock. First, the statistical data of the distribution of the shale pore structure was investigated. A stochastic method was used to generate a pore structure model with the distribution parameters of the actual rock to study the influence of these parameters on the induced cracks during hydraulic fracturing.
The steps for the formation of a multi-scale pore structure are as follows. First, a skeleton sample model containing large particles with diameters ranging from 0.1 to 2.0 mm was generated. Smaller particles with diameters ranging from 0.02 to 0.1 mm were then inserted to create a dense sample with a sufficiently low porosity. The total number of particles reached about 2.0 Â 10 5 after running the floater-elimination procedure. Following a previously reported method , the mechanical properties of the BPM sample were set, as shown in Table 1. The initial sample was an anisotropic rock model, and the pore distribution was not uniform. However, the pore sizes were not small enough. From the perspective of a multi-scale distribution, the distribution parameters of the pore structure were not consistent with the actual distribution parameters of shale. The second step was to modify the initial model locally. This step was achieved by inserting even smaller particles directly. The voids requiring size modification were chosen using a random function (if only the effect of the quantitative distribution pattern on fracture evolution is studied, the random function is a uniform distribution function in space; if the influence of spatial uniformity of multi-scale pores and its matching relationship on fracture evolution are studied, then the random function should be a non-uniform distribution), as shown in Figure 8. Based on the diagram in Figure 8(a), smaller particles were inserted directly to form smaller pores, as shown in Figure 8(b) (a 2D graph was used here only to illustrate the algorithm, and a 3D model was used in the actual calculations). If the orders of magnitude of the pore sizes are too different, the final BPM will contain a large number of particles, which will lead to a large required storage space and computing time. Another method is to adjust the parameters of the local domain in the seepage calculation (such as the apparent volume multiplier of the domain, aperture multiplier of the pipe, and microscopic permeability), the distribution parameters of the pore structure of actual shale are fit to ensure seepage equivalence. If the microparameters in the seepage calculation for Figure 8(a) are adjusted such that the pressures p 1 , p 2 , and p 3 of the adjacent FD are the same as those in Figure 8(b), the flow between the two FDs is the same, and they are considered to be seepage equivalent. This does not actually insert small particles. In practical calculations, it prevents large storage space and computing time requirements caused by a large number of particles in the BPM. However, when the number of pores is counted, it is counted based on Figure 8 (b). Section "Fluid-solid coupling numerical method" shows the detailed seepage algorithm.  The latter method was adopted in the modeling process. The generated samples had an equivalent porosity range of 3%-5%. Numerical models with two types of pore structures were generated based on the research needs. The first category was used to explore the effect of the quantitative distribution on the induced cracks. The actual pore-size distribution in shale may obey a unimodal lognormal distribution or bimodal lognormal distribution according to our literature survey (Bolton et al., 2000;Cao et al., 2016;Pan et al., 2017). So, it is only in quantity to fit the distribution parameters of the actual shale, the spatial position of making local modification is uniformly random. Second category is used to investigate the effect of spatial distribution. The local position of adjusting the microparameters is no longer a random uniform distribution, and it has spatial unimodal normal distribution or bimodal normal distribution. Then, the mean value and quadratic mean deviation of DUIs of the multi-scale pores are calculated according to the method in Section "Characterization method of matching relation of multi-scale, randomly distributed pores." Some samples that their mean value of DUIs are closer and their quadratic mean deviation of DUIs varies relative larger are selected from many of the obtained samples to carry out fracturing experiments.
To validate the model established above, comparisons of the breakdown pressure between the analytical solutions and the numerical simulation results are shown in Table 2. The analytical solution was calculated using a previously reported theoretical formula (Fjar et al., 2008) To simplify the calculation, it was assumed that the maximum and minimum horizontal stresses were equal. Based on the percentage error, our model can simulate the fracturing process with acceptable accuracy.

Characterization method of matching relation of multi-scale, randomly distributed pores
In shale or other reservoir rocks, various pore sizes are distributed. The number of relatively small and large pores and their different relative location relationships and spatial distributions will have different effects on the induced fracture distribution. In this study, we first analyzed the spatial distributions of the pores in different size ranges based on the DUI method . The DUI is a quantitative index used to describe the uniformity of the fracture distribution. Its value ranges from 0 to 1. The closer the value to 1, the more uniform the distribution. The DUI method can also represent the distribution uniformity of pores in space. In this study, the following algorithm based on the DUI was used to analyze the characteristics of multi-scale, randomly distributed pores. The DUIs of pores of different sizes were subsequently plotted as a histogram. The mean value and quadratic mean deviation of the DUIs were used to characterize the matching relations of different sized pores and the features of the spatial distribution of pores in the rock. The distribution of induced cracks was analyzed based on the fracturing experiment results for numerical samples with different parameters (namely, mean values and quadratic mean deviations). The mean value and quadratic mean deviation of the DUIs in our simulation were calculated using the following method. We divided the pore size into five intervals (for example, the statistical data of a numerical model are shown in Table 3 and Figure 9). The quadratic mean deviation of the DUIs can reflect the uniformity of the pore distribution between different size ranges. If it is zero in an extreme case, the DUI value of the pore distribution of different size ranges is the same, i.e. the matching degree of different sized pores is better. If the quadratic mean deviation of the DUIs is large, the distribution of pores in a specific size range is very non-uniform, and in another specific size range, it is relatively uniform. Because a larger mean value corresponds to a greater sum of the DUIs, the mean value can reflect the uniformity of the distribution of multi-scale pores in general.

Numerical experiment design
Fracturing experiments of samples with different quantitative distributions of pores. To study the effect of the quantitative distributive characteristics of the multi-scale pores on the induced cracks, three numerical samples with different quantitative distributions were established, designated as S-Q-1, S-Q-2, and S-Q-3. The pore distributions of samples S-Q-1 and S-Q-2 were unimodal lognormal distributions, and that of sample S-Q-3 was a bimodal lognormal distribution. Through statistical analysis, the pore structure distribution parameters of three numerical samples were calculated. The mean value of sample S-Q-1 was close to the first mean value of sample S-Q-3, and the mean value of sample S-Q-2 was close to the second mean value of sample S-Q-3. Their pore size distribution diagrams are shown in Figures 10 to 12.
To eliminate the influence of the horizontal stress ratio, all in-situ stresses were set to the same value in the fracturing experiments. The numerical experiments of the above three samples were conducted first under an in-situ stress of 0 MPa for two types of fracturing fluids. The viscosity of the first fracturing fluid was set to 100 mPa.s,, and the interfacial tension was set to 2.5 Â 10 À2 N/m. Its physical properties were similar to those of common fracturing fluids. The viscosity of the second fracturing fluid was set to 0.02 mPa.s,, and the interfacial tension was set to 1.75 Â 10 À4 N/m. Its physical properties were similar to those of supercritical CO 2 . The in-situ stresses were increased to 3 MPa, 6 MPa, and 9 MPa, and two kinds of fracturing experiments were conducted to compare the characteristics of the induced cracks. For simplicity, the variation of the physical properties of the fracturing fluid with the stress state was not considered. Fracturing experiments of samples with different spatial distributions of pores. Using the methods introduced in Section "Generation of numerical sample with multi-scale pores," 12 numerical samples with different spatial distributions of the pore structure were generated by setting different random parameters. The mean value and quadratic mean deviation of the DUI histograms of these samples are shown in Table 4. Due to deposition over a long geological timescale, nanoscale pores dominate fine-grained rock such as shale. In the numerical samples, we constructed, the relatively small pores account for an overall majority, and their spatial distributions were very uniform. Relatively large pores accounted for less, and their spatial distributions were relatively non-uniform, in accordance with the pore structure of the actual rock. Hydraulic fracturing simulations with two kinds of fracturing fluids were performed based on these 12 numerical samples to investigate the influence of the spatial distribution of pores on the induced cracks. All the in-situ stresses were maintained at 1 MPa in this group of fracturing experiments to eliminate the effect of in-situ stresses on the configuration and propagation of the induced fracture network.
Fracturing experiments under high in-situ stress ratio. To eliminate the impact of in-situ stresses on the propagation of induced cracks, all the in-situ stresses were set to the same value in the fracturing experiments described in Sections "Fracturing experiments of samples with  different quantitative distributions of pores" and "Fracturing experiments of samples with different spatial distributions of pores." However, if there was an in-situ stress difference, the significance of the effect of the inhomogeneous pore structure on the induced fracture network was unknown. To explore this, samples whose induced fracture networks were relatively complex were selected from the results of the fracturing experiments described in Sections "Fracturing experiments of samples with different quantitative distributions of pores" and "Fracturing experiments of samples with different spatial distributions of pores" (namely, samples S-Q-2 and S-S-12). Common hydraulic-fracturing experiments and supercritical CO 2 fracturing experiments were conducted separately for samples S-Q-2 and S-S-12 under high in-situ stress ratios (the vertical principal stress S v was 10 MPa, the minimum horizontal principal stress S h was 10 MPa, and the maximum horizontal principal stress S H was 17 MPa).

Results and discussion
To quantitatively study the distribution characteristics of the induced cracks, the DUI was used to quantify the complexity of the fracture geometry. In the distribution diagram of the   induced microcracks from our simulations, some microcracks coalesced into macrocracks, and others were scattered. We did not remove the isolated scattered microcracks when the DUI was calculated, because these scattered distributed microcracks may continue to coalesce into new macrocracks during the fracturing experiment, potentially generating a fracture bifurcation. The DUI not only can reflect the uniformity of the fracture geometry directly and effectively but can also be used to evaluate the fracture bifurcating ability. Thus, we used this quantitative indicator to analyze the fracturing performance.

Influence of pore-distribution characteristics on induced fractures
Influence of quantitative distribution. Although the bonds were broken after the fracturing experiments, the particles remained close together due to the existence of boundary stresses. Thus, the opening of the induced cracks could not be seen from the surface of the sample with the naked eye. In the figures of the results shown below, each small black marker in the diagram actually shows the location of one microcrack, which was recorded sequentially in the numerical simulations, and microcracks coalesced into macrocracks. All the small black markers representing the location of microfractures were projected onto a plane perpendicular to the wellbore to form the two-dimensional figures presented in this paper. Based on the fracture patterns of the numerical results, the macrocrack network was visualized on a cloud map of the microcracks, and the characteristics of the induced fracture network were analyzed. When the in-situ stresses were 0 MPa, 3 MPa, 6 MPa, and 9 MPa, the comparison of the final fracture patterns of the three samples is shown in Figures 13 to 24. Based on the converged macrocracks of sample S-Q-1, there was no evident difference between the results of the hydraulic fracturing with two kinds of fracturing fluids under 0 MPa in-situ stresses. The number of macrocracks from the second hydraulic fracturing process was greater than that of the first hydraulic fracturing process only when the in-situ stresses were relatively high. However, even if there was no evident visual difference in the macrocracks, the DUI of the microcracks from the second hydraulic fracturing process was larger than that from the first hydraulic fracturing process (i.e. there were relatively more microcracks scattered around the surrounding area, and they were evenly distributed). These microcracks also improved the permeability of the rock.
The pore sizes of sample S-Q-1 were smaller than those of sample S-Q-2 overall. The number of visually converged macrocracks after the second hydraulic fracturing process was         greater than that after the first hydraulic fracturing process. The DUI of the microcracks for the second hydraulic fracturing process was also greater than that for the first hydraulic fracturing process. This phenomenon suggests that the advantages of the second hydraulic fracturing process (which represented supercritical CO 2 fracturing) are more evident when the reservoir rock is fine-grained.
To quantify the differences between the induced cracks for the unimodal and bimodal lognormal distributions, the DUIs of the induced cracks in each sample were plotted, as shown in Figure 25.
As shown in Figure 25, the difference between the fracturing performances of the two kinds of hydraulic fracturing fluids was not evident when the in-situ stresses were relatively small, based on the DUI. When the in-situ stresses were relatively high, the difference between the two was clearer. This phenomenon can be explained by the relationship between the breakdown pressure and the in-situ stresses. When the in-situ stresses were small, the breakdown pressure was small. Thus, the seepage velocity in nanoscale pores was small enough to be ignored, and the influence on the induced cracks was very small.  When the characteristics of the pore structure were different, such as the differences between the unimodal and bimodal lognormal distributions, the patterns of the macrocracks were significantly different. The DUI of the microcracks for the bimodal lognormal distribution was relatively high, and that of the unimodal lognormal distribution was relatively small in the first hydraulic fracturing process. In the second hydraulic fracturing process, the DUI was large for the bimodal and unimodal lognormal distributions. This showed that the second fracturing fluid (fracturing fluid with ultra-low viscosity and ultra-low interfacial tension, such as supercritical CO 2 ) had better adaptability to various reservoir rocks with different pore structures.
Influence of spatial distribution. The results for different spatial distributions of the fracture geometry are presented in this section. When the mean value and quadratic mean deviation of the DUIs of pore structure were different, the final fracture pattern varied greatly, as shown in Figures 26 to 37.                The macrocracks of the second hydraulic fracturing process were more complex than those of the first hydraulic fracturing process. The dependence of the DUI of the induced microcracks on the mean value and quadratic mean deviation are shown in Figures 38  and 39. Based on the quantitative DUI analysis, the mean value was smaller when the microcracks were more uniform, whereas the quadratic mean deviation was larger when the microcracks were more uniform. This indicates that the bifurcation of the induced fracture was more significant when the uniformity of the pores distributed in different size ranges varied to a greater extent.

Influence of in-situ stresses
The above results were obtained with no differences in the in-situ stresses. We also conducted fracturing experiments with samples S-Q-2 and S-S-12 under high in-situ stress ratios (the vertical principal stress S v was 10 MPa, the minimum horizontal principal stress S h was 10 MPa, and the maximum horizontal principal stress S H was 17 MPa). We compared the results for the two samples to evaluate how significant the effect of the pore structure on the induced fracture network was. The final fracture patterns of samples S-Q-2 and S-S-12 are shown in Figures 40 and 41.
Based on the results presented above, regardless of the type of hydraulic fracturing performed, the macrocracks were all single symmetric wing-type cracks under high insitu stress ratio conditions. The only difference was that the second hydraulic fracturing process (used to simulate supercritical CO 2 fracturing) created a scattered distribution of microcracks, and there was almost no scattering in the distribution of microcracks around the converged macrocracks in the first hydraulic fracturing process (used to simulate common hydraulic fracturing). The features described in Section "Influence of pore-distribution characteristics on induced fractures" disappeared. This shows that stringent external conditions are required to maximize the advantages of supercritical CO 2 fracturing, such as the reservoir not having experienced large tectonic movement in its geological history and no folds or faults being present. The above characteristics will only occur under the ideal in-situ stress condition (in-situ stress ratio close to 1).

Conclusions
Reservoir rock contains a large number of multi-scale, unevenly distributed pores. Shale from different reservoirs in different geological environments may have different pore structures. The spatial distribution is inhomogeneous and may have unimodal or bimodal lognormal distributions. On the one hand, the inhomogeneity of the pore structure contributes to the anisotropy of the reservoir rock. On the other hand, it will result in the inhomogeneity of the seepage velocity and pressure field. This inhomogeneity of the seepage is superimposed on the anisotropy of rock's physical properties, which will have a certain effect on the distribution of the induced cracks. The effect of the distributive characteristics of multi-scale pores on the induced fracture was studied in this paper. The differences in the fracture distributions of different kinds of hydraulic fracturing processes, which were due to the large differences in the physical and mechanical properties of the fracturing fluid, were analyzed. It was revealed how microcracks connect micropores and coalesce into macrocracks. The results are summarized as follows.
1. When the pore structure of the reservoir rock had a unimodal logarithmic normal distribution, the induced fractures of the two kinds of hydraulic fracturing processes were not significantly different when the in-situ stresses were very small, based on the appearance of converged macrocracks. However, based on the distribution uniformity of microcracks, the DUI of the second hydraulic fracturing process (used to simulate supercritical CO 2 fracturing) was larger than that of common hydraulic fracturing. When the in-situ stresses were large, the differences in the fracture geometry between the two kinds of hydraulic fracturing processes were distinct. The distinction was particularly clear when the mean values of the DUIs of the pores with different scales were very small. 2. When the pore structure of the reservoir rock had a bimodal logarithmic normal distribution, the DUI of the induced microcracks increased with the increase in the in-situ stress in the second hydraulic fracturing process. However, there was no such trend in the first hydraulic fracturing process. This shows that the advantages of supercritical CO 2 fracturing become more significant when the in-situ stresses are large. 3. When the spatial distribution of the pore structure was inhomogeneous, the DUI of the induced microcracks increased with the increase in the quadratic mean deviation of the DUIs of pores with different scales and with the decrease in the mean value. 4. The previous effects were only efficient under an appropriate in-situ stress condition (namely, when the horizontal stress ratio was close to 1.0). To what extent the pore structure can affect the spatial distribution of the fracturing network depends on the in-situ stresses. When the boundary conditions were set to a high in-situ stress ratio (the vertical principal stress S v was 10 MPa, the minimum horizontal principal stress S h was 10 MPa, and the maximum horizontal principal stress S H was 17 MPa), basically all the fracturing results were single symmetric wing fracture types.
Finally, it is important to note that we only discussed the influence of the special physical and mechanical properties of supercritical CO 2 (such as the ultra-low viscosity, high density, and ultra-low surface tension) on the induced cracks. In the actual process of supercritical CO 2 fracturing, the transformation of its state, the thermal stress caused by endothermic and exothermic factors, and the special chemical action of supercritical CO 2 will have various effects on the induced cracks. These interactions were ignored in this study. It is a significant challenge to achieve a thorough understanding of the mechanisms of supercritical CO 2 fracturing.

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 receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the Foundation of Science Project of Beijing Municipal Education Commission (Grant No: SQKM201710016013).