Abstract
Natural fractures are the key factors controlling the enrichment of shale oil. It is of great significance to clarify the distribution of natural fractures to guide the selection of sweet spots for shale oil. Taking the Qing-1 Member shale oil reservoir in the northern Songliao Basin, China as an example, a new method considering the factors affecting fracture distribution was proposed to quantitatively predict the structural fractures. And the effect of natural fractures on shale oil enrichment was discussed. Firstly, the types and characteristics of fractures in shale oil reservoirs are characterized by using core and outcrop data. Combined with the experimental analysis, the influences of fault, mechanical stratigraphy, mineral composition and content, TOC, and overpressure on fracture intensity were clarified. Then, the number and density of fractures are quantitatively predicted according to the power-law distribution of fault length. Next, geomechanical simulation and fracture prediction were carried out on the model which was established with comprehensive consideration of the influencing factors of fracture distribution. Finally, the fracture distribution is evaluated comprehensively based on above prediction. The prediction results in this work are consistent with the core measurements.
1. Introduction
Different scales of natural fractures are developed in shale oil reservoirs [1–6]. These natural fractures are effective storage spaces and main seepage channels for shale oil reservoirs [7–11]. Natural fractures not only control the migration, accumulation, preservation, and single well productivity of shale oil, but also affect the fracturing measures and production results [12–17]. Therefore, the quantitative prediction of fracture distribution has an important guiding significance for shale oil exploration and development.
Many studies have confirmed that there are fractal characteristics between different scales of fractures [18–21]. The number and density of small-scale fractures can be predicted by extrapolating the power law distribution of the accumulated length or the accumulated density of large-scale fractures (faults) identified from seismic data, respectively [22, 23]. However, it is difficult to accurately predict the location of the small-scale fractures using this method. Maerten et al. [24] and Gong et al. [25] proposed a method to determine the number, development location and orientation of fractures based on the combination of the fractal growth model of faults and geomechanical simulation. This method assumes that the surrounding rock is homogeneous and mainly predicts the fractures formed by the disturbance stress field caused by fault activity. However, many studies have confirmed that the shale reservoir has a strong heterogeneity. The development degree of fractures in shale is also affected by the total organic carbon content (TOC), mineral composition, mechanical stratigraphy, abnormal pressure and other factors [26–32]. Therefore, heterogeneous geomechanical model can be established by integrating these influencing factors for fracture prediction [33–35].
Taking the shale oil reservoir of Qing-1 Member in the northern Songliao Basin, China, as a case study, the development characteristics and controlling factors of fractures in the shale oil reservoir were analyzed through fine fracture characterization of outcrops and cores. An evaluation method based on fractal theory and comprehensive geomechanical simulation was established to quantitatively predict natural fractures (mainly for structural fractures). Finally, the influence of natural fractures on shale oil accumulation was discussed.
2. Geological Setting
Songliao Basin is the largest sedimentary basin in Northeast China, with a length of 750 km and a width of 330-370 km (Figure 1(a)) [6, 36]. The tectonic evolution of the northern Songliao Basin underwent three stages, and formed three structural layers, namely, rift structural layer (K1h-K1yc), depression structural layer (K1d-K2n) and inversion structural layer (K2m-Q) (Figure 1(b)) [37]. Large-scale water invasion occurred during the sedimentary period of the First Member of Qingshankou Formation (Qing-1 Member), forming a large lake basin with wide distribution and water depth up to semi-deep lake to deep lake. A set of thick lacustrine dark organic-rich shales were deposited, which were the main source rocks in Songliao Basin (Figure 1(c)). The organic matter type of Qing-1 Member shale in the study area is mainly type I, with TOC of 1% -5%, and the frequency of samples with TOC greater than 2% is about 48%. The Ro of Qing-1 Member is mainly in 0.9% -1.6%, which is in mature-high mature maturity stage [38].

(a)

(b)

