Effect of rock brittleness on propagation of hydraulic fractures in shale reservoirs with bedding-planes

© 2020 The Authors. Energy Science & Engineering published by the Society of Chemical Industry and John Wiley & Sons Ltd. Hydraulic fracturing forms complex hydraulic fracture networks (HFNs) in shale reservoirs and significantly improves the permeability of shale reservoirs. Although rock brittleness is a major factor in determining whether a shale reservoir can be fractured, the relationship between HFNs and rock brittleness remains unclear. To investigate this relationship in a shale reservoir with bedding planes, this paper presents a series of hydraulic fracturing simulations based on a hydromechanically coupled discrete element model. In addition, we analyzed the sensitivity of the difference in rock brittleness to bedding-plane angle and density. The parameters used in the model were verified by comparing the simulated results with experimental results and a theoretical equation. The results showed that breakdown pressure and injection pressure increased with increasing rock brittleness. The tensile hydraulic fracture of a shale reservoir (THFSR) was always the most abundant type of hydraulic fracture (HF)—almost 2.5 times the sum of the other three types of HFs. The distribution of areas with higher fluid pressure deviated from the direction of the maximum principal stress when the angle between the bedding plane and maximum principal stress directions was large. Upon increasing this bedding-plane angle, the breakdown pressure and rock brittleness index first decreased and then increased. However, regardless of bedding angle, the relative proportions of the various types of HFs remained essentially constant, and the seepage area expanded in the direction of the maximum principal stress. Increased bedding-plane density resulted in a gradual increase in the total number of HFs, with significantly fewer of the THFSR type, and the large seepage areas connected with each other. This study thus provides useful information for preparing strategies for hydraulic fracturing. Disciplines Engineering | Science and Technology Studies Publication Details Chong, Z., Yao, Q., Li, X., Zhu, L. & Tang, C. (2020). Effect of rock brittleness on propagation of hydraulic fractures in shale reservoirs with bedding-planes. Energy Science and Engineering, This journal article is available at Research Online: https://ro.uow.edu.au/eispapers1/3878

tigate this relationship in a shale reservoir with bedding planes, this paper presents a series of hydraulic fracturing simulations based on a hydromechanically coupled discrete element model. In addition, we analyzed the sensitivity of the difference in rock brittleness to bedding-plane angle and density. The parameters used in the model were verified by comparing the simulated results with experimental results and a theoretical equation. The results showed that breakdown pressure and injection pressure increased with increasing rock brittleness. The tensile hydraulic fracture of a shale reservoir (THFSR) was always the most abundant type of hydraulic fracture (HF)-almost 2.5 times the sum of the other three types of HFs. The distribution of areas with higher fluid pressure deviated from the direction of the maximum principal stress when the angle between the bedding plane and maximum principal stress directions was large. Upon increasing this bedding-plane angle, the breakdown pressure and rock brittleness index first decreased and then increased. However, regardless of bedding angle, the relative proportions of the various types of HFs remained essentially constant, and the seepage area expanded in the direction of the maximum principal stress. Increased bedding-plane density resulted in a gradual increase in the total number of HFs, with significantly fewer of the THFSR type, and the large seepage areas connected with each other. This study thus provides useful information for preparing strategies for hydraulic fracturing.

K E Y W O R D S
bedding plane, hydraulic fracture network (HFN), hydromechanical coupling, rock brittleness, shale reservoirs

