The Investigation on Initiation and Propagation of Hydraulic Fractures in Shale Reservoir

Hydraulic fracturing is a necessary technique for shale gas exploitation. In order to have efficient stimulation treatment, a complex fracture network has to be developed, whereas with rich bedding planes and natural fractures, the mechanism of forming a fracture network is not fully understood and it is so tricky to predict propagation and initiation of hydraulic fracture. Therefore, in this paper, considering the strong anisotropy of shale reservoir, numerical simulation has been conducted to analyze fracture propagation and initiation on the basis of finite element and damage mechanics. Simulation results indicate that hydraulic fracture is not merely controlled by in situ stress due to strong anisotropy in shale. With plenty of bedding planes, hydraulic fracture tends to have initiation and propagation along the bedding plane. In particular, this influence becomes stronger with low strength and high development density of bedding planes. Additionally, in combination with natural fracture and bedding plane, the initiation point is usually on a natural fracture plane, causing relatively small breakdown pressure. In the process of fracture propagation, hydraulic fracture connects with natural fractures and bedding planes, forming dendritic bifurcation and more complicated paths. Numerical simulation proves that bedding plane and natural fracture are vital factors of hydraulic fracture. Compared to natural fracture, the bedding plane has a stronger impact on hydraulic fracture propagation. For the initiation of hydraulic fracture, natural fracture is the major effecting factor. The outcome of this study is able to offer theoretical guidance for hydraulic fracturing in shale.