(c)
3. Materials and Methods
In this study, strata and fractures were observed and characterized with 2 outcrops (cumulative formation thickness of 120 m) and 315 m long cores. The thickness and lithology of each set of mechanical strata were recorded, and the height, occurrence and spacing of fractures were measured. 72 samples were collected for related experiments to analyze the controlling factors of fractures. TOC tests (24 samples), X-ray diffraction whole-rock mineral analysis (37 samples), and triaxial compression mechanical tests (11 samples) were performed on these samples, respectively. The fault zone structures of two faults from the Niaohe outcrop were characterized in detail, and the width of their damage zones were determined. The influences of mechanical stratigraphy, mineral composition and content, TOC and abnormal pressure on the development intensity of natural fractures were analyzed.
Based on the analysis of the influencing factors of fracture development degree, the distribution of fractures was quantitatively predicted. Firstly, according to the fractal characteristics of different scales of fractures, the number and density of fractures were quantitatively predicted. Although the fracture length-cumulative frequency curve shows a lognormal distribution at a single scale, it was found that the fracture length-cumulative frequency curve shows a good power law distribution through different scales of fractures (Figure 2). Therefore, the number and density of small-scale fractures can be predicted by extrapolating the power-law distribution of cumulative length and cumulative density curves of large-scale fractures (faults), respectively. According to the plane distribution of fault density, the plane distribution of fault-associated fractures was quantitatively predicted.

Then, according to the influencing factors of fracture development, a heterogeneous geological model was established to simulate the palaeotectonic stress field during the fracture formation period. Combined with the actual fracture criterion of rocks in the study area (obtained by triaxial compression mechanics test), whether the rock had reached the failure state and the failure degree were judged. The failure degree was expressed by the failure rate, which is the ratio of the shear or tensile stress to the shear or tensile strength of the rock. By establishing the relationship between failure rate and fracture density of wells, the fracture density between wells were quantitatively predicted. Finally, combined with the plane distribution of fault-associated fractures, a comprehensive evaluation of fracture distribution was carried out.
4. Results and Discussion
4.1. Fracture Types and Characteristics
Three types of fractures, including structural fracture, bedding-parallel fracture and overpressure fracture, were identified in the Qing-1 Member shale oil reservoir (Figure 3), among which structural fractures are dominant (accounting for 77.5%). The fracture density mainly ranges from 0.5 m-1 to 2.0 m-1, with a maximum of 2.52 m-1. According to the mechanical properties and morphology, structural fractures can be divided into extension fracture and shear fracture. Extension fractures are the main ones in the Qing-1 Member shale oil reservoir, accounting for 75% of the structural fractures. They are controlled by mechanical stratigraphy and develop inside the mechanical layer. The extension fractures usually end at the upper and lower interface of mechanical layer. The fracture surfaces are flat and nearly perpendicular to or at large angles (>70°) to the beddings (Figure 3(a)). The fractures in the same set have good equal spacing, and the spacing is proportional to the mechanical layer thickness (Figure 4(a)). Most of the shear fractures are not controlled by the mechanical layer and cut through multiple sets of mechanical layers. The scale of the shear fracture is large, with height from a few meters to dozens of meters. The fracture surfaces are flat and smooth, with obvious striations (Figure 3(b)) and moderate to high dip angles.

(a)

(b)

(c)

(d)

(a)

