A Case Study on the Optimal Design of the Horizontal Wellbore Trajectory for Hydraulic Fracturing in Nong’an Oil Shale

A horizontal well with hydraulic fractures is key to forming a fracture network in oil shale for the generated hydrocarbon flows. By considering the influence of anisotropic strength, a prediction model is proposed for fracture initiation by studying different fracture initiation modes (FIMs) in oil shale: failure of the intact rock matrix and of the bedding planes. Through a case study on Nong’an oil shale, the influences of wellbore trajectory and bedding planes on the fracture initiation pressure (FIP), location (FIL), and FIM were analyzed and the induced changes in wellbore trajectory design were concluded. The preferred angle between the wellbore axis and the minimum horizontal principal stress was the same of 90◦ or 270◦, when the lowest required FIP corresponded to the failure of the intact rock matrix. However, when the angle corresponded to the failure of the bedding planes, the preferred direction of the wellbore axis was away from the fixed direction and not corresponding to the lowest required FIP due to the fracture morphology. The error between the theoretical and experimental results ranges from 7% to 9%. This research provides a framework for the design of horizontal wellbore trajectories in oil shale, for easier fracture initiation and more complex fracture networks.


Introduction
Hydraulic fracturing in the horizontal well is a promising approach to form the fracture network for the generated hydrocarbon flow in the in-situ retorting process [1]. Since the wellbore trajectory and many other factors, such as the in-situ stress regime, rock fabric and its physical properties, perforation parameters and so on [2][3][4][5][6][7], affect the fracture initiation and patterns together, and eventually influence the stimulated reservoir volume (SRV), the change of wellbore trajectory design due to the different fracture initiation mechanism in oil shale should be worked out. For oil shale, the shear slippage along the bedding planes and induced bedding fractures [8][9][10][11][12] indicates that its anisotropic strength contributes mostly to the fracture initiation, and thus changes the design of wellbore trajectory. bedding planes' dip angle and direction on the fracture initiation was concluded, and the FIP, FIL, FIM and the initial fracture morphology corresponding to different wellbore trajectories were confirmed. In addition, experiments were carried out on Nong'an oil shale to validate the accuracy of this fracture initiation model.

Stress State Around a Horizontal Wellbore in Oil Shale
Oil shales in Jilin Province, China are characterised by shallow buried depths and small thicknesses. The horizontal well is superior on account of its more workable seam and is thus conducive to lowering the production cost. Consequently, this research has focused on the stress state around the horizontal wellbore to analyse hydraulic fracture initiation in horizontal wells of oil shales. With the anisotropic strength of oil shale being induced by its bedding planes, it is first necessary to establish a hydraulic fracturing model of a horizontal well in oil shale. The simplified hydraulic fracturing model is shown in Figure 1, which is developed by ignoring the existence of natural fractures and other factors, and assuming that bedding planes with the same spacing, dip direction and dip angle are distributed in the oil shale. The X-Y-Z coordinate system in this figure refers to the in-situ stress coordinate system, in which (along X) represents the minimum horizontal principal stress, (along Y) represents the maximum horizontal principal stress and (along Z) represents the overburden pressure.

Failure Mechanism of the Intact Rock Matrix
This study has analysed the failure of oil shales through a comparison of the failures of the intact rock matrix and bedding planes. Through such comparative analysis, the change in the FIM of oil shale with horizontal wellbore trajectory (the drilling direction of a horizontal well) under different bedding planes (dip direction and dip angle) was determined, as well as the FIP and FIL. Consequently, different failure criteria were chosen for different FIMs in the computational formulas for different fracture initiations in oil shale.

Failure Criteria
Based on the lithology of oil shales and previous researches, the appropriate failure criterion for the intact rock matrix was determined. Colmenares and Zoback [37] verified the degree of fitting between the previous polyaxial test data (for Dunham dolomite, Solenhofen limestone, Shirahama Sandstone, Yuubari shale and Amphibolite from the KTB site) and seven different failure criteria, which included Mohr-Coulomb criterion, Mogi-Coulomb criterion, Drucker-Prager failure criterion and so on. Al-Ajmi and Zimmerman [38] estimated the collapse pressure and fracture pressure of Dunham Dolomite and Solenhofen Limestone with Mohr-Coulomb criterion, Mogi-Coulomb criterion and Drucker-Prager failure criterion and carried out a comparative analysis between the

Failure Mechanism of the Intact Rock Matrix
This study has analysed the failure of oil shales through a comparison of the failures of the intact rock matrix and bedding planes. Through such comparative analysis, the change in the FIM of oil shale with horizontal wellbore trajectory (the drilling direction of a horizontal well) under different bedding planes (dip direction and dip angle) was determined, as well as the FIP and FIL. Consequently, different failure criteria were chosen for different FIMs in the computational formulas for different fracture initiations in oil shale.

Failure Criteria
Based on the lithology of oil shales and previous researches, the appropriate failure criterion for the intact rock matrix was determined. Colmenares and Zoback [37] verified the degree of fitting between the previous polyaxial test data (for Dunham dolomite, Solenhofen limestone, Shirahama Sandstone, Yuubari shale and Amphibolite from the KTB site) and seven different failure criteria, which included Mohr-Coulomb criterion, Mogi-Coulomb criterion, Drucker-Prager failure criterion and so on. Al-Ajmi and Zimmerman [38] estimated the collapse pressure and fracture pressure of Dunham Dolomite and Solenhofen Limestone with Mohr-Coulomb criterion, Mogi-Coulomb criterion and Drucker-Prager failure criterion and carried out a comparative analysis between the calculation results and experimental data. Hashemi et al. [39] also compared the calculation results of the three failure Energies 2020, 13, 286 4 of 26 criteria with the results of a self-designed thick-walled hollow cylinder test. According to the above researches, the Mogi-Coulomb criterion proposed by Al-Ajmi and Zimmerman [40] was chosen for the failure of the intact rock matrix, as expressed by Equation (1): where τ oct is the octahedral shear stress and σ m,2 is the mean normal stress, which are defined as: where the parameters a and b are defined as: in which c O and ϕ O are the cohesion and the angle of internal friction of the intact rock matrix, respectively, and σ 1 , σ 2 , and σ 3 represent the first, second and third principal stresses, respectively.

Stress Analysis
In order to analyse the hydraulic fracture initiation of the intact rock matrix, it is first necessary to confirm the stress distribution around the horizontal wellbore. Therefore, we introduced the in-situ stress coordinate system (ICS) and horizontal borehole coordinate system (HCS) to calculate the stress components for the failure criterion of the intact rock matrix. It is assumed that σ h , σ H , σ V are along the North, East and down directions of the global coordinate system (GCS), respectively; thus, we have ignored the difference between the GCS and the ICS for simplified calculation.
In order to determine the stress state around the horizontal wellbore (σ HCS ), the stress state in the ICS is transformed into that in the HCS, as shown in Figure 2. It is assumed that the ICS is defined such that the X axis is along the direction of the minimum horizontal principal stress (σ h ), the Y axis is along the direction of the maximum horizontal principal stress (σ H ) and the Z axis is along the direction of the overburden pressure (σ V ). With the angle of 90 • being measured from the Z axis to the Z H axis, which is aligned with the axis of the wellbore, and the angle α H being measured between the X axis and the projection of the X H axis on the X-Y plane, the configuration of the wellbore with respect to the ICS has been defined. Here, α H refers to the angle between the drilling direction and the minimum horizontal principal stress, and also represents the wellbore trajectory of the horizontal well. The rotation matrix for the transformation of the stress components from the ICS to the HCS is derived by rotating the axes (X-Y-Z) of the ICS by an angle of α H about the Z axis and by an angle of 90 • about the Y H axis, as expressed in the following formula: Energies 2020, 13, 286 5 of 26 the X axis and the projection of the XH axis on the X-Y plane, the configuration of the wellbore with respect to the ICS has been defined. Here, refers to the angle between the drilling direction and the minimum horizontal principal stress, and also represents the wellbore trajectory of the horizontal well. The rotation matrix for the transformation of the stress components from the ICS to the HCS is derived by rotating the axes (X-Y-Z) of the ICS by an angle of about the Z axis and by an angle of 90° about the YH axis, as expressed in the following formula:  Figure 2 shows the relationship between the ICS and the HCS, which is used to calculate the stress distribution around a horizontal wellbore. The stress components acting on the horizontal wellbore in the polar coordinate system are calculated using Equation (7), which is stated by Bradely [41] for elastic and homogeneous formations: where σ rr , σ θθ and σ zz are the radial, tangential and axial stresses, respectively; τ θz , τ rθ and τ rz are the shear stresses; P m is the injected hydraulic pressure; P p is the formation pore-pressure; α is Biot's parameter of the formation; and ν is Poisson's ratio of the rock. By substituting Equation (6), the stress components acting on the horizontal wellbore can be expressed by: By calculating the eigenvalues of σ PCS , the principal stresses σ 1 , σ 2 and σ 3 acting at any point on the wellbore have already been confirmed. Further, Zhang et al. [9] and Wang [10] showed that the oil shale in Jilin Province was under the normal faulting stress regime. As this research mainly investigates the Nong'an oil shale in Jilin, σ 1 , σ 2 , and σ 3 can be expressed by: Energies 2020, 13, 286 6 of 26 By using the principal stresses calculated from Equation (9) and the Mogi-Coulomb failure criterion presented in Equation (1), the fracture initiation from the wellbore of the intact rock matrix can be evaluated.