Introduction
Shale gas is potential unconventional gas in the world. The success of shale gas development in America makes shale gas become a hotspot in petroleum engineering [1,2]. The porosity and permeability of shale are so low that they must be stimulated before having economical production [3]. Hence, in recent years, hydraulic fracturing has been applied to enhance shale gas production. Hydraulic fractures have different morphologies in various geological and engineering conditions. To have better improvement of shale gas production, a complex fracture network is favorable to conduct efficient stimulation treatment according to the field experi-ence of Barnett shale [4,5]. For having a complicated hydraulic fracture network, it is essential to have a deep understanding of hydraulic fracture initiation and propagation. In the beginning, an analytical method has been applied for hydraulic fracture initiation. Hubbert and Willis [6] firstly gave an analytical model for fracture initiation in impermeability formation. Then, Haimson and Fairhurst [7] established an analytical model of fracture initiation for vertical open borehole. For hydraulic propagation, analytical models, which have been built to depict hydraulic fracture, include the KGD model [8], PKN model [9], and P3D model [10]. Based on the above analytical methods, numerous studies on the mechanism of hydraulic fracture initiation and propagation have been conducted, analyzing influence factors such as perforation parameters, rock strength, and fracturing fluid [11][12][13].
To further investigate the hydraulic fracture in shale including multiple discontinuities, numerical simulation has become a major tool. In particular, finite element method (FEM), extended finite element method (XFEM), boundary element method (BEM), and discrete element method (DEM) are the most widely used in hydraulic fracture [14][15][16]. With a long history in numerical simulation, FEM has been well developed and is able to deal with hydraulic fracturing in an anisotropy medium with complex stress [17]. However, in FEM, remeshing is necessary for simulating fracture propagation, lowering the computational efficiency [18]. To overcome this limitation, XFEM removes the remeshing by isolating mesh and fracture, highly improving the computational efficiency [19]. BEM has been established for discontinuous interface in hydraulic fracture initiation and propagation. Since BEM merely has meshing on boundary or discontinuous interface, the number of computing elements decreases, saving computing resources [20,21]. The biggest limitation of BEM is that this method hugely relied on analytical solution. This could cause an obvious error when singular point occurs at the boundary variables [22]. Based on the assumption that rock is the assembly of particles and the continuity is not essential between particles, DEM has a better ability of simulating hydraulic fracture in formation with multiple discontinuities compared to FEM, XFEM, and BEM. However, for the DEM calculation, it is not easy to accurately determine input parameters in and between small particles.
Based on these numerical approaches, many scholars have investigated hydraulic fracture. Guo et al. [23] analyzed the influence of in situ stress state on natural fracture and hydraulic fracture by using a cohesive-zone finite-element model, concluding that small horizontal stress difference is beneficial for hydraulic fracture extending along natural fracture. Tang et al. [24] used three-dimensional displacement discontinuity method to develop a 3D fracture model and analyzed the interaction between hydraulic fracture and discontinuity, indicating that it is easier to have shear failure along structural plane with smaller elastic modulus, and Poisson has no obvious impact. Based on the discrete element method, Lee et al. [25] simulated the shear and tensile failure when hydraulic fracture contacts with natural fractures, concluding that tensile failure likely occurs with low approach angle. Keshavarzi and Mohammadi [26] numerically simulated hydraulic fracture propagation and found that the hydraulic fracture will penetrate through natural fracture at high angle of approach. Based on displace discontinuity method, Olson and Taleghani [27] simulated the multistage horizontal well treatment and discussed multiple fracture propagation. Suo et al. [28] established a coupled mixed-mode model based on the extended finite element method, discussing influences of elastic parameters, in situ stress, and fracturing fluid rate on hydraulic fracture geometry. A seepage-stress-damage model, established by Zheng et al. [29], shows that vertical stress difference is a key factor that determines whether the bedding plane can be initiated. By coupling the parallel bond model and smooth joint model, Chong et al. [30] conducted the research on the effect of injection rate, initial aperture, and permeability coefficient with variable anisotropy degrees. According to a 3D simulation of hydraulic fracture in layered formation from Tong et al. [31], a high-density bedding plane is able to enhance the complexity of the hydraulic network, but the height of hydraulic fracture is dramatically restricted by the high-density bedding plane. Except for numerical model, laboratory tests, such as true triaxial test, CT scanning, and acoustic emission, are applied to observe hydraulic fracture morphology and analyze the mechanism of fracture propagation. Based on the true triaxial test, Olson et al. [32] investigated the interaction between hydraulic fractures and natural fractures, indicating 3 forms of interaction between hydraulic fracture and natural fracture, which are hydraulic fracture bypassing, hydraulic fracture diverting, and the combination of bypass and diversion. Besides, Blanton [33] and Teufel and Clark [34] all found that the interaction between hydraulic fractures and natural fractures is associated with stress state and approach angle and pointed out that the weak plane with small friction coefficient or low cohesion is likely to modify the propagation direction. These conclusions are also proved by Zhou et al. [35] who analyzed the impact of weak plane strength on hydraulic fracture by building layered medium with paper sheets and cement block. Zhou et al. [36] and Jiang et al. [37] applied CT scan and true triaxial test to describe the morphology of hydraulic fractures in shale, showing that hydraulic fracture geometry is largely governed by in situ stress and natural fractures. In addition, combining specimen splitting, CT scanning, and acoustic emission (AE), Ma et al. [38] discussed the geometry of hydraulic fractures and influences of bedding plane on fracture initiation and propagation. Bo et al. [39] observed the micro fracture process of shale under the three-point bending load with SEM, showing that the branch cracks have the tendency of connecting with the main crack. Also, its propagation length is relied on crack propagation resistance and propagation energy. Three-point bending load is also used by Heng et al. [40] to analyze the influence of bedding plane on shale fracture toughness, showing that low consolidation of bedding plane causes small fracture toughness and hydraulic fracture tends to have bifurcation or diversion into bedding plane. Based on large-scale true triaxial test, Tan et al. [41] gave the classification of hydraulic fracture geometries, which are simple fracture, fishbone-like fracture, fishbone-like fracture with fissure opening, and multilateral fishbone-like fracture network. By combining AE event monitoring and CT scanning, Hou et al. [42] found that large horizontal stress difference is adverse to having a complex fracture network because this high stress difference is able to stop fracture diversion into natural fracture or bedding plane, likely forming a major fracture instead of a complex fracture network.
Even if there are so many studies on hydraulic fracture propagation in shale reservoir, it is still a challenging work and its mechanism still has uncertainty. The reason why the mechanism of hydraulic fracture propagation in shale 2 Geofluids is so complicated is strong anisotropy. Shale anisotropy comes from 3 aspects: bedding planes, natural fractures, and heterogeneity of rock elements. In particular, bedding planes and natural fractures ( Figure 1) have been considered as one of the vital influencing factors of hydraulic fracture. Influences of bedding plane or natural fracture on hydraulic fracture have been extensively studied, but the work of coupling bedding plane, natural fracture, and heterogeneity of rock elements is rarely conducted. Therefore, in this paper, considering the anisotropy of bedding plane, natural fracture, and heterogeneity of rock elements, RFPA has been used to establish a numerical model of hydraulic fracture by setting different mechanical properties for discontinuities (bedding plane and natural fracture) and giving mechanical parameters with Weibull distribution to rock elements. Mechanisms of fracture propagation and initiation have been fully investigated. This research gives a new approach for analyzing hydraulic fracture in shale, and its findings can offer a reference for hydraulic fracturing in shale reservoirs.

