The Influence of Bedding Planes and Permeability Coefficient on Fracture Propagation of Horizontal Wells in Stratification Shale Reservoirs

The complexity of hydraulic fractures (HF) significantly affects the success of reservoir reconstruction. The existence of a bedding plane (BP) in shale impacts the extension of a fracture. For shale reservoirs, in order to investigate the interaction mechanisms of HF and BPs under the action of coupled stress-flow, we simulate the processes of hydraulic fracturing under different conditions, such as the stress difference, permeability coefficients, BP angles, BP spacing, and BP mechanical properties using the rock failure process analysis code (RFPA-Flow). Simulation results showed that HF spread outward around the borehole, while the permeability coefficient is uniformly distributed at the model without a BP or stress difference. The HF of the formation without a BP presented a pinnate distribution pattern, and the main direction of the extension is affected by both the ground stress and the permeability coefficient. When there is no stress difference in the model, the fracture extends along the direction of the larger permeability coefficient. In this study, the in situ stress has a greater influence on the extension direction of the main fracture when using the model with stress differences of 6MPa. As the BP angle increases, the propagation of fractures gradually deviates from the BP direction. The initiation pressure and total breakdown pressure of the models at low permeability coefficients are higher than those under high permeability coefficients. In addition, the initiation pressure and total breakdown pressure of the models are also different. The larger the BP spacing, the higher the compressive strength of the BP, and a larger reduction ratio (the ratio of the strength parameters of the BP to the strength parameters of the matrix) leads to a smaller impact of the BP on fracture initiation and propagation. The elastic modulus has no effect on the failure mode of the model. When HF make contact with the BP, they tend to extend along the BP. Under the same in situ stress condition, the presence of a BP makes the morphology of HF more complex during the process of propagation, which makes it easier to achieve the purpose of stimulated reservoir volume (SRV) fracturing and increased production.


Introduction
In recent years, hydraulic fracturing has been extensively applied to increase the production rate and to realize long-term and stable yield in ultra-low-permeability shale reservoirs. Compared to the traditional fracturing volume techniques, the formation of fractures during shale formation is more complex owing to the existence of geological discontinuities such as natural fractures (NFs), bedding planes (BPs), and faults [1]. The ability to effectively con-trol the fracture formation as well as to make effective fractures remains a key problem in stimulated reservoir volumes (SRVs). A large number of studies have analyzed and researched fracture pressure under different conditions using different theoretical models [2][3][4][5][6][7]. The HF initiation is closely associated with the characteristics of shale rock mechanics and the heterogeneity of shale reservoirs. Chenevert and Mclamore [8][9][10] reported that the compressive strength of layered rock such as shale was a function of the confining pressure and the orientation of the anisotropic plane (BP or cleavage plane). Large quantities of NFs and BPs may change the initiation, propagation, and geometry of HF in a completely different way from those in isotropic and homogeneous media [11][12][13][14][15][16]. NFs and BPs may cause hydraulic fractures to extend along or through the structure surfaces [11,17,18]. Various forms of HF in conditions of different confining pressures, NFs, and crustal stress fields as well as the best boundary conditions for the production of many fracturing cracks have also been investigated [19][20][21]. The extension rule of HF after intersecting with NFs has been very clear [22][23][24][25], but when making contact with BPs of different bedding angles, the extension rule will be more complicated [26][27][28][29][30]. Under the effect of flow, the BPs will affect the fracture extension and result in fractures that develop in the direction of the original crack deviations. The complexity of HF geometries can be varied in layered formations because of the influence of BPs [31,32]. Thus, in order to analyze the borehole stability of horizontal wells and to optimize hydraulic fracturing design, it is necessary to further investigate the interaction mechanisms of HF and BPs considering the effect of coupled stress and flow. In addition, in order to obtain precise stress measurements and to enhance shale production, it is also important to explore the mechanism of HF and predict the geometry of the HF [33].
Numerical simulation experiments were carried out on the basis of the coupling flow-stress-damage model using a numerical code, which is referred to as rock failure process analysis (RFPA), which was developed by Tang C.A. [11,34,35]. It is known that the strength heterogeneities of rock mass are among the factors, and numerical experiments should be conducted to estimate how the formation of fractures is quantitatively influenced by the strength heterogeneities during the process of fracture network formation [36]. On the basis of a large number of studies conducted by our research group about the influences of the spatial location of NFs, the angle between NFs and the maximum principal stress, the length of NFs, the mechanical properties of NFs, the angle between BP and the maximum principal stress, BP spacing, and the mechanical properties of BP on initiation and propagation of the network cracks around the well [37][38][39], we further studied the effect of BP angle, BP spacing, BP compressive strength, and BP elastic modulus on the law of crack propagation considering permeability coefficient differences between the matrix and the BP of shale formation, and quantitatively analyzed initiation pressure and total break pressure under different BP angles. To implement the initiation of fractures and subsequent propagation, this study performs two-dimensional (2D) numerical simulations of the behavior of a cylindrical hole in the center of a shale reservoir subjected to different bedding inclination angles, BP spacing, BP mechanical properties, permeability coefficients, and an increasing injection pressure. From the results obtained from numerical experiments, we can acquire the initial pressure, total breakdown pressure, and fracture geometry at different conditions. In particular, it can be better implemented to the understanding of HF-BP interaction near wellbores by the numerical models when the main crack is closing to the BP and after contacting with the BPs.