| INTRODUCTION
Major advances in hydraulic fracturing technology have enabled a rapid increase in worldwide shale-gas production, 1,2 exceeding 20 billion m 3 in China's fuling gas field alone. Hydraulic fracturing creates a hydraulic fracture network (HFN) in a reservoir, thereby enhancing the reservoir's permeability and conductivity. 3,4 A pump injects high-pressure fluids into the reservoir to generate hydraulic fractures (HFs) that interact with and activate the bedding planes (BPs), generating massive HFNs in the reservoir. To date, significant time and effort have been devoted to theoretical and numerical modeling and physical experimentation aimed at elucidating the reciprocal relationship between HFs and BPs and to investigate HF characteristics in order to predict the resulting HFN. 5,6 Previous studies have shown that the following several factors contribute to creating complex HFNs in shale reservoirs. (a) The viscosity of the fracturing fluid: low-viscosity fracturing fluids, which have relatively less surface tension, can easily penetrate into micropores, and therefore these fluids spread over wider areas, resulting in more extensive HFNs. 7 (b) The abundancy of BPs in the reservoir 8 : Fracturing fluids permeate more extensively in reservoirs with abundant BPs, where the pore pressure distribution is such that HFs can initiate and propagate in various directions. (c) The difference between maximum and minimum principal stresses acting on the reservoirs 9 : As the ratio of maximum principal stress to minimum principal stress increases, HFs more easily traverse BPs, and this extended propagation results in long, narrow HFNs. However, when this ratio is smaller, HFs propagate in various directions, resulting in more complex HFNs. (d) Bedding planes: Some parameters such as BP angle and density can affect propagate track of HFN. (e) Inherent reservoir characteristics 10 : Reservoir characteristics such as the rock brittleness index B also affect HFN formation and can be crucial to fracturing efficiency.
To date, no consensus has been reached as to how rock brittleness affects HFN formation. In rock engineering, rock brittleness is among the most important mechanical properties of rock and an important factor in the determination of when rock failure is likely to occur. However, researchers have given various definitions of brittleness. Table 1 summarizes nine such methods proposed over the last two decades. These methods can be divided into three types, based on whether they define brittleness in relation to strength, 11,12 strain, 13,14 or modulus. 15 Mineral composition is the primary factor affecting the rock brittleness of shale reservoirs. 16,17 Numerous works have established models to predict B from the mineral composition, as summarized in Table 2.
In addition, field tests, test data, and mineralogical reports often result in different assessments of rock brittleness because these methods are based on different assumptions. For example, Obert and Duvall 18 claim that rock brittleness is characterized by the point where a rock sample fails as it reaches or slightly exceeds its yield strength, whereas Evans et al 19 define the brittle index as 1% of the deformation when the rock fails. In the Evans et al definition, if a rock fails at less than 1% deformation, then it is classified as brittle rock, and failure at 1%-5% deformation is classified as brittle-ductile rock. In another approach, Tarasov and Potvin 15 characterize rock brittleness based on the rock's ability to resist failure when subjected to the combined effect of its anisotropy and external loads.
The most common method of defining B was proposed by Rickman, 20 who combines both Young's modulus E and Poisson's ratio . These two indices can be obtained directly from field data and do not require physical experiments on rock core or analyses of the composition of the rock. Langenbruch and Shapiro 21 studied how this definition of rock brittleness B relates to rock failure caused by HF and applied an analytical stress wave solution to real shale reservoirs. Shimizu 22 Based on modulus 7 E is Young's modulus and M is postpeak modulus. 15 the connection between rock brittleness and the formation of HFNs remains unclear. In such reservoirs, rock brittleness depends not only on Young's modulus and Poisson's ratio but also on both the bedding angle and bedding density.
In this work, we used the discrete element method (DEM) to simulate Rickman's brittleness, and then compared the numerical data with data from physical experiments and a theoretical equation to identify the microparameters of shale reservoirs with different layer orientations. The resulting microparameters were then used to establish an HF model that considers BPs. This model describes how rock brittleness affects HFNs in terms of breakdown pressure, fluid pressure distribution, number and formation of HFs, and seepage area. Using this model, we also analyzed the sensitivity of HFN formation to bedding angle and bedding density.