(b)
Bedding-parallel fractures are near horizontal fractures formed by rupturing along bedding planes during diagenetic evolution under various geological processes, including sedimentation, diagenesis, tectonism, hydrocarbon generation pressurization, dissolution, etc. The bedding-parallel fractures are characterized by small size, small aperture, poor connectivity and strong heterogeneity (Figure 3(c)). The distribution of bedding-parallel fractures is mainly controlled by sedimentation, that is, the primary bedding is the foundation of the development of bedding-parallel fractures. In the low energy fine-grained sedimentary background with rich organic matter, the sedimentary rhythm developed and the sedimentary texture was rich, which lay a good sedimentary geological condition for the formation of bedding-parallel fractures. Bedding-parallel fractures can increase the storage space of shale and improve the seepage capacity of the reservoir to a certain extent, and have an important influence on the effect of hydraulic fracturing in the development stage. The overpressure fractures are short and wide lenticular and appear as veins filled with calcite (Figure 3(d)). Overpressure fractures do not contribute much to shale reservoirs.
4.2. Controlling Factors of Fracture Development
According to relevant statistical analysis, fault, mechanical stratigraphy, mineral composition and content, TOC and overpressure are the main factors controlling the development of fractures in the Qing-1 Member shale oil reservoirs.
Fault is an important factor controlling the development of fractures by controlling the local stress distribution of different structural positions [39, 40]. Figure 5(a) shows two small reverse faults developed on the outcrop in Bin County, Harbin, China. The throw of F1 is 45 cm with occurrence of 103°/65° (dip/dip angle) and the throw of F2 is 25 cm with occurrence of 271°/52° (dip/dip angle). These two small reverse faults have distinct “duality structure”, i.e., the fault zone consists of a fault core and its surrounding damage zones (Figure 5(a)). The fault core of F1 is 3-7 cm thick and consists of fault breccia (Figure 5(b)) and gouge (Figure 5(c)). The fault core of F2 is 1-3 cm thick and is mainly composed of fault breccia. Two types of damage zones were identified in these two faults, namely, linking damage zone and wall damage zone (Figure 5(d)). A series of fractures of different scales are developed in the damage zones. Fractures are densely developed near the fault, especially in the linking damage zone, and the fracture density can reach 23.8 m-1. Fracture density decreases with increasing distance from the fault. When the fracture density is consistent with the regional fracture density, it marks the end of the damage zone (Figure 5(d)). The damage zone widths of faults F1 and F2 are 0.97 m and 0.59 m, respectively, and their maximum fracture densities are 19.6 m-1 and 17.5 m-1, respectively, indicating that the fault scale has an important effect on the maximum fracture density and damage zone width. They are proportional to the fault displacement. When the displacement is high (after several hundred meters), the growth rate of damage zone width decreases.

(a)

(b)

(c)

(d)
The mechanical stratigraphy controlled the fracture nucleation, propagation, and geometry [41–45]. Fractures usually initiated (nucleate) in relatively brittle formations and terminated at the interface between brittle and ductile formations [46]. Fracture spacing statistics in the outcrops confirm that the fracture spacing and fracture density are correlated with the mechanical layer thickness (Figure 4). The fracture density in the mechanical layer depends on the brittleness of the rock, and the higher the brittleness, the more developed the fracture is (Figure 6(a)). The fracture development degree is inversely proportional to rock mechanical parameters such as compressive strength and Poisson’s ratio (Figure 6(b)), and is proportional to the elastic modulus (Figure 6(c)). Under the same stress condition, the shale with high mineral content such as quartz, calcite, dolomite and feldspar have strong brittleness and developed a natural fracture system, while the toughness shale with clay minerals is not conducive to the development of natural fractures. The fractures in the shale with high carbonate mineral content were usually fully filled, while the fractures in the shale with quartz and feldspar were weakly filled. The TOC is positively correlated with the degree of fracture development, that is, the higher TOC, the higher fracture density (Figure 6(d)). Due to the influence of undercompaction, hydrocarbon generation and clay mineral transformation and dehydration, the Qing-1 Member formed abnormally high fluid pressure, which reduced the effective stress and made the rock prone to fracture.

(a)

(b)

(c)

(d)
4.3. Fracture Prediction Based on Fractal Theory
Based on the fault interpretation results of Sun et al. [37], we analyzed the fractal growth model of the faults in Qing-1 Member. Based on the statistics of the fault length in the study area, the fault length-cumulative frequency (Figure 7(a)) and the fault length-cumulative density scatter diagrams (Figure 7(b)) were drawn in the logarithmic coordinate system. A good linear relationship was identified from the data points in the middle of Figure 7(a) (fault length between 2000 m -7000 m). However, data points with fault length less than 2000 m and greater than 7000 m deviate from this linear relationship, which is caused by the low quality of seismic data and some large faults beyond the scope of the study area, respectively [20, 25]. By fitting the relationship between the fault length and the cumulative frequency, the cumulative numbers of different scales of faults were predicted (Table 1). Using the same method, the relationship between fault length and cumulative density was established to predict the cumulative density of different scales of faults (Table 1). The prediction results show that the number of fractures increases rapidly with the decrease of fault scale. There are 195341706 fractures with length greater than 2 m in the study area, and the average fracture density is 1.91 m-1. The predicted fracture density is consistent with the fracture density obtained from the cores, which proves that the fractal theory can be used to predict the fracture density well. The plane distribution of fractures was determined according to the fault density distribution (Figure 8(a)).

(a)

(b)

(a)