Numerical Approaches of RFPA 2D -Flow
The RFPA 2D -Flow code can be used to simulate the progressive failure of heterogeneous and permeable rock material based on the finite-element theory and the statistical damage theory [35]. The four-node isoparametric element is applied as the basic element mesh. This coupled flow, stress, and damage (FSD) model in RFPA 2D -Flow has been validated in previous publications [18,35,40,41]. The progressive failure process of a quasibrittle material such as rock subjected to gradually increasing static loading can be simulated. The main governing formulations of the analysis are as follows: Geometric equation Constitutive equations : Seepage equation : Coupling equation Equations (1) to (4) are based on Biot's theory of consolidation [42], and equation (5) represents the effect of stress on the permeability, which is introduced to describe the dependency of permeability on stress and damage, where σ ji = stress, ρ = unit weight of rock,ε ij = strain, α = coefficient of pore water pressure, p = pore water pressure, λ = Lamé's coefficient, δ ij = Kronecher's delta function, G = modulus of shear deformation, Q = Biot's constant, K = permeability coefficient, K 0 = reference of permeability coefficient, a = a coupling parameter that reflects the influence of stress on the permeability coefficient, and ξð>1Þ = damage factor to account for the increased permeability of a material that is induced by a damage variable.
Continuum damage mechanics is used to describe the constitutive law of mesoscopic elements [43]. As shown in Figure 1(a), when the stress of the element satisfies the strength criterion (for example, the Mohr-Coulomb criterion), the element begins to be gradually damaged. With damage processes, the elastic modulus of the element may degrade step by step in elastic damage mechanics, which can be defined as follows: where D represents the damage variable, which is expressed as the ratio of microcracks, micropores, and defects in the material volume element. D = 0 indicates a state of reference, which indicates the absence of damage to the integrity of the material. D = 1 is equivalent to a complete loss of material. E and E 0 are the elasticity moduli of the damaged and undamaged materials, respectively. During elastic deformations, the rock permeability decreases as the cracks in the rock have a tendency to close under pressure, while it increases with the expansion and penetration of new fractures. A dramatic increase in rock permeability can be expected as a result of the generation of numerous microfractures. In other words, the permeability will increase significantly with damage to the rock. Upon reaching the peak load, the permeability may gradually drop again if the failed rock is further compacted, or the permeability may increase continuously as the failed rock is further extended [35].
When the tensile stress in an element reaches its tensile strength f t , the constitutive relationship illustrated in Figure 1(a) is adopted: The damage variable D can be described as follows [35]: where f tr is the residual tensile strength of the element and f t is the tensile failure strength of the element. ε t0 is the tensile strain at the elastic limit and is called the tensile threshold strain. When the value of the uniaxial tensile strain is ε t0 , the element begins to be damaged, but it does not immediately lose its bearing capacity. D decreases continuously as the degree of damage increases (0 < D < 1). ε tu is the ultimate tensile strain of an element at which the element completely loses its tensile load capacity. At this stage, the element will be in a full state of tensile fracture (damage), i.e., D = 1. In this case, the permeability can be described as follows: When the element is under uniaxial compression, the constitutive law is as shown in Figure 1(b). An element is considered to have failed in the shear mode when the compressive or shear stress has satisfied the Mohr-Coulomb failure criterion which is chosen as the second damage criterion [11,35,40,41,44]: where σ 1 ′ is the major effective principal stress, σ 3 ′ is the minor effective principal stress, ϕ is the friction angle, and f c is the compressive failure strength of the element. The damage factor under uniaxial compression is described as: where f cr is the residual tensile strength of the element, and ε c0 is the ultimate compression strain of the element. In this case, the permeability can be described as follows: The models consider the heterogeneity of material properties and the random distribution of defects, which is different from other mechanical software using the assumption of homogeneity. In fact, owing to the unequal distribution of defects in the rock medium, there is a large difference in their  3 Geofluids properties at the macroscopic and mesoscopic levels. Although the mechanical properties of the mesoscopic units are simple, the macroscopic nonlinear properties of material deformation can be reflected through the mesoscopic unit damage, and some complicated damage phenomena may still be described by their evolution. For heterogeneous rocks, the material's mechanical properties for different elements in RFPA are assumed to be randomly distributed throughout, conforming to the Weibull distribution [45]: where u is the mechanical property variable of the material element, such as elastic modulus, strength properties, or Poisson's ratio. u 0 is the corresponding average mechanical property, and m is the homogeneity index, i.e., a parameter defines the shape of the distribution function gðuÞ representing the degree of material heterogeneity; a larger value of m implies a more homogeneous material, and vice versa [35].
In RFPA-Flow, a given fixed loading is applied to the model incrementally in a quasistatic manner. Then, coupled flow-stress analysis is performed. The stress state of every element is then examined for failure before the next load step is implemented. The elastic modulus of each damaged element at every stress or strain level can be calculated using the above derivation of damage variable D as well as equation (6). Then, the analysis is restarted under the present boundary and loading condition in order to redistribute the stresses in the model without causing new damage. Finally, the increased external load (or displacement) is used as the input parameter for the analysis performed in the subsequent step.