NUMERICAL SIMULATION OF FLUID FLOW
The numerical model was implemented in the two-dimensional (2D) particle flow code (PFC 2D ). PFC 2D is a classical 2D DEM based on circular particles that contact each other through a simple contact logic instead of through complex constitutive models. In PFC 2D , a rock mass is modeled as an assembly of rigid circular particles. The rock matrix is represented by using the particle bond method (PBM), and the BP is represented by using the smooth-joint method (SJM).
In PBM, the contact between adjacent particles is modeled as a series of springs 23 (Figure 1). The particle motions must satisfy only Newton's second law of motion and Hooke's law, and complex constitutive relations do not need to be defined for the whole model because after deformation, the initially defined constitutive model may no longer be applicable. Moreover, adjusting the timing and parameters of a constitutive model of hydraulic fracturing throughout the calculation process is quite difficult. Therefore, among available numerical software programs, PFC 2D is the most suitable for simulating fracture initiation, propagation, and connection behaviors in a reservoir subjected to hydraulic fracturing. PFC 2D also represents BPs by the SJM, which is a common approach. Once a BP is identified, the original contact between two particles is replaced by a microscale slip surface, as shown in Figure 2, such that two particles can overlap or pass through each other. 24 The SJM algorithm is explained in the PFC 2D manual, and therefore, Figure 2 shows a simple overview of the SJM.
Based on the rock deformation and failure process applied in PFC 2D , hydraulic fracturing can be simulated by a viscous fluid flow network implemented in an independent subroutine. As shown in Figure 3, the polygon formed by the lines connecting the centers of adjacent particles is called the fluid domain, and the contacts between particles form the borders of the fluid domain, as shown by the blue polygon in Figure 3B.
Each fluid domain is a storage unit, and these units are connected to each other through flow channels, indicated by the yellow line in Figure 3B. According to the Poiseuille equation, fluid flow occurs when there is a pressure difference between different fluid domains. The fluid flows laminarly between parallel plates. Therefore, the volumetric flow rate, q, per unit width is 25,26 : (1) q = a 3 12 where µ is the dynamic viscosity of the fluid, in Pa·s; P 1 and P 2 are the fluid pressures in two neighboring domains, in Pa; and L and a are the length and aperture of the flow channel, respectively, in m. The accumulated fluid pressure within a fluid domain leads to a pressure difference between different fluid domains. The pressure within a given domain acts on the surfaces of the particles surrounding the domain. The change in fluid pressure, dP, is given by: where K f is the bulk modulus of the fluid, in Pa; dt is the time increment, in s; and V r is the pore volume, in m 2 (two dimension).
Although dV r is often ignored in complex hydromechanically coupled models, we include it here to obtain more realistic results from the simulations.
The value of dV r is related to the deformation of the fluid domain caused by fluid pressure, as shown in Figure 3C. Within the fluid domain, the component f i of fluid pressure is perpendicular to the border of the domain and acts on other particles. The component f i can be calculated from the fluid pressure P f , the chord length l i , and the unit thickness of the model. If the bond is broken, l i is defined as the distance between the centers of the two nearby boundaries of a domain.
According to experiments, the reservoir permeability is directly related to the in situ stress. Therefore, in PFC 2D , the aperture of the flow channel is considered to depend on the normal stress n exerted on the parallel plate, which can be described as follows: where a inf and a zeo are the apertures at infinite and zero n , respectively, in m, and the coefficient is the rate at which the aperture decays with increasing n , in Pa (usually, = −0.15 27 ). When the normal stress n = 0, the particles in the reservoir model just contact each other, and when n approaches infinity, the aperture approaches a inf . Equation 3 allows for the permeation of fluid into the reservoir without pre-existing fractures, which accurately simulates the leak-off phenomenon. However, in the equation, the aperture of the flow channel is calculated from the microscope flow rate instead of directly from the permeability of the reservoir. To remedy this problem, the apertures a zeo and a inf are calculated by using where k and V are the permeability and total reservoir volume in d and m 2 , respectively.
Changes in the fluid flow that occur after bond breaking induces fractures must also be considered within this model. If fractures appear in connected domains, the fluid flow within the flow channel changes instantaneously. Therefore, after the bond breaks, the fluid pressure P ′ f is assigned the average of the two pre-bond-break fluid pressures in the two domains: where P fA and P fB are the fluid pressures of the two reservoir domains before the bond breaks, in Pa.
However, if the volumes of the two domains differ significantly, the fluid pressure after the appearance of fractures induced by broken bonds cannot be calculated with Equation 6. Therefore, in the coupling model proposed herein, we use: where V fA and V fB are the volumes of the fluids in different domains, and V 0A and V 0B are the volumes of the fluids in the different domains under zero hydraulic pressure, all in m 2 . Figure 4 shows the modeling procedure used to simulate the hydromechanical coupling in PFC 2D . The left side of the figure shows the equations that describe the interactions between particles in shale reservoirs after every time step, and the right side describes the fluid migration in the rock reservoirs. The tensile and shear fractures can be determined by using Equations 7 and 8: where R refers to the average of the radiuses of the two particles, in m; max and max are the tensile and shear strength, in MPa, respectively; and A, I, and J are polar moment of inertia of the parallel-bond cross section. At the same time, these fracture positions (in the rock matrix or BPs) can also be determined based on the position at which the HF was initiated. Therefore, in an independently developed subroutine, we defined four types of HFs: tensile hydraulic fractures of a shale reservoir (THFSR), shear hydraulic fractures of a shale reservoir (SHFSR), tensile hydraulic fractures of a BP (THFBP), and shear hydraulic fractures of a BP (SHFBP).