(b)
4.4. Fracture Prediction Based on Geomechanical Simulation
Based on the above fracture development characteristics and controlling factors, the distribution map of TOC, formation pressure coefficient and lithofacies during the fracture formation period were selected to establish the geological model. Firstly, the TOC and formation pressure coefficient contour maps and lithofacies distribution map were drawn based on the relevant data of the study area. These maps were then standardized based on the correlation coefficient between them and fracture development degree. At last, these maps were superimposed with the Doublefox software to obtain the geological model.
According to the analysis of the palaeotectonic stress field evolution, fracture formation stages, and the fracture effectiveness, the structural fractures in the Qing-1 Member formed in three stages (sedimentary periods of Yaojia Formation, Nenjiang Formation and Mingshui Formation, respectively) [6, 37]. Among them, most of the structural fractures of the first two stages have been fully filled and are invalid fractures. Therefore, this study is mainly to predict the effective structural fractures of the third stage.
According to the orientation and magnitude of the principal stress obtained from the study of palaeotectonic stress field in the study area [37], the initial boundary conditions were set as follows: the maximum principal stress is N-S direction, and the minimum principal stress is E-W direction; the maximum principal stress is 65 MPa, and the minimum principal stress is 20 MPa. After a lot of inversion calculation, the optimal results were selected to determine the final boundary conditions of numerical simulation. Mechanical model includes the determination of mechanical properties and petrophysical parameters of the geological model. The rock mechanical parameters of different lithologies and lithofacies were determined by triaxial mechanical experiments.
According to the above geological model, rock mechanical parameters, and boundary conditions, the palaeotectonic stress field at the sedimentary period of Mingshui Formation in the study area was numerically simulated using the finite element software ANSYS14.0. The simulation results show that the magnitude of the tectonic stress field is different in different regions (Figure 9). In the northwest of study area, the maximum principal stress is relatively large, mainly distributed in the range of 62-67 MPa, and the minimum principal stress is distributed in the range of 20-22 MPa. Secondly, in the southeast of the study area, the maximum principal stress is mainly distributed in the range of 60-65 MPa and the minimum principal stress is distributed in the range of 18-20 MPa. In other areas, the maximum principal stress is relatively small, generally less than 60 MPa, and the minimum principal stress is 16-20 MPa.

(a)

(b)
Based on the numerical simulation of tectonic stress field, the extension failure rate and shear failure rate of the Qing-1 Member were calculated. The results (Figure 10) show that the shear failure rate of the Qing-1 Member is mainly distributed in the range of 0.5 - 1.5, while the extension failure rate is mainly distributed in the range of 1.2 - 2.5, generally greater than 1, indicating that the study area is dominated by extension fractures, which is consistent with the observed results of outcrops and cores. The plane distribution of fractures was predicted by establishing the relationship between total failure rate and fracture density. The prediction results (Figure 8(b)) show that the fracture density is mainly distributed in the range of 0.5 - 2.0 m-1, which is consistent with the core observation results. Fracture density is highly heterogeneous on the plane, and fracture development is controlled by lithofacies. Fractures are most developed in high TOC laminated siliceous shale, followed by medium TOC laminated siliceous shale, and poorly developed in massive argillaceous shale, which reflects the comprehensive effects of mechanical stratigraphy, TOC, mineral composition and content on fracture density.

(a)

(b)
4.5. Fracture Comprehensive Evaluation
As can be seen from Figures 8(a) and 8(b), the fracture density distributions predicted by the above two methods are greatly different. This is due to the consideration of different factors affecting fracture development by the two methods. Therefore, in order to obtain a more accurate fracture distribution law, both of these two prediction results should be taken into account.
Based on the analysis of fracture genetic types, development characteristics, controlling factors, and quantitative characterization of fracture parameters, combined with the numerical simulation of palaeotectonic stress field and the fracture prediction results of fractal method, a comprehensive prediction of fracture density distribution in the Qing-1 Member was made (Figure 11(a)). A comparative analysis of the fracture prediction results and the fracture density obtained from 15 cores (Figure 12) shows a good positive correlation between them, indicating that the fracture prediction results of this method are reliable. Based on this, a comprehensive evaluation of the fracture distribution was conducted, and the shale reservoirs in the study area were divided into four grades of fractured zones (Figure 11(b)). Grade I fractured zone is the most developed fracture zone with fracture density greater than 2.0 m-1; Grade II fractured zone is a relatively developed fracture zone with fracture density ranging from 1.0 to 2.0 m-1; Grade III fractured zone is a relatively undeveloped fracture zone with fracture density ranging from 0.5 to 1.0 m-1; Grade IV fractured zone is undeveloped fracture zone with fracture density less than 0.5 m-1.