Simulation Models of Hydraulic Fracturing
3.1. Characterizations of BP Realization. It is believed that the characterizations of BPs play a critical role on the response of stratification shale reservoirs to fluid injection [14,15]. The explicit representation of BPs with realistic characterizations is thus important in the numerical modeling. The geometric properties of BPs are often described by some statistical parameters, such as the BP angle distribution, spacing distribution, and length distribution. Combinations of these statistical characteristics and the mechanical parameters of the BP, such as the compressive strength and the elastic modulus, are essential for the fracture extension characterization of shale reservoirs.

Model Establishment.
The horizontal well in stratification shale reservoirs is drilled along the direction of minimum horizontal stress. Under the combined action of the maximum horizontal stress and the vertical stress, every vertical section is considered to be in the plain strain condition during the process of hydraulic fracturing. Figure 2 shows the geometry and the set-up of the simulation model. The model represents a two-dimensional (2D) vertical section of a stratification shale reservoir with inclined BPs.
The diameter of the injection hole is 0.15 m. The spacing of two adjacent parallel BPs is 0.05 m, which is determined from the average BP spacing of a real underground core. As shown in Figure 3, eight BP angle configurations were realized, which represented hydraulic fracturing of the horizontal well in the matrix without the BP, the hydraulic fracturing at different BP angles β of 0°, 15°, 30°, 45°, 60°, 75°, and 90°.
In the model, the injection goes through a horizontal wellbore in the center of the model. The increasing injection pressure is imposed on the wellbore at a constant rate. Simulation results of the pore water pressure at each step were calculated using the plain strain. The horizontal and vertical stress levels (σ H and σ V ) are, respectively, 45 MPa and 39 MPa (the horizontal stress σ H is derived from the field hydraulic fracturing data, and the vertical stress σ V is estimated using density logging data in Longmaxi shale, Sichuan Basin). The initial pore water pressure imposed in the well hole is 5 MPa, and a single-step increment is 0.5 MPa. The input material mechanical parameters for the numerical models shown in Table 1 are referred by laboratory shale experimental data according to equations (14) and (15), which were obtained from the Lower Silurian Longmaxi Formation in the Sichuan Basin of China [46,47]. Note that in Table 1, the porosity and density are considered to be uniformly distributed, ignoring the increase or decrease in the fracture aperture during rock mass deformation. The distribution histogram of the elastic modulus and the uniaxial compressive strength of the models are as shown in Figure 4: f Macro f micro = 0:2602 ln m + 0: where E micro and f micro represent the microscopic mean value of the elastic modulus and strength (input value of numerical calculation) when the Weibull distribution is assigned; E Macro and f Macro are the macroscopic elastic modulus and strength of the numerical sample, respectively; and m is the homogeneity index.  4 Geofluids According to the mechanical properties of shale reported by other researchers, the variation range of the elastic modulus anisotropy is 0.97-2.34, with an average of approximately 1.46 [47]. When combined with the analysis of the elastic modulus data of the Longmaxi shale derived from triaxial compression experiments [47], the maximum degree of anisotropy on the elastic modulus becomes approximately 1.25. It is known that anisotropic elasticity behavior has a relatively small impact on the stress distribution, especially when the degree of anisotropy is low (<1.25) [48,49]. The model used in this paper is a 2D model. Therefore, the mean value of elastic modulus E V in previous studies [47] is taken as the macroscopic elastic modulus of the model. According to the experimental data, the average macroscopic elastic modulus can be calculated, i.e., 42.51 GPa. Then, using equation (14), we acquire the microscopic elastic modulus with a value of 52.95 GPa. In order to facilitate the calculation, the elastic modulus of the matrix was determined to be 55 GPa in accordance with the law that the elastic modulus of the matrix is greater than the overall elastic modulus of shale.
According to the data of the direct shear experiment in the literature [47], the uniaxial compressive strength of intact rock and BP (weak plane) can be obtained by relying on the Mohr-Coulomb criterion after obtaining the cohesion and internal friction angle of the intact rock matrix and BP. Then, the formula is as follows: where σ C is the uniaxial compressive strength of intact rock or BP, C is the cohesion of the intact rock matrix or BP, and ϕ is the internal friction angle of the intact rock matrix or BP.