| PARAMETER CALIBRATION
The mechanical parameters for rock were calibrated by using laboratory-scale experiments in which shale specimens were subjected to confining pressure. Next, based on the calibrated parameters, we developed HF models with different conditions of rock brittleness B. Finally, the resulting models were applied, and the breakdown pressure induced by HF was compared with theoretical calculations.

| Verification of rock mass parameters
In the PFC 2D , the macroparameters of the simulation models must be calculated based on the microparameters of the particles. These macroparameters are usually obtained from laboratory-scale experiments and from in situ field data. 8,28 Natural fractures in the rock mass degrade strength and necessitate a larger model. However, the direct application of microparameters calibrated from the laboratory and field data will result in an inappropriate increase in the strength of the modeled rock mass. Therefore, the size effect must be considered when calibrating the parameters used in the model.
In this work, we investigated specimens with vertical BPs (SVBP) and specimens with horizontal BPs (SHBP). The parameters that describe the rock mechanics in the DEM were We tested 20 SVBP specimens and 22 SHBP specimens under six different confining pressures: 5, 10, 15, 20, 25, and 30 MPa. See Ref. 29 for the test procedure and results. Figure 5 shows the data from these experiments.
The results shown in Figure 5 reveal that, within a specimen group (SVBP or SHBP), Young's modulus did not fluctuate significantly under the different pressure conditions. However, the layer orientation clearly affected Young's modulus: The average Young's modulus for SHBP was 28.45 GPa, which was far greater than that for SVBP, at 21.47 GPa. In addition, because of the different layer orientations in the shale specimens, the cohesion of SHBP (37.05 MPa) exceeded that of SVBP (31.25 MPa), and the internal friction angle of SHBP (34.02°) exceeded that of SVBP (29.24°).
Both the matrix and BP of the shale reservoir model must be calibrated before use. In the present work, the calibration was conducted by trial-and-error tests and a sensitivity analysis to match the numerical simulation results with the experimental data. The parameters obtained from the proposed model, including the cohesion and internal friction angle, were calibrated by using the mechanical parameters obtained under different confining pressures. Young's moduli for SVBP and SHBP changed little in response to the different confining pressures, at not less than 21.47 and 28.45 GPa, respectively. The parameters of the BPs were then calibrated by matching the results of the numerical simulation with the experimental data. During this calibration process, the parameters of the rock matrix and BPs were tested several times. The specific calibration procedures and methods for obtaining the parameters of shale rock containing BPs are detailed in our previous work. 29 Figures 6 and 7; Table 3 compare the results of the numerical simulations with experimental data. The results show that after the particle parameters were calibrated, the mechanical parameters and failure modes obtained from the numerical simulation matched well with those observed experimentally under different confining pressures. Tables 4 and 5 list the calibrated parameters of the shale rock and BPs.