Methodology
In this research, numerical simulation software, named as rock failure process analysis (RFPA) [43,44], has been applied to simulate the hydraulic fracturing process. Currently, many scholars used RFPA to solve geotechnical engineering issues related to the failure mechanism of rock mass, showing its practicability [45,46]. In terms of finite element method and damage mechanics, this numerical simulation can express the failure process of heterogeneous materials, offering an effective description of microscopic damage mechanism and macroscopic failure. The characteristic of RFPA is illustrated:

Geofluids
(1) RPFA simulation is based on the evolution of failure of rock micro elements, resulting in rock macro failure. This failure process from micro to macro, consistent with real rock failure, cannot be easily achieved from normal simulation method, showing its unique and strong ability of simulating rock failure (2) RFPA has a good ability of expressing anisotropy in formation. For shale mechanical property, we normally assume that shale is composed of 2 parts, which are matrix and bedding plane. Correspondingly, different mechanical parameters are set for these 2 parts, expressing the anisotropy of shale. Furthermore, in the RFPA, a more systematical anisotropy could be illustrated, because it not only shows the anisotropy between matrix and bedding plane but also reveals the anisotropies of rock elements in matrix or bedding plane (each element has heterogeneity). All rock elements in RFPA have heterogeneity, which is expressed by assuming that the mechanical properties (elastic parameters and strength parameters of elements) conform to the Weibull distribution. Eventually, shale anisotropies from bedding plane, natural fracture, and rock elements are all considered in the simulation of RFPA, shown in Figure 2 (3) To sum up, hydraulic fracturing in shale formation is a process of gradually creating failure in anisotropy media. Since RFPA is able to remarkably exhibit failure process and anisotropy, it is reasonable to consider RFPA as an appropriate method for simulating hydraulic fracturing in shale    4 Geofluids

Constitutive Law.
It is generally recognized that the progressive degradation of rock properties results from fracture initiation, growth, and coalescence. Numerical simulation of the fracturing process must reflect the progressive degradation of rock subjected to loading. In this paper, the constitutive law of elastic damage mechanics has been used to express rock mechanics. According to the equivalent strain principle in damage mechanics, the constitutive law in damage can be derived from constitutive law in undamaged condition, shown as [47] where E and E o are the elastic modulus of the damaged and the undamaged rock, respectively, MPa; D represents the damage variable; σ is the effective stress, MPa; and ε is the effective strain, %.
For D = 0 and D = 1, rock is in undamaged condition and completely damage condition, respectively. When D is between 0 and 1, it indicates variable damage degrees. Initially, all elements are considered to be elastic, and their elastic properties are defined by Young's modulus and Poisson's ratio. The stress-strain curve of each element is considered to be linearly elastic until the given damage threshold is reached. The maximum tensile stress criterion and the Mohr-Coulomb criterion are chosen to find the damage threshold. Primarily, the tensile stress criterion is used to determine whether the element is damaged in    [48]. Correspondingly, the constitutive law and damage variable are illustrated in Figure 3. Figure 3(a) presents the damage constitutive relations with given specific residual strength, and the damage variable is written as Equation (3) [49]. When tensile strain of element reaches ε to , tensile damage occurs, increasing with strain growth. The complete damage (D = 1) is obtained under the limit tensile strain (ε tu ),   Table 1.  6 Geofluids reaching the state of tensile failure. The relation between damage variable and strain is expressed in Figure 3(b).
where σ to is the threshold of tensile strength, MPa; ε to is the threshold of tensile strain, %; σ tr is the residual tensile strength, MPa; ε tu is the limit tensile strain, %; λis the residual strength index; and σ 3 is the minimum principal stress of element, MPa.