Failure Criteria
Two failure modes of oil shales are failure of the intact rock matrix and failure of the bedding planes. Therefore, it is necessary to analyse the failure mechanism of the bedding planes. Failure of the bedding planes can be divided into shear slippage along the bedding planes and tensile failure across the bedding planes. The Mohr-Coulomb criterion was selected to describe the shear slippage along the bedding planes, as expressed by Equation (10): where τ and σ n are the shear stress and normal stress projected onto the surface of the bedding planes, c B is the shear strength or cohesion of the bedding planes, and µ B is the coefficient of internal friction of the bedding planes, where: where ϕ B is the internal friction angle of the bedding planes, and τ B xy , τ B xz are the two shear stress components obtained by projecting onto the surface of the bedding planes.
Meanwhile, the tensile failure criterion was selected for the tensile failure across the bedding planes, which meant that tensile failure occurred when the normal stress on the bedding planes was larger than its tensile strength, as expressed in Equation (12): where S t is the tensile strength of the bedding planes.

Stress Analysis
For fractures to initiate from the horizontal wellbore when failure of the bedding planes occurs, it is necessary to transform the stresses acting on the wellbore to stress components which are projected onto the surface of the bedding planes. According to the coordinate system settings of Lee et al. [42], we introduced the bedding plane coordinate system (BCS) and polar coordinate system (PCS). This paper assumes that the GCS completely coincides with the ICS, which makes no difference to the calculation of the stress components acting on the bedding planes [42]. Figure 3 was drawn to show the transformation of the stress components acting on the wellbore in the PCS to the stress components in the corresponding rectangular HCS.
In Section 2.1.2, Figure 2 shows the relationship between the HCS and the ICS. For the stress components acting on the bedding planes, the relationship of the stress components between the BCS and the ICS has to be established, as shown in Figure 4. With the angle (0 − 360°) being measured from the X axis to the projection of the XB axis on the horizontal plane and the angle being measured from the XB axis to the horizontal plane, the normal direction in line with the XB axis has been confirmed, as well as the bedding planes. The rotation matrix for the transformation of the stress components from the ICS to the BCS is derived by rotating the axes (X-Y-Z) of the ICS by an angle of about the Z axis and then rotating the axes (X'-YB-Z) by an angle of about the YB axis, which can be expressed by:  In Section 2.1.2, Figure 2 shows the relationship between the HCS and the ICS. For the stress components acting on the bedding planes, the relationship of the stress components between the BCS and the ICS has to be established, as shown in Figure 4. With the angle α B (0 − 360 • ) being measured from the X axis to the projection of the X B axis on the horizontal plane and the angle β B being measured from the X B axis to the horizontal plane, the normal direction in line with the X B axis has been confirmed, as well as the bedding planes. The rotation matrix for the transformation of the stress components from the ICS to the BCS is derived by rotating the axes (X-Y-Z) of the ICS by an angle of α B about the Z axis and then rotating the axes (X'-Y B -Z) by an angle of β B about the Y B axis, which can be expressed by: Energies 2020, 13, 286 7 of 26 Figure 3. The stress components acting on the wellbore.
In Section 2.1.2, Figure 2 shows the relationship between the HCS and the ICS. For the stress components acting on the bedding planes, the relationship of the stress components between the BCS and the ICS has to be established, as shown in Figure 4. With the angle (0 − 360°) being measured from the X axis to the projection of the XB axis on the horizontal plane and the angle being measured from the XB axis to the horizontal plane, the normal direction in line with the XB axis has been confirmed, as well as the bedding planes. The rotation matrix for the transformation of the stress components from the ICS to the BCS is derived by rotating the axes (X-Y-Z) of the ICS by an angle of about the Z axis and then rotating the axes (X'-YB-Z) by an angle of about the YB axis, which can be expressed by:  From Equations (6), (14), and (16), the shear stress and normal stress projected onto the surface of the bedding planes can be expressed in the form of Equations (17) and (18): Then, by using the failure criteria for shear failure along the bedding planes and tensile failure across the bedding planes, the fracture initiation of the bedding planes from the wellbore can be evaluated.

Hydraulic Fracture Initiation Analysis of a Horizontal Well in Nong'an Oil Shale
It is well-known that oil shales in Jilin Province show the prominent feature of shallow buried depth, such as the oil shale 70 m underground in Nong'an and the oil shale 400 m underground in Fu'yu. From 2013 to 2015, Jilin University carried out geological investigation, hydraulic fracturing and in-situ exploration pilot tests on Nong'an oil shale. Based on these researches, the conditions of occurrence and mechanical parameters of Nong'an oil shale have been obtained, as shown in Table 1 [9,10]. By considering the fracture initiation model developed in the Section 2, an example checking calculation about Nong'an oil shale was conducted, to record the FIP and FIL when failure of the intact rock matrix or bedding planes occurred. By setting different bedding planes (different α B and β B ), the influence of bedding planes on the hydraulic fracture initiation, which included the changes in the FIP, FIL, and FIM, was concluded, and the optimum horizontal wellbore trajectory for lower FIP and better initial fracture morphology was discussed.

Fracture Initiation of the Intact Rock Matrix
Combining Equations (1) and (9), the FIP at every point of a horizontal wellbore can be calculated when failure of the intact rock matrix occurs; this can be used to analyse the changes in the FIP and FIL as the drilling direction (α H ) varies. As can be seen from Equation (8), the stress distribution on the wellbore has a direct relation with the angle θ, which represents the location on the wellbore. At the point θ = 0 • or 180 • on the wellbore, the tangential and axial stresses reach the minimum values, whereas they are maximum at the point θ = 90 • or −90 • . Consequently, with the increase in the injected hydraulic pressure (P m ), σ θθ (σ θθ < 0) decreases gradually, and eventually exceeds the ultimate strength of the rock. Thus, failure of the intact rock matrix is observed in the form of hydraulic fracturing.
By analysing the formula used for predicting the failure of the intact rock matrix, the angle (α H ) between the drilling direction of the horizontal well and the minimum horizontal principal stress was found to be the main factor, which was irrelevant in the case of the bedding planes. Consequently, Figure 5 was drawn with the FIP at every point on the wellbore under different drilling directions when failure of the intact rock matrix occurred.
Energies 2020, 13, 286 9 of 26 By considering the fracture initiation model developed in the Section 2, an example checking calculation about Nong'an oil shale was conducted, to record the FIP and FIL when failure of the intact rock matrix or bedding planes occurred. By setting different bedding planes (different and ), the influence of bedding planes on the hydraulic fracture initiation, which included the changes in the FIP, FIL, and FIM, was concluded, and the optimum horizontal wellbore trajectory for lower FIP and better initial fracture morphology was discussed.

Fracture Initiation of the Intact Rock Matrix
Combining Equations (1) and (9), the FIP at every point of a horizontal wellbore can be calculated when failure of the intact rock matrix occurs; this can be used to analyse the changes in the FIP and FIL as the drilling direction ( ) varies. As can be seen from Equation (8), the stress distribution on the wellbore has a direct relation with the angle , which represents the location on the wellbore. At the point = 0° or 180° on the wellbore, the tangential and axial stresses reach the minimum values, whereas they are maximum at the point = 90° or −90°. Consequently, with the increase in the injected hydraulic pressure ( ), ( < 0) decreases gradually, and eventually exceeds the ultimate strength of the rock. Thus, failure of the intact rock matrix is observed in the form of hydraulic fracturing. By analysing the formula used for predicting the failure of the intact rock matrix, the angle ( ) between the drilling direction of the horizontal well and the minimum horizontal principal stress was found to be the main factor, which was irrelevant in the case of the bedding planes. Consequently, Figure 5 was drawn with the FIP at every point on the wellbore under different drilling directions when failure of the intact rock matrix occurred. In Figure 5, the points (r, ) in the curves correspond to the point on the wellbore at which the FIP of the intact rock matrix is equal to r MPa. The FIPs on the wellbore are minimum at the point = 0° or 180°, irrespective of the value of . In other words, the change in the drilling direction only influenced the FIP, and not the FIL, as the FIL was at = 0° or 180° on the wellbore. When | − 90°× | ( = 1,3) was kept the same, which meant that the angle between the drilling direction and the minimum horizontal principal stress was constant, the FIP curves on the wellbore were identical. With the decrease in | − 90°× |, the FIP at = 0° or 180° decreased, while the fracture initiation region increased, as shown in Figure 5.
Because the FIL of the intact rock matrix was at the point = 0° or 180°, Figure 6 was drawn for the FIP at = 0° or 180° for different drilling directions. This figure shows that the FIP of the intact rock matrix reaches the maximum value ( ) for = 0° or 180° and the minimum value ( ) for = 90° or 270°. Consequently, with the decrease in the angle between the In Figure 5, the points (r, θ) in the curves correspond to the point θ on the wellbore at which the FIP of the intact rock matrix is equal to r MPa. The FIPs on the wellbore are minimum at the point θ = 0 • or 180 • , irrespective of the value of α H . In other words, the change in the drilling direction only influenced the FIP, and not the FIL, as the FIL was at θ = 0 • or 180 • on the wellbore. When |α H − 90 • × i| (i = 1, 3) was kept the same, which meant that the angle between the drilling direction and the minimum horizontal principal stress was constant, the FIP curves on the wellbore were identical. With the decrease in |α H − 90 • × i|, the FIP at θ = 0 • or 180 • decreased, while the fracture initiation region increased, as shown in Figure 5.
Because the FIL of the intact rock matrix was at the point θ = 0 • or 180 • , Figure 6 was drawn for the FIP at θ = 0 • or 180 • for different drilling directions. This figure shows that the FIP of the intact rock matrix reaches the maximum value (P Rm f max ) for α Hmax = 0 • or 180 • and the minimum value (P Rm f min ) for α Hmax = 90 • or 270 • . Consequently, with the decrease in the angle between the drilling direction of the horizontal well and the maximum horizontal principal stress, the FIP of the intact rock matrix decreased.

Fracture Initiation of the Bedding Planes
The existence of bedding planes suggests that there are two different FIMs (failures of the intact rock matrix and bedding planes). Further, failure of the bedding planes implies shear failure along the bedding planes and tensile failure across the bedding planes. In this section, the FIP and FIL of the bedding planes have been calculated, as discussed in Section 2 with Mohr-Coulomb criterion for shear failure along the bedding planes and tensile failure criterion for tensile failure across the bedding planes.
Energies 2020, 13, 286 10 of 26 drilling direction of the horizontal well and the maximum horizontal principal stress, the FIP of the intact rock matrix decreased.

Fracture Initiation of the Bedding Planes
The existence of bedding planes suggests that there are two different FIMs (failures of the intact rock matrix and bedding planes). Further, failure of the bedding planes implies shear failure along the bedding planes and tensile failure across the bedding planes. In this section, the FIP and FIL of the bedding planes have been calculated, as discussed in Section 2 with Mohr-Coulomb criterion for shear failure along the bedding planes and tensile failure criterion for tensile failure across the bedding planes.
The change in the FIP at every point of the wellbore was drawn to simplify the fracture initiation analysis. Figure  Besides, the lowest point of each curve in Figure 7 represents the FIP and FIL of the bedding planes. (1) Based on analysis of the curve, Figure 7 reveals that the FIPs on the wellbore apparently exhibit periodic characteristics. The FIP at the point is equal to that at the point + 180°.
Consequently, there are two lowest points (( , ) and ( + 180°, )) for each curve, which means that there are two FILs on the wellbore with the FIP of . As a result, this study only analysed the FILs in the range −90-90° on the wellbore.
(2) In order to simplify the range of the influencing parameters, specific bedding planes were taken into consideration in Figure 7, which showed that the FIP curves on the wellbore of = The change in the FIP at every point of the wellbore was drawn to simplify the fracture initiation analysis. Figure 7 displays the FIPs on the wellbore (θ = 0 • − 360 • ) of the bedding planes with α H = 20 • , α B = 140 • or 320 • and β B = 10 • . In Figure 7a, the point A (50, 8.880) corresponds to that the FIP at the point θ = 50 • on the wellbore is equal to 8.880 MPa with α H = 20 • , α B = 140 • and β B = 10 • . Besides, the lowest point of each curve in Figure 7 represents the FIP and FIL of the bedding planes.
Energies 2020, 13, 286 10 of 26 drilling direction of the horizontal well and the maximum horizontal principal stress, the FIP of the intact rock matrix decreased.

Fracture Initiation of the Bedding Planes
The existence of bedding planes suggests that there are two different FIMs (failures of the intact rock matrix and bedding planes). Further, failure of the bedding planes implies shear failure along the bedding planes and tensile failure across the bedding planes. In this section, the FIP and FIL of the bedding planes have been calculated, as discussed in Section 2 with Mohr-Coulomb criterion for shear failure along the bedding planes and tensile failure criterion for tensile failure across the bedding planes.
The change in the FIP at every point of the wellbore was drawn to simplify the fracture initiation analysis. Figure  Besides, the lowest point of each curve in Figure 7 represents the FIP and FIL of the bedding planes. (1) Based on analysis of the curve, Figure 7 reveals that the FIPs on the wellbore apparently exhibit periodic characteristics. The FIP at the point is equal to that at the point + 180°.
Consequently, there are two lowest points (( , ) and ( + 180°, )) for each curve, which means that there are two FILs on the wellbore with the FIP of . As a result, this study only analysed the FILs in the range −90-90° on the wellbore.
(2) In order to simplify the range of the influencing parameters, specific bedding planes were taken into consideration in Figure 7, which showed that the FIP curves on the wellbore of = (1) Based on analysis of the curve, Figure 7 reveals that the FIPs on the wellbore apparently exhibit periodic characteristics. The FIP at the point θ x is equal to that at the point θ x + 180 • . Consequently, there are two lowest points ((θ c , P Bm f ) and (θ c + 180 • , P Bm f )) for each curve, which means that there are two FILs on the wellbore with the FIP of P Bm f . As a result, this study only analysed the FILs in the range −90-90 • on the wellbore.
(2) In order to simplify the range of the influencing parameters, specific bedding planes were taken into consideration in Figure 7, which showed that the FIP curves on the wellbore of α B = x and α B = 180 • + x were symmetrical and axisymmetric about the axis θ = 180 • . The planes share the same FIP and the sum of FILs is equal to 360 • . Irrespective of the value of β B , the FIP curves on the wellbore  First, the special cases were examined, in which the dip angle of the bedding planes was 0 • or 90 • . Figure 8a showed that the change in α B made no difference to the change in the FIP with α H when the dip angle was equal to 0 • (β B = 90 • ). This was because the bedding planes were parallel to the horizontal plane of the wellbore for β B = 90 • , which meant that the spatial relationship between the locations of the bedding planes and horizontal wellbore was not affected by α B ; therefore, no effect of the change in α B was observed on the stress state of the bedding planes. For tensile failure across the bedding planes, the FIP reached the maximum value (P Bm f max ) for α Hmax = 90 • or 270 • and the minimum value (P Bm f min ) for α Hmin = 0 • or 180 • . For shear slippage along the bedding planes, the FIP reached P Bm f max for α Hmin = 0 • or 180 • and P Bm f min for α Hmax = 90 • or 270 • . Besides, it was shear slippage which occurred, as its FIP was smaller. Figure 8g revealed that failure of the bedding planes did not occur when α H was equal to α B or α B ± 180 • with β B = 0 • . This was because the normal direction of the bedding planes was parallel to the axis of the horizontal wellbore, which meant that the normal stress projected onto the surface of the bedding planes was unable to induce failure of the bedding planes. In Wang's numerical simulation research, this phenomenon can be seen that when the wellbore passes through the weak surface vertically, the bedding planes has very little influence on the FIP [43]. Consequently, failure of the intact rock matrix occurred at the intersection of the bedding planes and horizontal wellbore, which resulted in transverse fractures in this case.
From On the other hand, when α B was equal to 90 • , α Hmax for shear slippage was equal to 0 • or 180 • with small β B (such as β B = 75 • ). For tensile failure, when |α B − 90 • | was equal to 90 • or 0 • , α Hmax and α Hmin remained constant, irrespective of the value of β B . When |α B − 90 • | decreased from 90 • to 0 • , α Hmin changed (decreased for α B < 90 • and increased for α B > 90 • ) less than 30 • and α Hmax changed (decreased for α B < 90 • and increased for α B > 90 • ) less than 20 • with the decrease in β B . Moreover, the changes in α Hmax and α Hmin were concentrated in the range of β B = 75 • − 60 • . In other words, the trend of the variation in the FIP with α H for tensile failure was mainly affected by α B . The change in β B had little influence on the trend, and the smaller the β B , the lower would be its influence. For shear slippage, when |α B − 90 • | was equal to 90 • , α Hmax and α Hmin remained constant, irrespective of the value of β B . When |α B − 90 • | ranged from 90 • to 0 • , α Hmin changed (increased for α B < 90 • and decreased for α B > 90 • ) less than 10 • and α Hmax changed (increased for α B < 90 • and decreased for α B > 90 • ) less than 30 • with the decrease in β B . The larger the |α B − 90 • | and the smaller the β B , the lower would be the change in α Hmax and α Hmin . However, when α B was equal to 90 • , turned from 0 • or 180 • to 90 • or 270 • (sudden change), whereas α Hmin changed from 90 • or 270 • to 0 • or 180 • (gradual change), with the decrease in β B . Consequently, the trend of the variation in the FIP with α H for shear slippage was also mainly affected by α B and the smaller the |α B − 90 • |, the lower would be its influence. Moreover, the trend was affected by β B ; the larger the β B , the greater would be its influence. Particularly, when |α B − 90 • | was small enough and β B large, the trend was more controlled by β B , which led to the sudden change in α Hmax and gradual change in α Hmin . It can be concluded that for fracture initiation analysis of the bedding planes, the drilling direction (α H ) corresponding to the minimum FIP is mainly affected by the dip direction of the bedding planes, while the dip angle has little influence. On the other hand, when the normal direction of the bedding planes is nearly parallel to the maximum horizontal principal stress and the dip angle is quite small, the drilling direction corresponding to the minimum FIP is more influenced by the dip angle. and = 180°+ were symmetrical and axisymmetric about the axis = 180°. The planes share the same FIP and the sum of FILs is equal to 360°. Irrespective of the value of , the FIP curves on the wellbore of = and = are the same as those of = 180°− and = 360°− .
Consequently, for fracture initiation analysis of the bedding planes, it was assumed that ranged from 0° to 180°, ranged from 0° to 90° and ranged from 0° to 360°.
3.2.1. FIP of the Bedding Planes , to analyse the changes in the FIP of the bedding planes with drilling direction ( ) for different bedding planes. In Figure 8, seven images for different dip angles have been presented in different colors, which represent different , and in different line styles representing different FIMs. Bedding planes not only change the trend of the variation in the FIP with α H , but also influence the FIP. Consequently, α B and β B mainly contribute to the differences in P Bm f max and P Bm f min , while α H − α B and β B are the main factors influencing the FIP. The previous research, which found the angle of a horizontal well with the principal stresses had a higher impact on the FIP than the rock anisotropy [30], and the experimental research about the influence of bedding plane angle on the FIP [44] also proved this observation. In Figure 8, the difference between P Bm f max and P Bm f min (∆P Bm f ) changed much more with β B than with α B . When β B changed from 75 • to 15 • , the FIP of shear slippage ranged from 5-9 MPa to 5-14 MPa, while that of tensile failure ranged from 9-12 MPa to 3-70 MPa. It was obvious that β B influenced the FIP of the bedding planes more than α B and that tensile failure across the bedding planes was more sensitive. For tensile failure, when |α B − 90 • | ranged from 90 • to 0 • , ∆P Bm f increased with the increase in P Bm f max and P Bm f min . When β B decreased, ∆P Bm f increased with the increase in P Bm f max and decrease in P Bm f min . For shear slippage, P Bm f max was affected by β B more than by α B , and gradually increased with the decrease in β B . With the decrease in |α B − 90 • |, P Bm f max gradually decreased when β B was quite large, which decreased first and then increased when β B was quite small. Meanwhile, P Bm f min was affected more by α B than by β B , and gradually increased with the decrease in |α B − 90 • |. With the decrease in β B , P Bm f min gradually decreased when |α B − 90 • | was quite large, which decreased first and then increased when |α B − 90 • | was quite small.
The influence of the bedding planes on both the trend of the variation in the FIP with α H and its minimum and maximum values led to changes in the FIM with α H . For Nong'an oil shale, the change in the FIM with α H of the bedding planes is shown in Figure 8. When β B was large enough (from Figure 8a to Figure 8d), shear failure along the bedding planes occurred, irrespective of the values of α B and α H . When β B became smaller (from Figure 8e to Figure 8g), there were cases with α H = α HT (α HT was around α B ± 90 • ), in which tensile failure occurred with the increase in the injected hydraulic pressure. Moreover, with the decrease in |α B − 90 • |, the range of α HT decreased. However, the existence of an inflection point (α BT ), which could be seen in Figure 8e, suggested that when |α B − 90 • | was less than |α BT − 90 • |, there were no cases with α H = α HT . Meanwhile, α BT increased with the decrease in β B . It can be concluded that β B has a greater influence than α B on the FIM of the bedding planes and should therefore be the priority factor in FIM analysis. The influence of α B reveals that the closer the normal direction of the bedding planes and the direction of minimum horizontal principal stress, the smaller would be the range of α HT , with the precondition that β B is small enough.

Fracture Initiation Location of the Bedding Planes
When failure of the intact rock matrix occurred, the FIL was at the point θ = 0 • or 180 • on the wellbore, as discussed in Section 3.1. This section analyses the change in the FIL when failure of the bedding planes occurs. Figure 9 shows the change in the FIL (from −90 • to 90 • ) with α H for different bedding planes.
planes is nearly parallel to the maximum horizontal principal stress and the dip angle is quite small, the drilling direction corresponding to the minimum FIP is more influenced by the dip angle.
Bedding planes not only change the trend of the variation in the FIP with , but also influence the FIP. Consequently, and mainly contribute to the differences in and , while − and are the main factors influencing the FIP. The previous research, which found the angle of a horizontal well with the principal stresses had a higher impact on the FIP than the rock anisotropy [30], and the experimental research about the influence of bedding plane angle on the FIP [44] also proved this observation. In Figure 8, the difference between and (∆ ) changed much more with than with . When changed from 75° to 15°, the FIP of shear slippage ranged from 5-9 MPa to 5-14 MPa, while that of tensile failure ranged from 9-12 MPa to 3-70 MPa. It was obvious that influenced the FIP of the bedding planes more than and that tensile failure across the bedding planes was more sensitive. For tensile failure, when | − 90°| ranged from 90° to 0°, ∆ increased with the increase in and . When decreased, ∆ increased with the increase in and decrease in . For shear slippage, was affected by more than by , and gradually increased with the decrease in . With the decrease in | − 90°|, gradually decreased when was quite large, which decreased first and then increased when was quite small. Meanwhile, was affected more by than by , and gradually increased with the decrease in | − 90°|. With the decrease in , gradually decreased when | − 90°| was quite large, which decreased first and then increased when | − 90°| was quite small.
The influence of the bedding planes on both the trend of the variation in the FIP with and its minimum and maximum values led to changes in the FIM with . For Nong'an oil shale, the change in the FIM with of the bedding planes is shown in Figure 8. When was large enough (from Figure 8a  ( was around ± 90°), in which tensile failure occurred with the increase in the injected hydraulic pressure. Moreover, with the decrease in | − 90°|, the range of decreased. However, the existence of an inflection point ( ), which could be seen in Figure 8e, suggested that when | − 90°| was less than | − 90°|, there were no cases with = . Meanwhile, increased with the decrease in . It can be concluded that has a greater influence than on the FIM of the bedding planes and should therefore be the priority factor in FIM analysis. The influence of reveals that the closer the normal direction of the bedding planes and the direction of minimum horizontal principal stress, the smaller would be the range of , with the precondition that is small enough.

Fracture Initiation Location of the Bedding Planes
When failure of the intact rock matrix occurred, the FIL was at the point = 0° or 180° on the wellbore, as discussed in Section 3.1. This section analyses the change in the FIL when failure of the bedding planes occurs. Figure 9 shows the change in the FIL (from −90° to 90°) with for different bedding planes.  the increase in and decrease in | − 90°|, the range of the FIL gradually shrank with = 90° or −90° as the centre. Further, there were only two initiation points, which in general deviated from = 0° or 180° in the cases with ≠ 0°. For shear slippage, when was equal to 90°, made no difference to the FIL (| | = 33°− 43°), and four initiation points was on the wellbore, as shown in Figure 9a. When was equal to 0°, the change in the FIL was so great that | | was equal to 0°, while | | was equal to 90°, except for the case with = 0°. Further, when | | was not  Further, there were only two initiation points, which in general deviated from θ = 0 • or 180 • in the cases with β B 0 • . For shear slippage, when β B was equal to 90 • , α B made no difference to the FIL (|θ T | = 33 • − 43 • ), and four initiation points was on the wellbore, as shown in Figure 9a. When β B was equal to 0 • , the change in the FIL was so great that |θ Tmin | was equal to 0 • , while |θ Tmax | was equal to 90 • , except for the case with α B = 0 • . Further, when |θ T | was not equal to 0 • or 90 • , there were four initiation points. When β B was not equal to 0 • or 90 • , with the decrease in β B , |θ Tmax | and |θ Tmin | gradually decreased until θ Tmin = 0 • .|θ T | reached the maximum of |θ Tmax | 90 • with four initiation points when α H was equal to α Hmax , as shown in Figure 9d. The α Hmax and α Hmin mentioned above referred to those in Figure 8 corresponding to P Bm f max and P Bm f min . Particularly, when α B was equal to 90 • , there exhibited θ T = ±|θ Tmax | ±90 • with α H = 90 • or 270 • and θ T = θ Tmin with α H = 0 • or 180 • , as shown in Figure 9c. In Section 3.2.1, it was revealed that when α B was equal to 90 • , α Hmax of shear slippage changed from 0 • or 180 • to 90 • or 270 • with the decrease in β B . Consequently, the influence of β B on the FIP was larger than that on the FIL, as the α H corresponding to θ Tmax and θ Tmin displayed no change with the decrease in β B when α B was equal to 90 • .
Bedding planes change the phenomenon that the FIL of hydraulic fracturing is at the point θ = 0 • or 180 • on the wellbore. This phenomenon could be observed in Hou's experimental researches, in which there is a deflection of the FIL away from the points θ = 0 • or 180 • [45]. Moreover, for shear slippage along the bedding planes, there is a proper drilling direction which changes the number of initiation points for the FIL from two to four. The trend of the variation in the FIL with α H is mainly influenced by α B . Further, the change in β B also contributes to the difference in the trend of the variation in the FIL with α H , though the influence is smaller than that on the trend of the variation in the FIP with α H . The larger the β B , the greater is the deviation in the FIL from θ = 0 • or 180 • . Generally, the dip angle is the main factor influencing the deviation in the FIL from θ = 0 • or 180 • .

Influence of the Bedding Planes on Fracture Initiation
Hydraulic fracture initiation in oil shale has two different modes: failure of the intact rock matrix and failure of the bedding planes, and three different failure criteria are used to evaluate the FIP of the intact rock matrix, tensile failure across the bedding planes and shear slippage along the bedding planes. When the hydraulic pressure reaches the minimum FIP of the three FIMs, hydraulic fracture initiates. Consequently, through fracture initiation analysis of the intact rock matrix and bedding planes in Section 3.2, the influences of α H (0 • − 360 • ), α B (0 • − 180 • ) and β B (0 • − 90 • ) on FIP were studied, as reported in the following sections. Figure 10 was plotted to confirm the relative priorities of the three FIMs for different bedding planes of Nong'an oil shale. equal to 0° or 90°, there were four initiation points. When was not equal to 0° or 90°, with the decrease in , | | and | | gradually decreased until = 0°. | | reached the maximum of | | ≠ 90° with four initiation points when was equal to , as shown in Figure 9d. The and mentioned above referred to those in Figure 8 corresponding to and . Particularly, when was equal to 90°, there exhibited = ±| | ≠ ±90° with = 90° or 270° and = with = 0° or 180°, as shown in Figure 9c. In Section 3.2.1, it was revealed that when was equal to 90°, of shear slippage changed from 0° or 180° to 90° or 270° with the decrease in . Consequently, the influence of on the FIP was larger than that on the FIL, as the corresponding to and displayed no change with the decrease in when was equal to 90°. Bedding planes change the phenomenon that the FIL of hydraulic fracturing is at the point = 0° or 180° on the wellbore. This phenomenon could be observed in Hou's experimental researches, in which there is a deflection of the FIL away from the points = 0° or 180° [45]. Moreover, for shear slippage along the bedding planes, there is a proper drilling direction which changes the number of initiation points for the FIL from two to four. The trend of the variation in the FIL with is mainly influenced by . Further, the change in also contributes to the difference in the trend of the variation in the FIL with , though the influence is smaller than that on the trend of the variation in the FIP with . The larger the , the greater is the deviation in the FIL from = 0° or 180°. Generally, the dip angle is the main factor influencing the deviation in the FIL from = 0° or 180°.

Influence of the Bedding Planes on Fracture Initiation
Hydraulic fracture initiation in oil shale has two different modes: failure of the intact rock matrix and failure of the bedding planes, and three different failure criteria are used to evaluate the FIP of the intact rock matrix, tensile failure across the bedding planes and shear slippage along the bedding planes. When the hydraulic pressure reaches the minimum FIP of the three FIMs, hydraulic fracture initiates. Consequently, through fracture initiation analysis of the intact rock matrix and bedding planes in Section 3.2, the influences of (0°− 360°), (0°− 180°) and (0°− 90°) on FIP were studied, as reported in the following sections. Figure 10 was plotted to confirm the relative priorities of the three FIMs for different bedding planes of Nong'an oil shale.
From previous discussion, it was obvious that the FIP of the intact rock matrix did not exhibit any relationship with the bedding planes and that the FIP of the bedding planes would not be affected by for = 90°. Figure 10 displays the FIP of the intact rock matrix in the form of black solid lines, and the different line styles (solid and dashed lines) with different colors (red, green, blue, orange, magenta and dark yellow) represent shear slippage along the bedding planes and tensile failure across the bedding planes for different . In Figure 10g, the excessively large calculation results of the FIPs were ignored for better observation. At the moment the increasing hydraulic pressure reaches the minimum FIP, hydraulic fracture initiates. As shown in Figure 10a, fracture initiation of the bedding planes occurs for = , while that of the intact rock matrix occurs for = . Moreover, the possibility of shear slippage along the bedding planes was higher for Nong'an oil shale when failure of the bedding planes occurred. With the decrease in , tensile failure across the bedding planes gradually occurred, as observed in Figure 10 in the form of dashed lines. Moreover, when | − 90°| decreased, the range of (for which tensile failure across the bedding planes occurred) decreased as well. (1) Four different FIP relationships for different FIMs with different bedding planes and drilling directions From previous discussion, it was obvious that the FIP of the intact rock matrix did not exhibit any relationship with the bedding planes and that the FIP of the bedding planes would not be affected by α B for β B = 90 • . Figure 10 displays the FIP of the intact rock matrix in the form of black solid lines, and the different line styles (solid and dashed lines) with different colors (red, green, blue, orange, magenta and dark yellow) represent shear slippage along the bedding planes and tensile failure across the bedding planes for different α B . In Figure 10g, the excessively large calculation results of the FIPs were ignored for better observation.
At the moment the increasing hydraulic pressure reaches the minimum FIP, hydraulic fracture initiates. As shown in Figure 10a, fracture initiation of the bedding planes occurs for α H = α HC , while that of the intact rock matrix occurs for α H = α HAC . Moreover, the possibility of shear slippage along the bedding planes was higher for Nong'an oil shale when failure of the bedding planes occurred. With the decrease in β B , tensile failure across the bedding planes gradually occurred, as observed in Figure 10 in the form of dashed lines. Moreover, when |α B − 90 • | decreased, the range of α H (for which tensile failure across the bedding planes occurred) decreased as well. For better analysis of the change in the FIM, the FIP relationship between the intact rock matrix and the bedding planes was divided into four types, based on the variations in their maximum and minimum FIPs (P Rm f max , P Rm f min , P Bm f max and P Bm f min ) with α H : (I): P Rm f max > P Bm f max and P Rm f min < P Bm f min ; (II): P Rm f max > P Bm f max and P Rm f min > P Bm f min ; (III): P Rm f max < P Bm f max and P Rm f min > P Bm f min ; (IV): P Rm f max < P Bm f max and P Rm f min < P Bm f min .
(1) Four different FIP relationships for different FIMs with different bedding planes and drilling directions In the cases with any α B in Figure 10a,b, |α B − 90 • | = 0 • or 30 • in Figure 10c,d, the FIP relationship belonged to type I. The distribution of α HC was centred on α H = 0 • or 180 • only for α B = 0 • or 90 • . With the decrease in |α B − 90 • | (α B 0 • ), the centre point of the α HC distribution was closer to α H = 0 • or 180 • and the distribution range would be smaller. Besides, with the decrease in β B , the distribution range became larger. In particular, the change in the FIP with α H for different FIMs showed the same trend in Figure 10a, leading to the result that the distribution of α HC was centred on α H = 0 • or 180 • .
In the cases with |α B − 90 • | = 60 • or 90 • in Figure 10c and |α B − 90 • | = 60 • in Figure 10d, the FIP relationship belonged to type II. α HC revealed two different distributions. One was that α HC varied from 0 • to 360 • , which meant that failure of the bedding planes occurred in any drilling direction. The other was that α HAC was distributed on both sides of α Hmax (corresponding to P Bm f max ) of the bedding planes. This was because with the decrease in |α B − 90 • |, α Hmax of the bedding planes showed the tendency to move outward from α H = 0 • or 180 • , which led to the possibility that around α Hmax of the bedding planes, the FIP of the intact rock matrix was lower. Consequently, the smaller the values of |α B − 90 • | and β B , the more likely would be the occurrence of the second distribution with a wider distribution of α HAC .
In the cases with |α B − 90 • | = 90 • in Figure 10d, |α B − 90 • | = 60 • or 90 • in Figure 10e-g, the FIP relationship belonged to type III. The distribution of α HAC exhibited a similar trend as the second distribution of type II. In other words, α HAC was distributed on both sides of α Hmax of the bedding planes. The distribution range of α HAC increased with the decrease in |α B − 90 • | and β B .
In the cases with |α B − 90 • | = 0 • or 30 • in Figure 10e-g, the FIP relationship belonged to type IV. The distribution of α HC was centred on α H = 0 • or 180 • only for α B = 90 • . With the decrease in |α B − 90 • | (α B 0 • ), the centre point of the α HC distribution was closer to α H = 0 • or 180 • , and the distribution range was smaller. Moreover, the distribution range of α HC increased with the decrease in |α B − 90 • | and increase in β B .
(2) Influence of bedding planes on the FIP, FIL and FIM of oil shale The existence of α HC proves that bedding planes greatly change the FIP of oil shale. Without the bedding planes, the FIP of the intact rock matrix reached the maximum for α H = 0 • or 180 • and the minimum for α H = 90 • or 270 • . Besides, fracture initiated at the point θ = 0 • or 180 • on the wellbore. However, when the influence of bedding planes is considered, hydraulic fractures can initiate through two different modes (failure of the intact rock matrix and failure of the bedding planes). The change in the FIP with drilling direction becomes different from the rock without regard to the bedding planes, as does the change in the FIL. When α H was equal to α HC , failure of the bedding planes occurred, with lower FIPs and different FILs, which did not correspond to θ = 0 • or 180 • , as shown in Figures 9  and 10. The distribution range and interval length of α HC were controlled by the relationship between P Rm f max , P Rm f min , P Bm f max , and P Bm f min . Consequently, bedding planes not only reduce the FIP, but also change the FIL.
The influence of bedding planes on FIP is mainly determined by the effect of α B and β B on the trend of the variation in the FIP with the drilling direction, which leads to the difference in the relationships between P Rm f max , P Rm f min , P Bm f max and P Bm f min . Consequently, the distributions of α HC and α HAC change for different bedding planes (different α B and β B ). From an analysis of the FIP curves, it was found that irrespective of the value of β B , P Bm f min was larger than P Rm f min for small |α B − 90 • |. However, when |α B − 90 • | was large enough, P Bm f min was larger than P Rm f min first and then lesser than P Rm f min with the decrease in β B . Consequently, there existed an inflection point α BC . It was P Bm f min ≥ P Rm f min for |α B − 90 • | ≤ |α BC − 90 • |, which changed from P Bm f min > P Rm f min to P Bm f min < P Rm f min for |α B − 90 • | > |α BC − 90 • | as P Bm f min decreased with the decrease in β B .
In the cases with |α B − 90 • | ≤ |α BC − 90 • |, the FIP relationship was type I and type IV. The transformation between two FIP relationships was mainly attributed to the influence of β B on P Bm f max .
As P Bm f max increased with the decrease in β B , the relationship changed from P Bm f max < P Rm f max to P Bm f max > P Rm f max . Consequently, there was an inflection point β BC . It was P Bm f max < P Rm f max for β B > β BC and P Bm f max > P Rm f max for β B < β BC . With the decrease in |α B − 90 • |, the inflection point β BC decreased. The distribution of α HC in the cases of type I and type IV was around α H = 0 • or 180 • . Therefore, the relationship between P Bm f max and P Rm f max had a small effect on the distribution of α HC , and mainly affected the interval length.
In the cases with |α B − 90 • | ≤ |α BC − 90 • |, the FIP relationship was type I, type II and type III. The transformation from type I to type II was mainly attributed to the change in P Bm f min , which decreased with the decrease in β B and eventually became lower than P Rm f min . In the cases of type I, α HC was mainly distributed around α H = 0 • or 180 • , with α B affecting the central point of its distribution and β B changing its interval length. In the cases of type II, α HAC mainly distributed around α Hmax of the bedding planes, with α B affecting the central point of its distribution and β B changing its interval length. The distributions of α HC and α HAC almost interchanged in this transformation. The transformation from type II to type III was mainly attributed to the change in P Bm f max , which increased with the decrease in β B , varying from P Bm f max < P Rm f max to P Bm f max > P Rm f max . There also existed an inflection point β BC . It was P Bm f max < P Rm f max for β B > β BC and P Bm f max > P Rm f max for β B < β BC . With the decrease in |α B − 90 • |, the inflection point β BC decreased. In the cases of type II and type III, α HAC was mainly distributed around α Hmax of the bedding planes. For cases with some specific bedding planes, failure of the bedding planes occurred in any drilling direction, such as the case with |α B − 90 • | = 90 • or 60 • in Figure 10c. In this transformation, mainly a change in the interval length was observed. Consequently, the two transformations revealed that the relationship between P Bm f min and P Rm f min was the main reason for the change in the α HC distribution, which was controlled by α B and β B . According to FIP analysis of the bedding planes in the earlier section, it was found that P Bm f min was more affected by α B , while P Bm f max was more affected by β B . Consequently, α B was the main factor influencing the distribution of α HC . In the analysis on oil shale hydraulic fracturing, the dip direction is the first parameter to confirm the main region of α HC and α HAC , and the dip angle mainly controls the interval length.

Influence of the Bedding Planes on the Design of Horizontal Wellbore Trajectory
Bedding planes changed the trend of the variation in the FIP with drilling direction, as well as the maximum and minimum FIPs and the corresponding FIMs. Consequently, for hydraulic fracturing in oil shale with different bedding planes, the drilling direction and FIM which corresponded to the minimum FIP were different, and the horizontal wellbore trajectory design also needed to be changed. In this section, it was found that the influence of bedding planes on the distribution of α HC and α HAC , as well as on the maximum and minimum FIPs, was mainly attributed to the influence on the relationship between P Rm f max , P Rm f min , P Bm f max and P Bm f min , which led to different optimum drilling directions. Consequently, it was necessary to determine the FIP relationship between different FIMs in order to analyse the optimum drilling direction.
For the bedding planes which ensured that the FIP relationship belonged to type I or type IV, the FIP of oil shale reached the minimum value (P m f min ) for α H = 90 • or 270 • , which referred to P Rm f min of the intact rock matrix. When the drilling direction of the horizontal well corresponded to α H = 90 • or 270 • , the hydraulic pressure required was the easiest to achieve. In these cases, fractures tended to initiate across the intact rock matrix and extend with the increase in the hydraulic pressure; the preferred initial fracture morphology was transverse fractures.
For the bedding planes which ensured that the FIP relationship belonged to type II or type III, the distribution range of α HC was much more than that of α HAC , which meant that the drilling directions for the failure of the bedding planes were much more than those for the failure of the intact rock matrix. Moreover, there were cases with specific bedding planes, in which the failure of the bedding planes occurred in any drilling direction. The decrease in β B induced an increase in P Bm f max , which gradually became larger than P Rm f max , thus leading to the wider distribution of α HAC . The FIP of oil shale reached the minimum value P m f min , which referred to P Bm f min of the bedding planes, for α H = α Hmin (the drilling direction corresponding to P Bm f min ). Moreover, with the decrease in β B , failure of the bedding planes gradually changed from shear slippage along the bedding planes to tensile failure across the bedding planes, and the larger the |α B − 90 • | the more probable would be tensile failure. Consequently, when the drilling direction of the horizontal well corresponded to α H = α Hmin (the drilling direction corresponding to P Bm f min ), the hydraulic pressure required was the easiest to achieve. In these cases, fractures tended to initiate along the bedding planes and extend with the increase in the hydraulic pressure; the preferred initial fracture morphology was bedding fractures.
It is well known that hydraulic fractures prefer to activate the weak planes in the propagation of hydraulic fractures, due to their much lower tensile and shear strength [46][47][48]. This could lead to a complex fracture network. However, in the stage of hydraulic fracture initiation, the activation of the bedding planes is unfavored to the complex fracture network and better simulated reservoir volume [24,34]. For the bedding planes which ensure that the FIP relationship belongs to type II or type III, the drilling direction, which corresponds to the minimum FIP (P m f min ), is more conducive to the formation of bedding fractures than transverse fractures. Consequently, the horizontal well drilling direction needs to be changed for the preferred fracture morphology, which leads to the different optimal directions. If only hydraulic pressure is required, for the bedding planes which ensure that the FIP relationship belongs to type I or type IV, the optimal drilling direction should correspond to α H = 90 • or 270 • . For the bedding planes which ensure that the FIP relationship belongs to type II or type III, the optimal drilling direction should correspond to α H = α Hmin (the drilling direction corresponding to P Bm f min ). However, if the transverse fractures are also required, the optimal drilling direction should be kept the same as that for the former cases but changed to that corresponding to the minimum FIP of the intact rock matrix for the latter cases. In the latter cases, the drilling direction was the boundary of the α HC distribution range.
It is worth noting that, although the research above targets Nong'an oil shale, the developed fracture initiation model has broader applicability for the other sedimentary rocks, which shares the same assumptions stated in the first paragraph of the Section 2. As the hydraulic fracture initiates when the injected water pressure reaches the minimum FIP, by calculating the FIPs of different FIMs, the FIP, FIL, FIM, and the initial fracture morphology can be estimated. If we can obtain the accurate physical parameters in Table 1 of the target formation, the hydraulic fracture initiation laws in the formations, as well as the influence of bedding planes, can be analysed as well.

Experimental Validation
To validate the developed fracture initiation model, experiments on Nong'an oil shale were carried out. In this section, the experimental equipment, the preparation of samples and the test procedure are explained briefly.

Experimental Setup and Sample Preparation
The experimental equipment in Figure 11 is used to simulate the hydraulic fracturing in reservoir, consisting of four parts (the true tri-axial assembly, hydraulic pressure unit, high-pressure water pump, and data processing and control system). Cubic blocks of 200 mm are positioned in between the pressurized pistons, which simulate the in-situ stresses in three different directions. These pistons are controlled by the hydraulic cylinders in the tri-axial assembly and powered by the hydraulic pressure unit. In addition, the high-pressure water pump is used to inject the hydraulic fluid into the sealed The 200 mm cube samples were prepared using the pre-prepared cement and Nong'an oil shale of 100 × 200 × 200 mm, as shown in Figure 12. A 120 mm hole of φ 12 mm was drilled at the centre of the samples, and the upper section of the hole (80 mm) was enlarged to φ 16 mm. High-strength adhesive glue was used to seal the simulated wellbore of 14 mm diameter. In addition, in Figure  12 corresponded to that in the fracture initiation model and was set as the variable of 60°, 30°, and 0° in the experiments. The dip direction of bedding planes was set the same of = 0°, and the physical and mechanical properties of the oil shale samples were presented in Table 1.

Test Procedure
Experiments were conducted with the in-situ stresses of = 6 MPa, = 8 MPa, = 10 MPa. Additionally, the simulated wellbore was set in the direction of the maximum horizontal principal stress, which meant that the drilling direction of the horizontal wellbore was along the maximum horizontal principal stress ( = 90°). The fracturing fluid with a viscosity of 5 mPa/s was pumped into the sealed wellbore at an injection rate of 10 mL/min until a fracture initiated and propagated. The hydraulic pressure was recorded, and the sample was split open along the fracture plane to observe the fracture morphology after the experiments. The 200 mm cube samples were prepared using the pre-prepared cement and Nong'an oil shale of 100 × 200 × 200 mm, as shown in Figure 12. A 120 mm hole of ϕ 12 mm was drilled at the centre of the samples, and the upper section of the hole (80 mm) was enlarged to ϕ 16 mm. High-strength adhesive glue was used to seal the simulated wellbore of 14 mm diameter. In addition, β B in Figure 12 corresponded to that in the fracture initiation model and was set as the variable of 60 • , 30 • , and 0 • in the experiments. The dip direction of bedding planes was set the same of α B = 0 • , and the physical and mechanical properties of the oil shale samples were presented in Table 1. The 200 mm cube samples were prepared using the pre-prepared cement and Nong'an oil shale of 100 × 200 × 200 mm, as shown in Figure 12. A 120 mm hole of φ 12 mm was drilled at the centre of the samples, and the upper section of the hole (80 mm) was enlarged to φ 16 mm. High-strength adhesive glue was used to seal the simulated wellbore of 14 mm diameter. In addition, in Figure  12 corresponded to that in the fracture initiation model and was set as the variable of 60°, 30°, and 0° in the experiments. The dip direction of bedding planes was set the same of = 0°, and the physical and mechanical properties of the oil shale samples were presented in Table 1.

Test Procedure
Experiments were conducted with the in-situ stresses of = 6 MPa, = 8 MPa, = 10 MPa. Additionally, the simulated wellbore was set in the direction of the maximum horizontal principal stress, which meant that the drilling direction of the horizontal wellbore was along the maximum horizontal principal stress ( = 90°). The fracturing fluid with a viscosity of 5 mPa/s was pumped into the sealed wellbore at an injection rate of 10 mL/min until a fracture initiated and propagated. The hydraulic pressure was recorded, and the sample was split open along the fracture plane to observe the fracture morphology after the experiments.

Test Procedure
Experiments were conducted with the in-situ stresses of σ V = 6 MPa, σ h = 8 MPa, σ H = 10 MPa. Additionally, the simulated wellbore was set in the direction of the maximum horizontal principal stress, which meant that the drilling direction of the horizontal wellbore was along the maximum horizontal principal stress (α H = 90 • ). The fracturing fluid with a viscosity of 5 mPa/s was pumped into the sealed wellbore at an injection rate of 10 mL/min until a fracture initiated and propagated. The hydraulic pressure was recorded, and the sample was split open along the fracture plane to observe the fracture morphology after the experiments.

Experimental Results
In the calculation of fracture initiation pressure corresponding to the experiments, β B was set as 60 • , 30 • and 0 • corresponding to the experiments, with the conditions of σ V = 6 MPa, σ h = 8 MPa, σ H = 10 MPa, α H = 90 • , α B = 0 • , P p = 0 MPa. The theoretical results and the experimental results were shown in Figure 13.
Energies 2020, 13 Figure 13.  Figure 13 presents the theoretical and experimental results of the FIPs under the giving conditions. The error between the theoretical results and the experimental results ranges from 7% to 9%. Since the FIP will increases as the viscosity and injection rate increase when they are at relatively low values [11,48,49], the errors of 7% to 9% may attribute to the giving conditions in the experiments (viscosity of 5 mPa/s and injection rate of 10 mL/min). By comparing the theoretical and experimental results, it can be found that this developed fracture initiation model can predict the fracture initiation pressure in oil shale well. In addition, since the main fracture morphology observed around the sealed wellbore in these experiments was bedding fractures, the fracture morphology for = 60° was taken as an example in Figure 14. Figure 14 shows the main fracture plane along bedding planes and the secondary fracture induced by natural fracture. The theoretical results showed that shear slippage along bedding planes occurred under these experimental conditions, which corresponded to the main fracture morphology observed. Although the fracture initiation position in the experiments is hard to be confirmed by observing the fracture morphology, the estimation of the fracture initiation model in the fracture morphology can be validated to some extent.   Figure 13 presents the theoretical and experimental results of the FIPs under the giving conditions. The error between the theoretical results and the experimental results ranges from 7% to 9%. Since the FIP will increases as the viscosity and injection rate increase when they are at relatively low values [11,48,49], the errors of 7% to 9% may attribute to the giving conditions in the experiments (viscosity of 5 mPa/s and injection rate of 10 mL/min). By comparing the theoretical and experimental results, it can be found that this developed fracture initiation model can predict the fracture initiation pressure in oil shale well. In addition, since the main fracture morphology observed around the sealed wellbore in these experiments was bedding fractures, the fracture morphology for β B = 60 • was taken as an example in Figure 14. Figure 14 shows the main fracture plane along bedding planes and the secondary fracture induced by natural fracture. The theoretical results showed that shear slippage along bedding planes occurred under these experimental conditions, which corresponded to the main fracture morphology observed. Although the fracture initiation position in the experiments is hard to be confirmed by observing the fracture morphology, the estimation of the fracture initiation model in the fracture morphology can be validated to some extent.  Figure 13.  Figure 13 presents the theoretical and experimental results of the FIPs under the giving conditions. The error between the theoretical results and the experimental results ranges from 7% to 9%. Since the FIP will increases as the viscosity and injection rate increase when they are at relatively low values [11,48,49], the errors of 7% to 9% may attribute to the giving conditions in the experiments (viscosity of 5 mPa/s and injection rate of 10 mL/min). By comparing the theoretical and experimental results, it can be found that this developed fracture initiation model can predict the fracture initiation pressure in oil shale well. In addition, since the main fracture morphology observed around the sealed wellbore in these experiments was bedding fractures, the fracture morphology for = 60° was taken as an example in Figure 14. Figure 14 shows the main fracture plane along bedding planes and the secondary fracture induced by natural fracture. The theoretical results showed that shear slippage along bedding planes occurred under these experimental conditions, which corresponded to the main fracture morphology observed. Although the fracture initiation position in the experiments is hard to be confirmed by observing the fracture morphology, the estimation of the fracture initiation model in the fracture morphology can be validated to some extent.

Research Prospect
As discussed above, the horizontal wellbore trajectory design in oil shale need to consider the influence of bedding planes, which affects the FIP, FIL and FIM. However, nature fractures distributed in oil shale would make the fracture initiation more complex than what considers in this paper and would make great difference on the fracture propagation as well. In addition, in real world cases, the undulation of wellbore and the local variations in dip in oil shale are inevitable, which require additional optimization on the wellbore trajectory, as the intersection between the bedding planes and the wellbore will be changed. By simplifying the undulation of horizontal wellbore to the intersection change, the calculation results of the preferred wellbore trajectory range for several major bedding planes need to be integrated for fracture complexity and easier fracture initiation. The final wellbore trajectory will induce the multiple fracture initiation, as the injected water pressure is higher than the minimum FIP required for most well sections. Therefore, future work needs to consider the influence of stress shadow effects on the local stress distribution [50] and the influence of natural fractures on the fracture propagation [51,52].

Conclusions
This study set different failure criteria for different FIMs of oil shales: failure of the intact rock matrix and failure of the bedding planes (shear slippage along the bedding planes and tensile failure across the bedding planes) to analyse the variation in the FIP and FIL for different FIMs. With a case study on Nong'an oil shale, the changes in the FIP, FIL, FIM, and the wellbore trajectory design due to the influence of bedding planes were discussed. In addition, experiments for the validation of the fracture initiation model were conducted. Through this research, the following conclusions were drawn.
(1) A fracture initiation model, which considers the anisotropic strength induced by bedding planes, has been developed for horizontal wells in oil shale. The error between the theoretical and experimental results for the fracture initiation pressure in Nong'an oil shale ranges from 7% to 9%.
(2) When failure of the intact rock matrix occurred, the FIP reached the minimum value (P Rm f min ) with the angle between the drilling direction and the minimum horizontal principal stress of 90 • or 270 • (α H = 90 • or 270 • ). When failure of the bedding planes occurred, the trend of the FIP changing with drilling direction was mainly controlled by bedding planes' dip direction (represented by the angle between the normal direction of the bedding planes and the minimum horizontal principal stress, α B ), with little effect of dip angle (90 • − β B ). The smaller the β B , the lesser was its influence. Particularly, for small |α B − 90 • | and large β B , the trend of shear slippage was more affected by β B . Moreover, both α B and β B affected the maximum and minimum FIPs, and β B was the main influencing factor to which tensile failure was sensitive. Consequently, α B − α H and β B were the main factors influencing the FIP, and the optimal design of the horizontal wellbore trajectory in oil shale should take the change in the bedding planes into consideration.
(3) When failure of the intact rock matrix occurred, the FIL was at the point θ = 0 • or 180 • on the wellbore. On the other hand, when failure of the bedding planes occurred, it was away from the point θ = 0 • or 180 • . Moreover, the dip direction of the bedding planes was the main factor affecting the trend of the variation in the FIL changing with drilling direction, and the dip angle mainly controlled the degree of deviation. Particularly, when shear slippage along the bedding planes occurred, there existed α H which turned the two initiation points to four on the wellbore. Consequently, for the anisotropic rock-like oil shale, the FIL changes with the FIM, and is not fixed at θ = 0 • or 180 • .
(4) For Nong'an oil shale, the bedding planes made a difference to the changes in the FIP, FIL and FIM with drilling direction. Failure of the bedding planes occurred for α H = α HC , while failure of the intact rock matrix occurred for α H = α HAC . The optimal design of the horizontal wellbore trajectory revealed two different directions, depending on the purpose. If only hydraulic pressure was required, the optimal drilling direction corresponded to α H = 90 • or 270 • for the cases with the minimum FIP corresponding to failure of the intact rock matrix, whereas it was α H = α Hmin (the drilling direction corresponding to P Bm f min ) for the cases with the minimum FIP corresponding to failure of the bedding planes. However, if transverse fractures were also required, the optimal drilling direction was kept the same for the former cases and changed to that corresponding to the minimum FIP of the intact rock matrix, which was the boundary of the α HC distribution range, for the latter cases.
In this paper, we mainly completed the research about the change in the design of wellbore trajectory due to the influence of bedding planes on the fracture initiation and calculated the minimum required FIP for different bedding planes and wellbore trajectories. The further researches need to consider the influence of stress shadow effects and natural fractures on the local stress distribution and fracture propagation in oil shale.