| Models of rock mass with differing brittleness
The rock brittleness B is a crucial factor in determining whether an unconventional tight reservoir is worth fracturing. Numerous methods have been proposed to identify the rock brittleness. The most commonly applied is Rickman's method, 20 in which Young's modulus E and Poisson's ratio of the rock mass are used to determine the rock brittleness B, as follows.
where, E max and E min are the maximum and minimum Young's modulus, respectively, and max and min are the maximum and minimum Poisson's ratio, respectively.
From Equations 9-11, an increased Young's modulus and reduced Poisson's ratio will result in a high brittleness index, and an increased Young's modulus clearly reduces the likelihood of rock failure. The relationship between the rock brittleness index and the ratio of normal stiffness to shear stiffness is established by Young's modulus E and Poisson's ratio of the rock mass. The sensitivity analysis of the parameters calibrated with the simulation model revealed that, holding the other parameters constant, an increase in the ratio of normal to shear stiffness left the model strength unchanged. However, such increases can lead to an increased Young's modulus E and also a decreased Poisson's ratio ( Figure 8). Hence, an increase in the ratio of normal to shear stiffness can improve the brittleness index. This further confirms the suitability of Rickman's brittleness formula for this study, in that the elastic modulus and Poisson's ratio used in the formula are directly related to the ratio of normal to shear stiffness.
According to Rickman, in Equations 9 and 10, E max = 58.6 GPa, E min = 5.5 GPa, max = 0.37, and min = 0.16. To simplify the calculation, Young's moduli were converted at 1 F I G U R E 6 Peak strength observed experimentally compared with that obtained by simulation 29 psi = 6.895 kPa, which was determined based on the dynamic values obtained from the P-and S-polarized waves of log data acquired from boreholes in the Barnett shale reservoir. 30 In this work, the ratio of normal to shear stiffness of particles and that of parallel bonds were adjusted based on the calibrated parameters in conjunction with various E and (but the same E max , E min , max , and min ) to obtain different rock brittleness values corresponding to diverse working conditions (Table 6). Six DEM models were thus developed: Three SVBP models had brittleness values of 26, 35, and 59%, respectively, and three SHBP models had brittleness values of 39, 49, and 70%, respectively, as calculated from Equations 9-11.
The macroparameters obtained from the DEM depend on the ratio of total model size to average particle radius R m .   When R m > 50, the effect of particle radius variation on the model parameters is extremely limited, and these macroparameters tend to remain stable. 31 The model established in this study was 100 × 100 mm 2 , and the particle radius and the ratio of total model size to average particle radius were in the ranges of 0.3-0.5 mm and 200-333, respectively. Therefore, the parameters obtained from the DEM could be considered reliable. Figure 9 illustrates the size and loading condition applied to the shale reservoir model of HF. The model was developed as follows. A wall with infinite stiffness was generated to surround a 100 × 100 mm 2 area, and the PBM was used to bond together randomly distributed particles. The particleto-particle contact was replaced with the SJM. A 5-mm-diameter borehole was established at the center of the model for injecting the fracturing fluid, and the maximum stress S max and minimum stress S min were applied in the vertical and horizontal directions, respectively. The quantities S max and S min used in the shale reservoir model were calculated by using:

| Setting up the model and verifying the hydraulic parameters
where z is the depth in km. In this work, we used z = 1.5 km, which is the depth of the shale rock in Pengshui County, China. Therefore, according to Equations 12 and 13, the S max and S min in this work were 25.6 and 15.6 MPa, respectively.
Because of the size effect, the calibrated microparameters could not be applied to the shale reservoir containing BPs. Therefore, we used the same scale from the experiment in this model. The large ratio (>200) of the model size to the particle size in this work revealed clearly how the HFN was formed. The model simulated a continuously injected fluid (water) into the modeled borehole until an HF reached a model border, at which point the duration of the injection was recorded, and then, this duration was applied to all subsequent operations.
To further verify the reliability of the model, the breakdown pressure obtained from a numerical simulation of an isotropic HF model (without BPs) was compared with the calculation results. Haimson   following regression equation for calculating breakdown pressure: The calibrated parameters discussed above were used in the numerical simulation to obtain the breakdown pressure, and the initial pore pressure P 0 was set to zero. In addition, the tensile strength of the isotropic shale reservoir was set at 8.51 MPa; S max = 20 or 10 MPa, and S min ranges from 5 to 10 MPa. The hydraulic properties were tuned based on repeated tests. The breakdown pressures obtained from the numerical simulation and from the experiment revealed that the errors between the numerical simulation and the theoretical value were acceptable ( Figure 10). Table 7 lists the verified hydraulic properties.