Element Damage in Shear
Mode. The constitutive law in the previous section only considers the situation where the element is damaged in the tensile mode. In order to reflect the damage of element under compressive and shear stress, the Mohr-Coulomb criterion is chosen to determine damage threshold, expressed as [50] where ϕ is the internal friction angle of element, degree; σ 1 is the maximum principal stress of element, MPa; and σ c is the uniaxial compressive strength of element, MPa.
According to the constitutive law in Figure 4(a), damage variable in shear mode is shown as Equation (5) [51]. When compressive strain reaches maximum compressive strain, element begins to be damaged. The damage variable in this process is shown in Figure 4(b). It can be seen that the damage variable increases after the threshold of compressive strain and the residual strength still exists after shear failure.   7 Geofluids where σ co is the threshold of compressive strength, MPa; σ cr is the residual compressive strength, MPa; and ε co is the threshold of compressive strain, %.       9 Geofluids strength are required input parameters. For shale reservoir with parallel bedding plane, we divide shale sample into 2 parts, i.e., bedding plane and rock matrix. After transformation from micro scale to macro scale, all macro parameters are shown as follows.
For input parameters in Table 1, density, porosity, and permeability are directly acquired from rock physical test. Rock mechanical parameters in Table 1, such as uniaxial compressive strength of matrix and bedding plane, tensile strength of matrix and bedding plane, internal friction angle of matrix and bedding plane, elastic modulus, Poisson's ratio, and residual strength index, are obtained by trial method, depicted in Figure 5. According to the method in Figure 5, we firstly acquire shale uniaxial compressive strength and tensile strength using uniaxial test and Brazilian splitting test. Then, a series of input parameters (rock mechanical parameters in Table 1) are set for establishing RFPA simulation of uniaxial test and Brazilian splitting test. Consequently, uniaxial compressive strength and tensile strength from RFPA simulation can be acquired. Based on the comparison between experimental data and simulation results, input parameters are determined when experimental data has the best agreement with simulation results.

Model of Injection Wellbore in Formation.
The model is 2 m × 2 m square. In this paper, take horizontal well on the direction along minimum principal stress as an example ( Figure 6). The in situ stress and injection parameters are determined by oilfield data (logging, drilling, and geological data), shown in Table 2. In the boundary of this simulation, at the vertical and horizontal direction, stresses are constant as σ v and σ H , respectively. Also, displacements at boundary are zero. In homogeneity condition, simulation result is shown in Figure 7. Hydraulic fracture expresses typical dual-wing fractures, extending along the direction of maximum in situ stress, which is consistent with outcomes of other researches [52][53][54], proving the practicability of this model.

Model Verification. The verification has been conducted from 2 aspects:
(1) Verification can be confirmed from the comparison between experiment and simulation, shown in Figure 5 in Section 2.2. Uniaxial compressive strength and tensile strength are selected as parameters for verification. The reason why we choose these parameters is that shear failure mode and tensile failure mode both exist during hydraulic fracturing. In other words, hydraulic fracture is the result of rock shear failure and tensile failure. Uniaxial compressive strength and tensile strength represent those two types of failure modes. Thus, they can be regarded as useful parameters for validating whether this model is appropriate for the simulation of shale hydraulic fracture. Comparison results are shown in Figure 8. Also, relative errors between experimental data and simulation results are given in Figures 9  and 10. It is shown that simulation results have a good agreement with experimental data by using these input parameters. Relative errors are 6.47% (uniaxial test and its simulation) and 6.63% (Brazilian splitting test and its simulation), proving that input parameters of this model are appropriate for simulation  10 Geofluids (2) In addition, the model has been verified based on hydraulic fracturing curve of horizontal wellbore.
Since the pressure of hydraulic fracturing curve is associated with hydraulic fracturing initiation and propagation in the formation, hydraulic fracturing curve can be used to verify this model. According to FMI image from a vertical well (Figure 11(a)), we acquired a part of shale reservoir that has approximately 12.5 parallel bedding planes per meter. Correspondingly, at this part of reservoir, a horizontal well has been drilled and its simulation model has been built by RFPA, shown in Figure 11(a). Also, hydraulic fracturing curve of this horizontal wellbore has been given in Figure 11(b). Breakdown pressures of oilfield data (from hydraulic fracturing curve) and simulation are illustrated in Figure 11(c). It is found that breakdown pressures of oilfield data and simulation are 57.6 MPa and 54.3 MPa, respectively. The difference is merely 3.3 MPa, proving the practicability of this simulation