(a)

(b)

4.6. Effect of Natural Fractures on Shale Oil Accumulation
Exploration practices have confirmed that shale oil enrichment is closely related to natural fractures [47]. By comparing the relationship between daily oil productivity and fracture density (Figure 13), the high-productivity wells are located in the Grade I and II fractured zones. However, the relationship between the oil productivity and fracture density is also very complicated. Four phenomena in exploration practices can reflect this complexity. (1)No fracture, no oil enrichment. When the fracture density is less than 0.5 m-1, the oil productivity is low or not at all. Even with fracturing treatment, it is difficult to obtain industrial oil flow(2)Low fracture density, requiring fracturing. In the study area, some wells have low developed natural fractures with density of 0.5-1.0 m-1, but most of the structural fractures and bedding-parallel fractures contain oil. Although the natural oil productivity is low, industrial oil flow can be obtained after fracturing(3)Fractured reservoirs with moderate fracture density. As early as 1999, more than six wells produced industrial oil flow in the shale of Qing-1 Member [16, 47]. These shales were considered fractured reservoirs because they have low porosity and permeability but have acquired natural productivity(4)High fracture density but low productivity. There are two possible reasons for this phenomenon. On the one hand, although fractures are highly developed, most of them are filled with minerals and are ineffective fractures. On the other hand, due to the high degree of fracture development, oil leakage channels were formed, and the oil migrated to other formations along the fracture system. For example, the predicted fracture density of well SY1 is as high as 2.86 m-1, but the productivity was only 1.70 t/d. Through statistical analysis of fracture prediction results and oil productivity, with the increase of fracture density, the daily oil productivity shows an upward trend, but with the further increase of fracture density, the daily oil productivity shows a downward trend (Figure 13), which may be due to the oil leakage caused by over-developed fractures

5. Conclusions
There are three types of fractures in shale oil reservoirs: structural fracture, bedding-parallel fracture and overpressure fracture, among which structural fractures are dominant. Mechanical stratigraphy, mineral composition and content, TOC, abnormal high pressure, and fault are the main factors controlling the development of fractures in shale oil reservoirs. The number and density of fractures can be quantitatively predicted according to the power law distribution of fault length. The plane fracture density distribution can be quantitatively predicted by geomechanical simulation with comprehensive consideration of fracture controlling factors. The accuracy of fracture prediction can be improved by combining these two prediction results.
Shale oil enrichment is closely related to natural fractures. When natural fracture density is low, there is no connectivity between the fractures, making it difficult to achieve natural productivity without fracturing. With the increase of fracture density, fractures in the mechanical layer are connected, but the fracture density is not high enough to cut through the roof and floor, that is, there is no vertical passage, and a fractured shale reservoir with natural productivity is formed. When the fracture density exceeds a certain value, the fractures are connected vertically and oil is transported to the overlying formations. Therefore, high fracture density may be detrimental to shale oil enrichment. Evaluation of subsurface fracture effectiveness and connectivity becomes the “golden key” for shale oil sweet spot selection.
Data Availability
The data involved in this paper are all in the text of the manuscript.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Authors’ Contributions
Lei Gong: Conceptualization, Writing—review & editing; Shuai Gao: Methodology, Original draft preparation; Bo Liu: Formal analysis; Jianguo Yang: Writing—review & editing; Xiaofei Fu: Investigation; Fei Xiao: Writing—review & editing; Xiaocen Su: Data curation; Rongzhi Fu: Visualization; Qi Lu: Investigation.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (grant numbers 42072155, 41902150, and U20A2093), Natural Science Foundation of Heilongjiang Province (grant nos. JJ2020LH0129 and TD2019D001), Young Innovative Talents Training Program for Universities in Heilongjiang Province (grant number UNPYSCT-2020147), China Postdoctoral Science Foundation (grant number 2018M631908), and Fund of Northeast Petroleum University (grant numbers HBHZX202001 and 2018QNL-12). The authors thank Jianguo Huang and Bingyang Lyu at Northeast Petroleum University for their constructive help.