SVBP AND SHBP CASES
During hydraulic fracturing, the rock brittleness B strongly affects HFN formation, and the interaction between HFs and the BP further complicates this interaction. This section investigates the formation of the HFN as a function of rock brittleness by considering injection pressure, type of HF, HF evolution, fluid pressure, seepage area, and other factors. Injection pressure is the average of the pressure measured around the injection hole by measurement circles. Figure 11 shows the injection pressure as a function of B. Based on this curve, the fluid injection can be divided into the following three stages. Stage I is the pressure-increase stage, in which the injection pressure increases from zero to the breakdown pressure in 2 seconds. In the first second, B exerts a relatively small effect on the injection pressure, although the breakdown pressures of shale reservoirs with different B differ from each other. A larger B requires a larger breakdown pressure, which was especially evident in SHBPs with B = 70%, where the breakdown pressure increased to 44.5 MPa. For the same B (B = 35%), the breakdown pressure of SHBP (29.7 MPa) was greater than that of SVBP (22.1 MPa), indicating that the BP orientation also affects the breakdown pressure.
In Stage II of the injection process, the injection pressure decreased by less than 21% when B was small and decreased by more than 35% when B was relatively large. Stage II lasted approximately 60% of the fracturing time. In Stage III, the injection pressure stabilized and remained constant independent of B.
The HF type and distribution are two key factors for evaluating HFN propagation. Figure 12 shows the formation of HFN for different B, and Figure 13 shows the relationship between injection time and the accumulated number of HF types. Because both S max and the SVBP have vertical orientations, HF in the SVBP generally propagated in the S max direction, as shown in Figure 12A-C. However, the BP caused the HFs to coalesce with each other during propagation, and therefore, a few HFs propagated in the S min direction. However, in the SHBP, the HF propagation was affected by both the vertical S max and the horizontal BP, resulting in HFs that radiated in all directions, as shown in Figure 12D-F. Therefore, the HF propagation is influenced not only by S max but also by the BP orientation.
For all rock brittleness B and all BP orientations, the most numerous type of HFs was THFSR. Although the number of instances of SHFBP exceeded that of THFSR at the beginning of the injection process, the latter surpassed the former  Figure 13C. During the entire injection process, the numbers of SHFSR and THFBP were small and grew slowly. Therefore, in this work, we mainly explored how THFSR and SHFBP evolved. The presence of the BP and different values of rock brittleness B resulted in large differences in the HF evolution. In the SVBP, increasing B from 26% to 59% resulted in THFSR decreases from over 2000 to approximately 1100, with SHFBP dropping to fewer than 300 ( Figure 13A,C). In the SHBP, a similar evolution occurred when B increased from 35% to 70% ( Figure 13D-F). Therefore, in all conditions, the growth of B can result in a significant decrease in THFSR and a small decrease in SHFBP, with SHFSR and THFBP numbers remaining essentially unchanged.
In addition to the HF distribution and evolution, the fluid pressure distribution given by the model is also relevant to the analysis of how B affects the HFN. Figure 14 shows the fluid pressure distribution in shale reservoirs with different values of B. In SVBPs, the average fluid pressure increased as B increased from 26% to 59%, and the fluid pressure was distributed along the S max direction ( Figure 14A,C). At the same time, because the BP was parallel to S max , the BP could be more easily activated when subjected to the same fluid pressure. The horizontal stresses acting on the vertical BP combined with each other (Figure 14D-F).
The SHBP case differed significantly from that of the SVBP. When B was relatively small (B = 35%), the BP was orientated perpendicular to S max , and the fluid pressure propagation was affected by the orientation of both S max and the BP. Therefore, the area with high fluid pressure was oriented at approximately 42° with respect to S max . With B increased to 49%, this angle was significantly reduced to approximately 16°. However, when B increased to the maximum of 70%, the fluid pressure was distributed such that areas with higher fluid pressure were almost parallel to S max . These results are explained by the shale reservoir having a small Young's modulus E and a large Poisson ratio when B is small, which allows the propagation of HF and fluid pressure in the shale reservoir to be more easily affected by the BP. Seepage area is another important parameter in HFN analysis, and Figure 15 shows how B was found to affect the seepage area in this study. Regardless of the orientation of the BP, the largest seepage area was always around the wellbore. With increasing B, the size of the area in which significant particle displacement occurred (greater than 5.00 × 10 −4 m) gradually decreased. In the SVBP with B = 59%, the BP was parallel to S max and the shale reservoir was stronger, and therefore, the fracturing fluid propagated to the borders of the model before generating a large seepage area ( Figure 15C). In the SHBPs, the fluid flow was relatively large because the direction of the BP was perpendicular to that of S max , and the HFs coalesced with each other in the normal direction at the ends of the BP, facilitating the flow of the fracturing fluid in the horizontal direction. For the reservoir models with different B, the particles in the model barely moved, and the fracturing fluid did not flow through the middle sections of the left and right sides (see FF-1 and FF-2 in Figure 15F) because the HFs mainly propagated in the direction of S max .