The Influence of Bedding Plane Occurrence on Hydraulic
Fracture. Bedding plane is a typical structure in shale reservoir, and its influence on hydraulic fracturing has been noticed in lots of researches. Thus, a model with bedding plane has been built for figuring out the propagation mechanism of hydraulic fractures in shale, shown in Figure 12. Results of this model are shown in Figure 13. It can be seen that with different included angles between bedding plane and maximum in situ stress, hydraulic fracture morphology changes. In homogeneity condition, propagation is totally controlled by in situ stress, extending along maximum in situ stress. In contrast, due to the impact of bedding plane, hydraulic fracture has the tendency of propagating along bedding plane. In high included angles (90 degrees and 75 degrees), propagation direction has slight inclination to bedding plane, but the influence of in situ stress is still on domination. Also, when influences of in situ stress and bedding plane have different directions, propagation of hydraulic  11 Geofluids fracture has two preferred directions. For stress effect, propagation direction should be along maximum in situ stress. For strength effect, propagation tends to be along weak strength part. As a result, hydraulic fracture propagates in the middle of bedding plane and in situ stress, showing bifurcation and relative larger width of hydraulic fracture. With small included angle, impacts from in situ stress and bedding plane are in similar direction. The propagation path is close to straight line with small width (Figures 13(e) and 13(f)). Meanwhile, when the first fracture point at the wall of borehole appears, pressure in wellbore can be acquired, referred as breakdown pressure. With bedding plane, there are 2 initiation types, i.e., initiation along rock matrix and initiation along bedding plane. It indicates that the breakdown pressure declines with decreasing included angle, shown in Figure 14. That is because initiation along bedding plane is likely to occur when the included angle is low.

The Influence of Bedding Plane Strength on Hydraulic
Fracture. The strength of bedding plane is different due to variable geological conditions. Therefore, in the condition of having bedding plane with variable strengths (Table 3), fracture propagation has been simulated, shown in Figure 15. With high bedding plane strength, the propagation path is in the direction between bedding plane and maximum in situ stress. In low bedding plane strength condition, the weak plane effect is stronger. Under the same stress condition, rock fracture tends to occur along weak strength area (i.e., bedding plane). Thus, the propagation path is directly along the bedding plane (Figure 15(b)).
When the strength of the bedding plane increases, the influence of the bedding plane declines and the impact of in situ stress gradually rises. Since influences of the bedding plane and in situ stress are in conflict, sporadic cracks appear around the main path of propagation, shown in Figures 15(c) and 15(d). Besides, it is noticed that the breakdown pressure declines with decreasing strength of bedding plane ( Figure 16). For strength 2 and strength 3, the difference of breakdown pressures is small, because both of them have a similar initiation mode, which is initiated from rock matrix, whereas for bedding plane having strength 1, the breakdown pressure has obvious decline due to the initiation along the bedding plane.
3.3. The Influence of Development Density of Bedding Plane on Hydraulic Fracture. In shale reservoir, the number of bedding plane is uncertain, which will cause a different impact on hydraulic fracture [56]. The fracture propagation      Figure 17.  13 Geofluids the bedding plane. Additionally, surrounding bedding planes are initiated, forming sporadic cracks around the main extension path and the bifurcation at the front area of hydraulic fracture (Figure 17(a)). With small development density (space distance is 80 mm), propagation of hydraulic fracture is simultaneously affected by in situ stress and bedding plane. The propagation direction has tendency of switching between in situ stress and bedding plane direction, shown in Figure 17(c). Furthermore, with increasing development density of bedding plane, breakdown pressure decreases, shown in Figure 18. But in the 20 mm and 40 mm space distance, breakdown pressures are almost the same. That is because both initiations happen along the bedding plane.
3.4. The Influence of Natural Fractures and Bedding Plane on Hydraulic Fracture. As mentioned above, except for bedding plane, shale reservoir contains many natural fractures, having an obvious impact on fracture initiation and propagation [57]. Natural fractures near the wellbore region have random distribution. Therefore, a model with bedding plane and random natural fractures has been built. Besides, under the combination of natural fracture and bedding plane, simulations in variable strengths of natural fracture, development densities, and angles of bedding plane have been conducted. In this simulation, liquid leakage has been ignored and all simulations are limited to the perspective of rock mechanics. Figure 19 demonstrates the hydraulic fracture propagation with variable strengths of natural fractures. Strengths of natural fractures are given in Table 4. On the path of hydraulic fracture propagation, natural fractures are gradually opened and extended. Particularly, during the propagation of hydraulic fracture, it is much easier for natural fractures to be open and extend due to their low strengths. Consequently, natural fractures connect with each other and bedding plane, forming dendritic bifurcation and having more complex fracture distribution. Furthermore, with bedding plane and natural fracture, breakdown point at the wall of borehole could be at rock matrix, bedding plane, or natural fracture. With decreasing strength of natural fracture, it is more likely to have initiation along natural fracture, causing low breakdown pressure, shown in Figure 20.
In addition, with decreasing development density of bedding planes, the impact of bedding plane on hydraulic fracture is clearly reduced, shown in Figure 21. It is obvious that the influence of bedding plane becomes stronger with decreasing space between bedding planes. With small space, more bedding planes exist in the formation. The weak strength effect along bedding plane direction is becoming stronger. In that case, hydraulic fracture has a greater tendency of propagating along bedding plane (the weak strength direction). On the other hand, when space distance increases, hydraulic fracture propagation is close to the horizontal path because in situ stress gradually takes control of hydraulic fracture propagation. For hydraulic fracture initiation, all initiation points are on natural fracture. However, breakdown pressure still has slight increment with decreasing development density of bedding planes, shown in Figure 22. That is because more weak planes could reduce the whole strength of formation, thus leading to the decline of breakdown pressure.
When bedding plane and natural fracture both exist, the variation of bedding plane angles is still an important factor, shown in Figure 23. At lower bedding plane angle, propagation direction is highly affected by the bedding plane. The propagation path is almost on the bedding plane. With increasing angle, the impact of the bedding plane wanes and propagation of hydraulic fracture mainly relied on in situ stress. It is noted that, even if bedding plane and natural fracture are both typical weak planes, the initiation point at the wall of borehole is highly controlled by natural fractures. With different bedding plane angles, the initiation point is always at natural fracture. Under the same initiation mode, breakdown pressure still changes with bedding plane angle. It is easier to break down the formation at low angle, shown in Figure 24. That is because bedding plane and in situ stress are both beneficial for hydraulic fracture to be opened along the vertical direction, giving rise to relatively small breakdown pressure in lower bedding plane angle.  14 Geofluids Furthermore, under certain structural plane characteristics, hydraulic fracture in variable in situ stress conditions has been simulated, shown in Figure 25. In this simulation, we set different ratios of σ H to σ v (Ri = σ H /σ v ) in situ stress conditions. It is shown that the influence of in situ stress is becoming stronger when stress difference grows. With high stress difference (Ri = 1:6), hydraulic fracture propagation is controlled by in situ stress, forming a single major fracture (like Figure 25). In contrast, with lower stress difference (Ri = 1:2), propagation direction is on the bedding plane. Besides, more branching fractures occur. Breakdown pressures with different in situ stress conditions are calculated, shown in Figure 26. Since all initiation points are at the natural fracture, breakdown pressure has little changes (<3 MPa) with variable stress differences.