Geofluids
The uniaxial compressive strength of intact rock is about 61.89 MPa and is acquired using linear regression curves of the triaxial compressive strength with confining pressure from the literature [46], while a value of approximately 62. 35 MPa was obtained using the Mohr-Coulomb criterion based on the data obtained from the direct shear experiment in the literature [47]. Then, we can obtain the microscopic mean values of strength f micro (input value of numerical calculation) using equation (15), which are 200.18 MPa and 201.66 MPa, respectively. In order to facilitate the calculation, a value of 200 MPa was used during the calculation and analysis of the model. When the values of the uniaxial compressive strength of the intact rock and BP were calculated according to equation (16), the strength ratio between the intact rock and BP was also acquired, and ranged from about 1.89-3.98. This strength ratio is the source of the strength reduction ratio in Section 3.3. The ratio of the compression stress to the tensile stress is 10, which is recommended by the system. The value of Poisson's ratio is based on the triaxial experimental data in the literature [46], while the value of the internal friction angle is based on the direct shear experiment in the literature [47]. In order to study the influence of a single variable, Poisson's ratio and the internal friction angle of the BP and matrix are assumed to have the same value. Other parameters use the values recommended by software.
The permeability coefficient K is also known as the hydraulic conductivity coefficient. In an isotropic medium,   6 Geofluids it is defined as the unit discharge of the unit hydraulic gradient, indicating the difficulty of the fluid through the porous skeleton: where K is the permeability coefficient; k is the permeability of porous media, which is only associated with the properties of solid skeletons; η is the viscosity; μ is the kinematic viscosity; γ is the unit weight; and g is the acceleration due to gravity. The larger the permeability coefficient, the more permeable is the rock. The permeability coefficient value ranges for different types of rocks are shown in Table 2. After unit conversion of the permeability coefficient of shale in Table 2, the permeability coefficient applicable to the model can be obtained as a reference of the permeability coefficient in the model.

Numerical Procedure.
In most cases of triaxial mechanics experiments, the elastic modulus increases with the increase of compressive strength. Refer to Sun and Zhang's definition of strength ratio when he was studying and analyzing HF propagation in shale gas reservoirs [50,51]; the strength reduction ratio is defined as follows [52]: where R is the strength reduction ratio; f cm and f c0 are, respectively, the compressive strength of the BP and that of the matrix; and E m E 0 are, respectively, the elastic modulus of the BP and that of the matrix. The reduction ratio R, which is adopted in simulation experiments without special instructions, is 0.25. The focus of this work is to study the interaction between the HF and the BP, and to explore the effect of the stress difference, permeability coefficient, BP angle, BP spacing, and the mechanical parameters of the BP on the initiation and extension of HF. A series of comparative analyses were performed to investigate the interactions between the HF and the BP. In addition to a qualitative evaluation of simulation results, the models' responses were compared and evaluated in terms of specific indices during injection. These indices include the (1) initiation pressure, which is defined as the pressure when the rock is damaged, and (2) total breakdown pressure, which is defined as the pressure when the rock mass is almost completely covered by HF and the opening BP, i.e., the rock model will be completely destroyed in the next calculation. We believe that the formation has completed SRV fracturing under this pressure.

Results and Discussion
4.1. Effects of Permeability Coefficient on Fracture Initiation and Propagation. Numerous laboratory and field studies have been conducted to show that not only is HF branched and nonplanar fracture growth possible but they are also fairly common [53]. Four groups of matrix models without BPs were compared, and the permeability coefficient and stress values are as shown in Figure 5. As the injection pressure increases, the models of stress difference at 0 MPa  Figure 5(b). However, the HF values of the models without a stress difference spread along the direction of the larger permeability coefficient as shown in Figure 5(d) (K v = 0:00008 m/d, K h = 0:00005 m/d). This indicates that the permeability coefficient has a small impact on the fracture initiation, and it mainly affects the fracture extension morphology. In Figures 5(e) and 5(f), when the permeability coefficient is uniformly distributed, the extension of HF of the model with a stress difference of 6 MPa is in the direction of the maximum principal stress. In Figures 5(g) and 5(h), when the permeability coefficient varies in the vertical and horizontal directions, the extension of HF of the model with a stress difference of 6 MPa is still mainly in the direction of the maximum principal stress. The simulation results show that in the model without a BP, when the stress difference is 6 MPa and the permeability coefficient K v = 0:00008 m/d, K h = 0:00005 m/d, the stress difference is the main influencing factor on the fracture extension. Therefore, the influence of the stress difference on the extension direction of HF should also be considered when studying the influence of the permeability coefficient. This is because the permeability coefficient of a BP far exceeds that of the matrix; in addition, the strength of the mechanical parameters of the BP is relatively low, so the fracture system of the stratification shale will become more complex after making contact with BPs, which warrants further study.