HFN PROPAGATION
In addition to E and , the angle and density of the BP also strongly influence the rock brittleness B. Therefore, this section presents the sensitivity analysis conducted to determine the sensitivity of HFN propagation to the angle and density of the BP. In this analysis, the parameters of the shale reservoir model were taken from Tables 4, 5, and 7.

| Impact of BP angle
The angle between the BP and the vertical is denoted as or the "BP angle." For a constant ratio of normal stiffness to shear stiffness of particles bonded in parallel, a change in the BP angle will cause fluctuations in the rock brittleness B. Figure 16 shows the rock brittleness, breakdown pressure, and total HF variations in response to .
As the BP angle increased from 0° to 90°, B first decreased and then increased. When = 30°, B dropped to a minimum (24.1%), whereas when increased to 90°, B increased to over 39% (almost twice the minimum value). The breakdown pressure of the shale reservoir followed the same trend. From = 30° to 90°, the breakdown pressure increased by a factor of approximately 1.8 because the value of B (an intrinsic property) for a shale reservoir with BPs is directly related to E and , and the breakdown pressure correlates positively with the tensile strength. These results were consistent with those of an experimental investigation, 33 where as increased, E and the tensile strength first decreased and then increased, whereas followed the opposite trend, reaching a minimum at = 30° and a maximum at = 90°.
In this study, as the BP angle increased, the quantity of HFNs generated first increased and then decreased after peaking at = 15° and 30°. Figure 17 shows how the quantity and percent of each type of HF were related to . Although the total quantity of HFs in the shale reservoirs varied strongly with , this variation was mainly caused by the variation in the quantity of THFSR. The number of the other three types of HFs remained almost constant. In particular, the number and percent of SHFSR and THFBP remained small for all BP angles. Because the THFSR was the most numerous, the percent of THFSR also remained essentially constant (>60%). Figure 18 shows the HF distribution, fluid pressure, and seepage area for different BP angles . Under the influence of After the BP was activated, mutual penetration occurred between HFs in the direction perpendicular to the BP because of the influence of S max . At = 0°-45°, the fluid pressure was less affected by the BP, and the area with higher fluid pressure expanded parallel to the orientation of S max . For ≥ 60°, the area with higher fluid pressure expanded at an angle to the orientation of S max . For the same , a smaller B led to a larger angle, as discussed in Section 4. At = 0°-30°, the areas in the model with larger particle flows (>5.00 × 10 −4 m) were close to each other. In the direction of S max , the areas with larger fluid flows connected with each other, resulting in better permeability of the shale reservoir. With increasing ( = 45°-60°), the size of the area with the larger particle flow decreased dramatically, and the particles flowed mainly through the areas around the injection hole. At very high BP angles ( = 75°-90°), the areas with larger particle flows started to expand into other areas and were not limited to the area around the injection hole. However, the permeability of the reservoir remained poorer than with small . Regardless of the BP angle , the seepage area always expanded in the direction of S max .