Conclusion
In this research, according to the characteristic of shale reservoir, influences of the bedding plane and natural fractures on the hydraulic fracture around injection wellbore are analyzed. By using the rock failure process analysis (RFPA), hydraulic fracture initiation and propagation have been obtained. The following conclusions have been acquired.
(1) For shale reservoir with bedding plane, hydraulic fracture propagation is not only dependent on in situ stress but also related to the bedding plane. In high included angle, the propagation direction has slight inclination to the bedding plane, but in situ stress still has a larger impact on hydraulic fracturing. In large included angle, due to the contradiction between influences of in situ stress and bedding plane, bifurcation is clearly observed in propagation. With low degree of included angle, hydraulic fracture straightly extends along the bedding plane. Also, breakdown pressure declines with decreasing included angle because of initiation along the weak plane in low included angle (2) For low bedding plane strength condition, the propagation path is a straight line along the bedding plane. When bedding plane strength is larger, its influence on hydraulic fracture becomes small; the propagation path has the tendency of diverting to maximum horizontal in situ stress. Besides, sporadic    16 Geofluids cracks appear around the main path of propagation due to the conflicting impact between bedding plane and in situ stress. In addition, the influence of the bedding plane is stronger in high development density, and surrounding bedding planes are opened, showing obvious bifurcation and a more complex fracture network (3) With natural fracture and bedding plane, the whole propagation direction is still mainly affected by the bedding plane, meaning the propagation direction is still close to the bedding plane, especially with high development density and low angle. In the propagation, natural fractures are opened and further extended, connecting with each other and bedding planes, forming dendritic bifurcation, contributing to have a complex fracture network. Furthermore, the initiation point is normally at the natural fracture plane, leading to low breakdown pressure. Although the bedding plane can affect breakdown pressure, natural fracture is still the key factor of hydraulic fracture initiation

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
All authors declare that no conflict of interest exists.