Effects of BP on Fracture Initiation and Propagation. The
BPs in layered sedimentary rocks influence the HF growth because of changes in rock properties and in situ stresses associated with the layers. Offsets in the fracture pathways have been documented in man-made HF that have been mined and mapped [53]. During the hydraulic fracturing process, a shear zone and open zone will appear in the BP,   Figure 5: Effect of the permeability coefficient on HF morphologies at different stress difference models without BP. The color shadow indicates the relative magnitude of the pore water pressure field (K h is the horizontal permeability coefficient of the matrix; K v is the vertical permeability coefficient of the matrix). 8 Geofluids changing the flow path of the fracturing fluid and affecting the extension path of the fracture subjected to the fracture tip stress field and the original in situ stress from previous studies [1,54]. Using numerical simulations based on a 2D boundary element model, Zhang et al. studied the potential mechanism of fracture deflection and propagation in hydraulic fracturing as well as the subsequent fluid invasion at the friction bedding interface. The growth of fluid-driven fractures along BPs will alter fracture growth and fluid in every direction, and will affect the overall fracture behavior. In that study, the branching of the hydraulic fracture, which is initially perpendicular to the bedding contact, is controlled by the frictional coefficient of the interface, elastic properties of the layers, remote stress condition, and injected fluid viscosity. The natural system that they modeled is a bedded sedimentary rock containing a single fluid-driven fracture confined to one layer. HF propagating initially perpendicular and towards the BPs can be deflected in the BP to create two daughter branches in the interface postintersection. The difference in the variation of the injection pressure was analyzed when the parent fracture is located in a rigid layer or soft layer [53]. Chuprakov et al. used numerical modeling to quantify the physical mechanisms of the mechanical activation of a natural fault due to contact with an HF [55]. An analysis of the total stress state induced along the NF is fulfilled numerically for different stages of hydraulic fracturing (HF tip approaching, coalescence, and fluid infiltration along the NF) [55]. To study the influences of the interaction of HF and preexisting NFs on the complex fracture nets, Gu et al. developed a criterion to determine whether a fracture crossed a frictional interface at nonorthogonal angles [56]. Zhang et al. believed that there was a potential state between the two extremes of interaction of HF and a BP in sedimentary rocks, which is when HF penetrates through the BP and HF may be arrested or blunted at the bedding contact, i.e., the fracture and the fluid flow were deflected into the BP and were divided into two branches. If there are flaws on the interface, potential reinitiation of a new fracture from one flaw will leave a step-over at the BP [53]. In this paper, we consider the injected fluid to be water, and the fracture can be allowed to penetrate to the adjacent layer or to induce a new fracture at some location along the interface, which does not occur in the study by Zhang et al. Upon completion of the open hole, the injection pressure is in the well-bore, and cracks can initiate and extend in the rock around the wellbore, so there will be many HF. According to the model results in Figures 5(g) and 5(h), when the horizontal bedding coefficient is not equal to the vertical bedding coefficient, and the stress difference is 6 MPa, the influence of in situ stress on fracture extension is more significant. If a single BP is added on the basis of this model in Figure 5(g) for hydraulic fracturing, there is the need for further study to determine whether the expansion morphology of HF will be affected. The simulation results show that the addition of a single BP does not affect the initiation of HF, but it has a certain effect on the extension of HF. In the extension stage of fractures, if the fracture makes contact with a single BP, it will expand along the BP; if not, it will mainly extend along the direction of the maximum principal stress ( Figure 6). The numerical simulation results are in agreement with the results obtained by the above scholars.  [53]. Simulation results showed that the changes in BP angles were closely related to the damaging and cracking processes in the stratification shale. Finally, fracture geometrical morphologies during the process of HF at different bedding angle models are not identical. The HF intersected and then deflected into the BP. According to the classical theory of rock mechanics, under the hypothesis that the strata rock is continuous, homogeneous, and isotropic, the initiation and extension of the induced crack is always along the maximum principal stress orientation [3,4]. However, based on the numerical simulation results (Figures 7 and 8), the initiation of the heterogeneous and anisotropic medium is not along the absolute maximum principal stress orientation, but it is close to the maximum principal stress orientation. Because mechanical properties of elements are heterogeneous, units having borehole walls with weaker mechanical properties may be preferred to achieve the extent of damage, leading to deviations in the induced crack from the direction of maximum principal stress. When the wellbore fluid column pressure  Figure 6: Interaction between single BP and HF. 9 Geofluids gradually increases to a certain value, because of the strong heterogeneity, the crack is not straight but is bent outward, and it forms the branch away from the borehole area.
As the permeability coefficients of the matrix and stratification are low, from the results of the numerical simulation (Figure 7), the impact of a weak BP on crack initiation and     10 Geofluids extension is greater than that of the maximum principal stress at this time: (i) When β = 0°, the HF cracks and extends along the BP through the largest tensile stress area (ii) When the value of β ranges from 15°to 75°, the HF cracks and extends along the symmetrical BPs through the wellbore where fractures are initiated. Under the action of the tensile stress, two major symmetric fractures are formed, and the fractures are mainly straight cracks (iii) When β = 90°, fracture initiation and propagation mainly extend along the BP through the hole, and largely deflects from the direction of the maximum principal stress on the right side. Deflection phenomenon may be the result of mutual influence resulting from the angle of the BP, the heterogeneity distribution of a weak unit, and the maximum principal stress Based on the simulation results of the high permeability coefficient, the following conclusions are gained (Figure 8). Compared to models with the low permeability coefficient, the models with high permeability coefficients for the BP and rock matrix have different fracture morphologies (initial seam and extension directions) at a high BP angle owing to the stress field. When the BP angle is small (0°~45°), fractures initiate and propagate along the BP direction, and the main control factor is the BP. When the BP angle is high (60°~90°), the initial seams and extension direction are influenced by the BP and the principal stress.
When β = 60°, the HF is initiated at the intersection of the wellbore and BP, as with the models with bedding angle values ranging from 15°to 45°. The lower branch of the HF extends along the BP because the place where the fracture is initiated is the tangent point between the weak BP and the left side of the wellbore. Then, the adjacent BP is opened with many microfractures after they are infiltrated by fluid. On the right side of the wellbore, before the upper branch of the HF is deflected to the adjacent BP direction, it first extends along the direction of maximum principal stress under the effect of the stress difference ( Figure 9). In Figures 9(a)-9(d), the fluid column pressures in the wellbore are 56 MPa, 57.5 MPa, 58.5 MPa, and 60 MPa, respectively. With the increase of the pore water pressure, the large tensile stress zone was formed at the tip of the main fracture along the BP [1,[11][12][13]. Partial BPs at the tensile stress zone are ripped, thus producing many small fractures that are parallel to the main fractures. Under the influence of the stress difference, compared with the main fractures, the small fractures tend to the direction of the maximum principal stress [57]. In this model, we can observe the phenomenon that the fracture diffuses to the adjacent BP.
When β = 75°, the lower branch HF cracks along the direction of the maximum principal stress through the wellbore, and the BP direction was opened simultaneously; however, the upper branch HF cracks in the direction of the maximum principal stress. Under the combined action of a high permeability coefficient and the maximum principal stress, the HF cracks and first extends along the direction of the maximum principal stress, and it then deflects to the intersecting BP, and the adjacent BPs are ripped. This is because the crack point is at the border of the BP and wellbore rock, and the BP is a weak surface, and outspreading along the BP needs less energy. The HF of the model with a BP angle of 75°occur at the phenomenon where HF go through the BP, and HF are longer than that of the model with a BP angle of 60°. With the increase in the pore water pressure, many small fractures appear near the maximum principal stress compared with those at a bedding angle of 60°. The secondary cracks appear on the clockwise side of the main crack. These phenomena prove that for different BP angles under the action of the same stress field, HF and an opened BP in hydraulic fracturing are different, and for the model at a BP angle of 75°, it is easier to produce a complex fracture network compared with that at a BP angle of 60°.
When β = 90°, the fracture morphology is totally different from that under the low permeability coefficient. In this condition, the influence of the BP on the fracture is very small and is mainly affected by the maximum principal stress. Therefore, the value of the permeability coefficient in the study of a fracture network of bedding shale formation cannot be ignored.
In addition to the qualitative evaluation of the models' responses, some indices are chosen to describe the injection process. In Table 3, K and K ′ are the permeability coefficients of the matrix and BP, respectively. P f 1 and P f 1 ′ are the initiation pressures acquired by the pore water pressure, P f 2 and P f 2 ′ are the initiation pressures acquired by acoustic Step = 103, 5 MPa-0.5 MPa (a) Step = 106, 5 MPa-0.5 MPa Step = 108, 5 MPa-0.5 MPa

(c)
Step = 111, 5 MPa-0.5 MPa 11 Geofluids emission (AE), and P f 3 and P f 3 ′ represent the total break pressures acquired by the pore water pressure. The details regarding the method to acquire the initiation pressure and total breakdown pressure can be obtained from the studies by Yang et al. [18]. The initiation pressure and fracture pressure show a linear increasing trend as the BP angle increases, and its rate of increase is basically the same except for the model with β = 15° (Figures 10 and 11). This illustrates that the model whose BP is parallel to the direction of the maximum principal stress will fail most easily. The smaller bedding angle is advantageous to the fracture extension in the direction of the maximum principal stress during the hydraulic fracturing process. For the same BP angle under the same stress state, the models with the larger   Acquired by pore water pressure (K 2 ) Acquired by acoustic emission (K 1 ) Acquired by pore water pressure (K 1 ) Figure 11: Initiation pressures acquired by pore water pressure nephogram or acoustic emission nephogram under different BP models. 12 Geofluids permeability coefficient have lower initiation pressures and total breakdown pressures. The changes in the permeability coefficients of the BP and rock matrix affect the propagation pattern of the HF with respect to the rock mass, initiation pressure, and total breakdown pressure [32,51]. The increase in the permeability coefficient weakens the influence of the BP on HF, and it enhances the impact of the maximum principal stress. AE monitoring has been used to reveal the spatial distribution and hypocenter mechanisms of AE events induced by rock failure [1,58]. Ning et al. reported that shear and tensile events were induced in hydraulically connected regions, and shear events also occurred around BPs that were not hydraulically connected based on the analysis of the hypocenter mechanisms [1]. In this study, the initiation pressure obtained using the microseismic monitoring of an AE is greater than that of the pore water pressure nephogram (Table 3), and this is because when a minor injury occurs, the AE is not detected. AE under different BP angle models at the final step was shown in Figure 12. AE takes place mainly near the BPs, and the phenomenon is consistent with the pore water pressure nephogram.

Effects of BP Spacing.
In the following studies, without special indication, the model adopts the low permeability coefficient K 1 . In order to investigate the impact of BP spacing on fracture initiation and propagation, we varied the spacing between two adjacent parallel BPs from 0.05 m to 0.10 m, and other parameters were the same as those in the previous model with the low permeability coefficient, K 1 . For a BP spacing of 0.05 m, the extension of the induced joint is limited in BPs. Compared to the BP spacing of 0.05 m (Figure 7), when BP spacing is 0.10 m, we obtain the following results (Figures 7 and 13). When β is between 0°and 15°, the initiation and extension of fractures easily occur on two adjacent different BPs rather than the single one throughout the wellbore. At β = 30°, lower branch HF extend along the BP. The upper branch cracking angle is along the direction of the maximum principal stress and finally extends along the intersecting BP. At β = 45°, straight cracks extend along the BP throughout the borehole. When the β values are 60°and 75°, the fracture makes more BPs open and connected. At β = 90°, the HF deflects and diffuses to form more branches far from the wellbore, but its extension direction is still primarily that of the maximum principal stress. With the increase in the BP spacing, the influence of stress increases, while the influence of the BP decreases. Far from the wellbore area, as the stress concentration decreases, adding high development BPs result into other BPs opened around the induced joint and forming more branch fractures [40,41].

Geofluids
Under a pore water pressure of 55 MPa (Figure 14), the models with different BP compressive strengths yielded different fracture lengths. As the compressive strength of the BP increases, the length of the fracture gradually decreases.
When the BP strength is less than 50 MPa, fracture extends along the two different adjacent BPs. The induced joint is relatively straight. When the BP strength exceeds 100 MPa, the fracture deflects along the direction of the maximum principal stress and then extends along the intersecting BP at the fracture tip. As the bedding strength increases, the straight fractures along the BP evolved into curving fractures in an irregular manner ( Figure 15). Therefore, the weaker the compressive strength of the BP, the more easily will the induced joints extend along the BP. With the augmentation of the compressive strength of the BP, the induced joint extension gradually deviates from the BP direction to the direction of the maximum principal stress, and some of the BPs may be opened during the development process. As shown in Figure 16, the change in the stiffness of the BP has almost no effect on the final failure mode of the model. Fractures initiated and extended along two symmetrical BPs tend to the direction of the maximum principal stress. Owing to the low strength of the BP, the BP is always the main control factor in the evolution of the fracture.

Effects of Strength Reduction Ratio.
For strength reduction ratios of 0.25, 0.5, and 0.75, the numerical models for simulation analysis at β = 0°, 15°, 30°, 45°, 60°, 75°, and 90°w ere established in order to determine the comprehensive effects of the compressive strength and elastic modulus on fracture initiation and propagation.
Numerical simulation results indicate that there are various forms of fractures at different BP angles and different reduction ratios (Figures 7 and 17). For a reduction ratio of 0.25, fractures extend mostly along the BP. The specific analysis can be seen in Section 4.2.1. Under a high reduction ratio, bifurcate extension occurs easily at the area far from the wellbore. Secondary cracks are formed in the main fracture edge, and these secondary cracks continue to branch, forming more multistage secondary cracks. A complex fracture network is induced by the branching of microcracks around the tip of the main fracture and pores, depending on the uniformity of the surrounding stress field [36].
Eventually, the main fracture and secondary cracks are intertwined, the fracture network system is formed, and a wide range of effective communication within the reservoir has been realized [32,40,41,44].
When β changes from 0°to 45°, under a reduction ratio of 0.5, fractures are initiated along two symmetrical BPs tending to the maximum principal stress. For a reduction ratio of 0.75, fractures are initiated between the direction of the horizontal maximum principal stress and the BP with a smaller deflection angle. Figure 13: HF morphologies under different BP angles; BP spacing is 0.10 m. 14 Geofluids When β = 60°/75°, at the lower reduction ratio of 0.25, the fracture initiation direction is not entirely consistent with the BP, and fractures propagate to neighboring BPs (Figures 7(f) and 7(g)). At the higher reduction ratios of 0.5 and 0.75, the fracture is initiated and propagates along the direction at a certain angle with the BP. Their failure mode develops into the model without the BP (Figures 17(i)-17(l)).
When β = 90°, for reduction ratios of 0.5 and 0.75, the fracture is initiated and propagates to form a large deflection angle along the direction of the horizontal maximum    Based on the number of loading steps, as the reduction ratio increases, the total breakdown pressure increased (Figures 17(a)-17(j)). When β is larger than 75°, the total breakdown pressure remains nearly unchanged at the reduction ratios of 0.5 and 0.75.
As the reduction ratio increases, there is a smaller difference between the mechanical properties of the BP and the mechanical properties of the matrix, and the maximum principal stress becomes the main factor that affects fracture propagation. The initial cracking position of the fracture is always approximately in the direction of the maximum horizontal principal stress, and the direction of the crack extension varies slightly with the BP.

Conclusions
In this study, the interaction of HF and BP in the hydraulic fracturing process has been explored by performing numerical simulation experiments. The results can provide the basis for the analysis of horizontal well borehole stability and the optimization of hydraulic fracturing designs.
(i) When there is no stress difference in the models without the BP and the permeability coefficient is uniformly distributed, HF spreads outward around the hole; however, HF extends along the direction of the larger permeability coefficient in the rock in the case where the permeability coefficient is not uniformly distributed. Whether the permeability coefficient is uniformly distributed or not, HF spreads along the direction of the maximum principal stress when the models without the BP have a stress difference of 6 MPa (ii) HF geometry formations for different BP angle models are not identical. When the permeability coefficients of the matrix and BP are low, the impact of the BP on fracture initiation and extension is greater than that of the maximum principal stress. When the permeability coefficients are large, for low BP angles (0°~45°), fracture initiation and propagation are along the BP direction and the dominant factor is the BP. For high BP angles (60°~90°), the influencing factors are the BP and the principal stress. The model with a higher BP angle may no longer propagate along the BP (iii) The initiation pressure and total breakdown pressure show a linear increasing trend as β increases, and its rate of increase is basically the same. The model whose BP is parallel to the direction of the maximum principal stress most easily fails. For the same BP angle, the initiation pressure and total breakdown pressure at low values of the permeability coefficient are larger than those with high values (iv) With the increase in either the BP spacing, the compressive strength of the BP, or the reduction ratio, the influences of the BP on fracture initiation and propagation decrease, and the cracks more easily bifurcate or deflect towards the direction of the Step = 108, 5 MPa-0.5 MPa (a) Step = 107, 5 MPa-0.5 MPa (b) Step = 106, 5 MPa-0.5 MPa (c) Step = 106, 5 MPa-0.5 MPa

(d)
Step = 105, 5 MPa-0.5 MPa (e) Figure 16: Effects of BP elastic modulus on final fracture morphology. 16 Geofluids maximal horizontal principal stress. The change in the stiffness of the BP has almost no effect on the final failure mode of the model

Data Availability
The data used to support the findings of this study are included within the article.

Conflicts of Interest
The authors declare no conflicts of interest.