| Impact of BP density
In addition to the BP angle , the BP density also strongly affects the rock brittleness, leading to differences in HFN formation in the shale reservoir models. Figure 19 shows the rock brittleness, breakdown pressure, and total number of HFs as a function of BP density. These results indicated that increased BP density dramatically decreased rock strength, as reflected by the decrease in Young's modulus and the increase in Poisson's ratio. B and breakdown pressure decreased logarithmically with increasing BP density. When the BP density exceeded 1.24 m −1 , the BP in the shale reservoir was fully developed and the breakdown pressure gradually decreased to S min . In addition, the cumulative number of HFs increased exponentially with increasing BP density. Figure 20 shows the various types of HF as a function of BP density. Although the total number of HFs increased exponentially with increasing BP density, the number of THFSR remained almost constant, or even decreased slightly. However, the number of SHFBPs increased exponentially because with increasing BP density, the initiation and propagation of HFs intensified and many HFs started to interact with the BP before propagating into the matrix. When the BP density was 1.86 m −1 , there was slightly fewer THFSR than SHFBP, and a sharp increase in the number of SHFBP caused the percent of THFSR to decrease exponentially. In addition, the numbers and percentages of SHFSR and THFBP were small and remained essentially constant. Figure 21 shows the distribution of HF, fluid pressure, and seepage area as a function of BP density. With small BP density ( = 0.31 or 0.62 m −1 ), HFs propagated almost throughout the shale matrix, producing a clear HFN that consisted mainly of THFSRs. During the propagation, only a small number of SHFBP became active, and because the shale reservoir was at maximum strength, the fluid pressure was generally high. Fluid flowed mainly through the seepage areas in the HFN and very little in other areas. No fluid flowed through the areas far from the injection hole, on either side of the model, because HFs in those areas did not fully coalesce with BPs.
With increasing BP density ( ≥ 0.93 m −1 ), the quantity of activated BPs in the HFN gradually increased. Over the injection time, most fracturing fluid flowed into the BP, resulting in a gradual decrease in the number of HFs in the shale matrix with increasing BP density. The fluid pressure and seepage area also became increasingly concentrated at the central axis around the injection hole and parallel to S max , and the average fluid pressure also dropped gradually. These results occurred because increasing BP density degraded the strength of the shale reservoir and caused a large quantity of fracturing fluid to flow gradually to peripheral areas around the injection hole. When the BP density reached its maximum value of 1.86 m −1 , the HFs propagated almost simultaneously in both the BPs and the shale matrix. Therefore, the number of SHFBP approached that of THFSR. Although the fluid pressure was mainly concentrated around the central axis of the injection hole, the aperture through which the fracturing fluid flowed was much larger than when = 0.31.

| CONCLUSIONS
In this work, we used the rock brittleness B, which is based on Young's modulus and Poisson's ratio, to construct six shale reservoir models containing horizontal and vertical BPs with varying brittleness. The parameters of the matrix and the BPs in the shale reservoir were calibrated with respect to experimental data, and the reliability of the fluid parameters of the model was verified by comparing the breakdown pressure obtained from the models with that obtained from theoretical calculations. The six shale reservoir models were then used to determine how rock brittleness affects HFNs based on several factors, including the distribution of breakdown pressure and fluid pressure, the distribution and number of HFs, and the seepage area. We also analyzed the sensitivity of the HFN formation to the angle and density of BPs because these factors also affect the rock brittleness. The results motivated the following conclusions.
1. The shale reservoir models showed three stages of HFN formation. The breakdown pressure and the magnitude of the fluid pressure increased with rock brittleness. In the SHBP with B = 70%, the breakdown pressure reached approximately 44.5 MPa. Although the number of THFSR decreased as B increased, the THFSR always accounted for the largest percent of HF types, whereas the number and percent of SHFSR and THFBP remained almost constant. The number of THFSR was highest, at over 2000, in shale reservoirs with B = 26%; this was 2.5 times the sum of the other three HF types. 2. During the injection process, the direction in which the area with higher fluid pressure was distributed varied with changes in the angle between the BP direction and the direction of S max , with the distribution directed parallel to S max with larger B, and deviating from S max with smaller B.
For the SHBP with the minimum B (B = 35%), the distribution direction changed by 42°. The size of the area with high fluid flow decreased with increasing B. Regardless of the BP orientation, the seepage area was always largest around the wellbore. 3. With increasing BP angle , B and the breakdown pressure first decreased and then increased. The number of HFs generated was minimal at = 15° and maximal at = 30°, with the latter case resulting in better permeability of the shale reservoir. Regardless of the BP angle, the percentages of the various types of HFs remained essentially constant, and the seepage areas always extended in the direction of S max . 4. The rock brittleness B and the breakdown pressure gradually decreased with increasing BP density . During this process, although the total quantity of HF increased gradually, the number of THFSR was no longer the largest. At = 1.83 m −1 , the numbers of SHFBP and of THFSR were almost the same. In addition, with increasing , the fluid pressure gradually concentrated around the central axis, and the areas with larger fluid seepage gradually connected with each other.