Fracture Initiation Model of Shale Fracturing Based on Effective Stress Theory of Porous Media

Unconventional oil and gas are important resources of future energy supply, and shale gas is the focus of the development of unconventional resources. Shale is a special kind rock of porous medium, and an orderly structure of beddings aligned in the horizontal direction where causing the strong elastic anisotropy of shale is easy. A new model has been established to calculate the fracture initiation pressure with the consideration of mechanical characteristics of shale and the anisotropic tensile strength when judging rock failure. The fracture initiation model established in this paper accurately reflects the stress anisotropy and matches well with the actual situation in porous media. Through the sensitivity analysis, the results show that σv/σH , Ev/EH , υv/υH , m/s, and A/B have a certain impact on the tangential stress when the circumferential angle changes, and there is a positive relationship between the initiation pressure and the above sensitive factors except for A/B. The results can provide a valuable and effective guidance for the prediction of fracture initiation pressure and fracture propagation mechanism under special stratum conditions of shale.


Introduction
Unconventional oil and gas are important resources of future energy supply.And shale gas is the focus of the development of unconventional resources.Shale gas can only be profitable by hydraulic fracturing of horizontal wells with the current volume fracturing techniques.Fracture initiation of horizontal well hydraulic fracturing is a technical difficulty related to fracturing technology.
As for fracture initiation, a lot of studies have been conducted on this topic in recent years [1][2][3][4][5][6][7][8].There are several studies that demonstrated the initiation of transverse fractures by laboratory observations [9,10] and field observations [11].Fairhurst [12] revealed the possibility of a fracture initiating as transverse and in the longitudinal direction.However, it is shown [13] that the axial stress is not a good predictor of transverse fracture initiation, because it remains constant during hydraulic pressurization.Based on the minimum strain energy density criterion, Ayatollahi and Sedighiani [14] took the effect of T-stress into account, and they studied and analyzed the effect of T-stress on the critical mode stress intensity factor of brittle and quasi-brittle materials.However, for all cases in their experiments, the material was assumed to be homogeneous and isotropic with linear elastic behavior.By true triaxial hydraulic fracturing experiments, Peng et al. [15] observed the diversion and the following propagation of hydraulic fractures both horizontally and vertically after initiation from a directional wellbore with a skew angle.They mentioned that due to strong heterogeneity of coal rocks, primary and secondary cracks were developed, which would lead to discontinuity and anisotropy of mechanical characteristics in coal.However, they did not make a further study of this problem in view of this limitation.Fallahzadeh et al. [16] modeled various scenarios of vertical and horizontal wells and in situ stress regimes; in addition to experimental studies, analytical solutions were developed to simulate the mechanism of fracture initiation in perforated boreholes in tight formations.Furthermore, it was found that stress anisotropy influences the fracturing mechanism in a perforated borehole and affects the geometry of the initiated near-wellbore fracture.Considering the existence of bedding planes, cracks and other structural planes are the preconditions of the simulated reservoir volume of shale formations.In order to analyze the effect of bedding planes on the propagation of hydraulic fractures in shale formation, Heng et al. [17] carried out three-point bending tests of notched cylindrical specimens with different bedding orientations, based on the distribution of stress field around a crack tip of anisotropic materials.In addition to the experimental methods, theoretical research and numerical simulation methods have also been applied to analyze the fracture initiation mechanism.Lecampion [18] proposed an approximation of the mixed criteria which accounts for solving a single nonlinear equation.He held that a higher mean compressive far-field stress tends to make the failure more strength-driven, while a higher differential far-field compressive stress promotes an energy-driven tensile failure.But they are solely interested in tensile fracture initiation and did not investigate shear failures which may occur in the beginning depending on the stress condition.Tunsakul et al. [19] investigated the failure behavior of underground gas storage caverns under high pressure.The observation and analysis results revealed that the lateral earth pressure coefficient at rest has a strong influence on the position of this initiation point, while the depth of the cavern has an insignificant effect.However, the findings in this research are based on limited testing conditions with specific rock, and therefore a broader set of studies is needed to enhance the reliability.Xie and Min [20] realized that there are several distinctive features of enhanced geothermal system stimulation compared with common hydraulic treatments in the hydrocarbon reservoirs.Based on the hydroshearing concept, they established generic models to estimate the location of the shearing onset, the required injection pressure, and the overall shearing growth direction during enhanced geothermal system hydraulic stimulation.They adopt the Coulomb failure criterion to define the shear strength of a single rock joint but neglected the cohesion of fracture.Mohr-Coulomb and Hoek-Brown [21] applied the shear failure strength criteria to the borehole stability analysis associated with a homogeneous linear poroelastic model, which ignored the influence of the intermediate principal stress on rock failure.Based on borehole stress solutions, Zhang et al. [22] analyzed the influencing factors, such as Poisson's ratio and shear stress, but the solutions were all derived from the homogeneous linear poroelastic theory.Sun et al. [23] thought that shear failure, plastic deformation, and the coupling effect near the fracture tip play an important role in the fracture initiation and propagation, and a new coupled model is presented.But they assumed that the rock density and the fluid viscosity and density are uniformly distributed in the formation without the consideration of rock anisotropy.An orderly structure of shale beddings aligned in the horizontal direction caused the strong elastic anisotropy of shale.Based on the rock transverse isotropic constitutive relation and seepage and deformation coupling numerical method, Li et al. [24] established a three-dimensional finite element numerical model for a horizontal well with perforating completion.Considering the elastic mechanics anisotropy difference of shale, the sensitivity analysis of perforation parameters on the fracture initiation pressure and location of the horizontal well was proposed.With both the numerical simulations and experiments, Huang et al. [25] discussed the initiation mechanisms of natural fractures during hydraulic fracturing.They concluded that the lateral stress coefficient plays a critical role in determining the stress magnitude and orientation around the fracture tip and predicted that there is a high probability for developing inclined fractures in the coals with strong heterogeneity.But no detailed reasons were given.Gong et al. [26] analyzed variation rules of fracture in initiation pressure and fracture starting point of hydraulic fractures in a radial well; based on fluid-solid coupling effect and the maximum tensile-stress criterion, they simulate and study local stress accumulation situation caused by the vertical and radial section in the drilling process and the fracturing section in the hydraulic fracturing process through the finite element method.By using Abaqus finite element calculation software, Guo et al. [27] established the fracture initiation models for a 3D single-stage threecluster perforation and a single-cluster perforation (containing natural fracture).The results show that the initiation pressure of the open-hole perforation is far below that of the casing perforation.However, the calculation model assumes that the shale reservoir is homogeneously linear and is different from the accrual situation.
To calculate the stress distribution around the borehole and judge the failure of the wellbore rock, the previous research work about fracture initiation of horizontal well fracturing basically assumed that shale is an isotropic homogeneous rock.In fact, the shale rock has obvious anisotropic mechanical properties.The influence of the anisotropy should be considered when calculating the effective stress distribution around the borehole and confirming the tensile strength of the rock.

Borehole Stress Distribution of Transverse Isotropic Shale
Because shale gas is mainly developed by horizontal wells or extended reach wells and there is a certain angle between the direction of the borehole axis and the geographic coordinate system or the stress coordinate system, therefore, the coordinate system transformation should be performed before calculating the borehole stress distribution, and the process is transforming the far-field in situ stresses from a geographic coordinate system to a borehole coordinate system.Firstly, the far-field in situ stress distribution is transformed from the principal stress coordinate system to the geographic coordinate system as shown in 2 Geofluids Figure 1(a).The stress transformation relationship can be listed as follows: where σ p is the far-field stress tensor under the principal stress coordinate system, σ g is the far-field stress distribution tensor under the geographic coordinate system, α pg is the azimuth of the maximum horizontal principal stress, and β pg is the angle between the direction of σ v and Z g -axis.
Secondly, the far-field in situ stress distribution is transformed from the geographic coordinate system to the borehole coordinate system as shown in Figure 1(b).The stress transformation relationship is shown as follows: The calculation method is based on the theoretical works of Lekhnitskii [28], Amadei [29], Vahid and Ahmad [30], and Lu et al. [31], and the induced stress distribution around the wellbore can be calculated according to the analytical work of Amadei [29] as follows: where Re is the real component of an imaginary number, Φ' i is the derivative of the stress analytical functions, and z i , μ i , and λ i are the complex numbers and can be calculated using the following formulations.Under the borehole coordinate system, z i is where a is the radius of the borehole, θ is the angle measured counterclockwise from the axis X b direction In the above expressions, β ij can be calculated by And a ij can be defined as For transverse isotropic shale formation as shown in Figure 2, the rock has the properties as follows: And the functions l 3 μ i and l 2 μ i are as follows:   4 Geofluids In (3), Φ' i can be calculated by where p w is the fluid pressure in the borehole.And the parameters Δ and η i can be calculated using Therefore, the borehole stress distribution can be solved according to the above theoretical model.Because shale is a porous medium, combined with the theory of effective stress, the effective stress around a borehole can be determinated by total stress and pore pressure as follows: where p p is the pore pressure and α ij is Biot's coefficient tensor; it describes the influence of pore pressure on the elastic solid matrix for the anisotropic case.
According to the results of Shao [32] and Tan et al. [33], Biot's coefficient α ij in three directions can be expressed as where K s is the bulk modulus of the solid phase and M ij is the stiffness matrix.M ij can be calculated using Therefore, the effective stress in three directions around a borehole can be defined as In fact, researchers are more concerned with the tangential stress around a borehole when fracturing, so we can convert the stress components as follows:

Mechanical Model of Fracture Initiation
For compressive strength of layered shale, the variation of uniaxial compressive strength in different inclination angles of the bedding is established using the empirical method by Hoek and Brown [21].And the theoretical calculation method of compressive strength is revised and consummated by Rao et al. [34], Ramamurthy et al. [35], Saroglou and Tsiambaos [36], Ismael et al. [37], and so on.However, there are few studies on the anisotropic tensile strength of shale.The existing fracture initiation model also does not consider the impact of anisotropic tensile strength of shale.In fact, the tensile strength of shale has obvious anisotropic characteristics, and it has an obvious effect on the fracture initiation.

Analysis of Anisotropic Tensile Strength of Layered Shale.
The stress state of any point T x, y in the disc can be obtained from (19), as shown in Figure 3.

19
where P is the loading pressure (N), D is the diameter of the rock sample (mm), and t is the thickness of the sample (mm).
In the disc center (point o and r 1 = r 2 = 0 5 D and θ 1 = θ 2 = 0), the horizontal and vertical tensile stresses can be expressed as The initial experience failure criteria were proposed by Hoek [32] as follows: where σ ci is the uniaxial compressive strength of the rock (MPa) and m and s are material constants of rock mass.
According to (19), (20), and ( 21), the shear stress is 0 in the disc center, and the stresses σ x and σ y are the minimum principal stress σ 3 and the maximum principal stress σ 1 , respectively.By combining (22), the term for experience failure criteria can be expressed as Because of the tensile strength σ t = σ 3 , ( 23) can be converted to Considering that both m and s are positive for the rock material constants, σ t and σ ci are opposite, and the root is obtained from (24) as follows: The empirical formula for the variation of the compressive strength of different inclination angles of the bedding for intact rocks is established by Jaeger [38] and Donath [39]: where A and B are constants, β m is the inclination angle of the bedding when the sample uniaxial compressive strength is minimum ( °), and β is the inclination angle of the bedding when the uniaxial compression test is conducted ( °).
To investigate tension failure along the inclination angle of the bedding, we can substitute (26) into (25).This yields an equation for the tensile strength of layered shale: Equation ( 27) reflects the characteristics of anisotropic tensile strength of shale, and we can further study the fracture initiation problem of a shale gas well by using it.

The Initiation Pressure Calculation
Model.The fracture initiates where and when the effective tangential stress magnitude equals to the rock tensile strength.It can be expressed as Combining (27), the term for tensile failure from the wellbore can be expressed as Figure 3: Forces acting on a point in the disc.
6 Geofluids Equation ( 29) is absolutely a new model, not only considering the mechanical characteristics of shale when calculating borehole stress but also considering the anisotropic tensile strength when judging rock failure.The fracture initiation model established in this paper is more consistent with the actual situation.

Model Application.
There is a horizontal well in a transverse isotropic shale formation, and the rock mechanical parameters and in situ stress are given in Table 1.The principal stress coordinate system is the same as the geographic coordinate system, and the direction of σ H is along with the X g -axis (north).The horizontal well is drilled along the principal stress of σ h .In this situation, we calculate and analyze the variation of the borehole stress in transverse isotropic shale formation.
In order to observe the relationship between the tangential stress and circumferential angle based on the theoretical model, the relationship curve of the two is shown in Figure 4, and the theoretical model data are all based on Table 1.It can be seen that tangential stress changes periodically with the increase of circumferential angle, and the following analysis will be based on this model.Besides, the transformation curve of the anisotropic tensile strength with the circumferential angle is shown in Figure 5.The tensile strength varies at different circumferential angles, and the fracture is initiated at 85 °instead of 90 °with the initiation pressure 81.24 MPa, which shows the fact that the rock tensile strength near the wellbore has the characteristics of anisotropy.Therefore, based on the model established in this paper with the consideration of transverse anisotropy, the sensitivity factors will be analyzed in the following.1, and other parameters are kept unchanged except for σ v .The ratios of σ v and σ H are, respectively, 1.0, 1.1, 1.2, 1.3, and 1.4 by increasing σ v .The relationship curves between the tangential stress and circumferential angle of different σ v /σ H are shown in Figure 6.Regardless of the ratio of σ v and σ H , the five curves have the same variation at different circumferential angles but have different ranges of the change.With the increase of σ v /σ H , the maximum tangential stress increases, as well as the rate of change.

Analysis of
The relationship curve between initiation pressure and σ v /σ H as shown in Figure 7 proves that σ v /σ H can influence the initiation pressure, and there is a negative correlation relationship between them.With the increase of σ v /σ H from 1.0 to 1.4, the initiation decreased from 77.25 MPa to 54.86 MPa, about 30% reduction, which is    7 Geofluids not a small value.It shows that the stronger the in situ stress anisotropy is, the stronger the tangential stress anisotropy will be, and the fracture is easier to initiate.The value of σ v /σ H equals 1.0 when the in situ stress is assumed to be isotropic, while the stress anisotropy of shale is more and more obvious with the increase of σ v /σ H . Therefore, the analysis of σ v /σ H is of great importance for field construction.

Sensitivity
Analysis of E v /E H .The following calculation data are all based on Table 1 and will be no longer restated.Make E v = 30 MPa unchanged and decrease E H so that E v /E H , respectively, would equal to 1.0, 1.1, 1.2, 1.3, and 1.4.As shown in Figure 8, at the same circumferential angle with the increase of E v /E H , the tangential stress decreases.Although the change law of these curves is the same as well as the circumferential angle of the inflection points, they cannot coincide when these curves are translated along the Y-axis.In other words, the distances between any two curves are unequal at any circumferential angle.
It can be seen from Figure 9 that the initiation pressure decreases with the increase of E v /E H ; the two have negative correlation, but not linear relationship.The initiation pressure decreases to 74.14 MPa when the E v /E H is 1.4, and the rate of decline is getting smaller and smaller.It shows that the fracture will be easy to initiate when the anisotropy of the elastic modulus of shale became stronger.19 and reduce υ H to make υ v /υ H = 1 0, 1.1, 1.2, 1.3, and 1.4.The change law under the four conditions is shown in Figure 10.At the same circumferential angle, the tangential stress increases with the increase of υ v /υ H . Besides, the curves of different υ v /υ H vary at different rates, and in the vicinity of the curve inflection points, the four curves have the largest tangential stress difference.

Sensitivity Analysis of υ
It is obvious from Figure 11 that with the increase of υ v /υ H , the initiation pressure decreases from 80.99 MPa to 79.83 MPa, but the range is very small.Compared with the influence of E v /E H , υ v /υ H has little effect on initiation pressure.8 Geofluids m/s, keep s equal to 0.11 MPa and increase m to make the ratio of m and s equal to 1.0, 3.0, 5.0, 7.0, and 9.0, respectively.It can be seen from Figure 12 that the values of the peak and trough in one cycle are not equal; this also proves the characteristics of tensile strength anisotropy.At any circumferential angle, the higher the m/s is, the smaller the tensile strength is.Initiation pressure and m/s show a nonlinear decreasing relationship in Figure 13.It indicates that the ratio of m and s can affect the anisotropy of tensile strength and further influence the initiation pressure.It also indirectly shows that the inherent properties of rock mass can influence the shale anisotropy.

Sensitivity Analysis of A/B.
Similarly, keep B equal to 48.4 MPa and increase A to make the ratio of A and B equal to 1.0, 1.5, 2.0, 2.5, and 3.0, respectively.The five curves in Figure 14 have exactly the same change laws; the four can be completely coincident by translating along the Y-axis.Besides, the distance between any two curves is equal.
A and B are both constants of the equation for the tensile strength of layered shale, which can characterize the anisotropic tensile strength of shale.Different from the change law mentioned earlier in this paper, the initiation pressure increases from 78.8 MPa to 84.82 MPa with the increase of A/B, and the two show a linear correlation in Figure 15.By observing the variation of initiation pressure with A/B, a more accurate prediction of fracture initiation can be achieved.

Conclusion
(1) A new theoretical model was established to calculate the fracture initiation pressure in this paper with the consideration of mechanical characteristics of shale and the anisotropic tensile strength when judging rock failure.
(2) The sensitive factors σ v /σ H , E v /E H , υ v /υ H , m/s, and A/B have a certain impact on the tangential stress and tensile strength when the circumferential angle changes.
(3) Due to the limitations of the study conditions, this paper does not carry out experimental research, which will also be our future study direction.9 Geofluids (4) Through the newly theoretical model established in this paper, more accurate fracture initiation pressure can be obtained, which is of great significance for judging and controlling the initiation of the fracture in the actual construction process.

Nomenclature σ p :
The far-field stress tensor under the principal stress coordinate system (MPa) σ g : The far-field stress distribution tensor under the geographic coordinate system (MPa) α pg : The azimuth of the maximum horizontal principal stress ( β pg : The angle between the direction of σ v and Z g -axis ( The far-field stresses under the borehole coordinate system (MPa) α bg : The azimuth of the borehole axis ( β pg : The deviation angle of the borehole axis (

°)
Re: The real component of an imaginary number (dimensionless) Φ' i : The derivatives of the stress analytical functions (dimensionless) z i : μ i and λ i complex numbers (dimensionless) E x , E y , E z , E h , and E v : The radius of the borehole (m) θ: The angle measured counterclockwise from the axis X b direction p w : The fluid pressure in the well (MPa) p p : The pore pressure (MPa) α ij : Biot's coefficient tensor (dimensionless) K s : Bulk modulus of the solid phase (GPa) M ij : Stiffness matrix (dimensionless) σ ci : The uniaxial compressive strength of the rock (MPa) m and s: Material constants of rock mass (dimensionless) σ x : Minimum principal stress (MPa) σ y : Maximum principal stress (MPa) A and B: Constants (dimensionless) The inclination angle of the bedding when the sample uniaxial compressive strength is minimum β: The inclination angle of the bedding when the uniaxial compression test is conducted.
,b τ xy,b τ xz,b τ yx,b σ yy,b τ yz,b τ zx,b τ zy,b σ zz,b , R bg = cos α bg cos β bg sin α bg cos β bg sin β bg −sin α bg cos α bg 0 −cos α bg sin β bg −sin α bg sin β bg cos β bg , 2 where σ b is the far-field stress distribution tensor under the borehole coordinate system, σ xx,b , σ yy,b , σ zz,b , τ xy,b , τ xz,b , and τ yz,b are the far-field stresses under the borehole coordinate system, α bg is the azimuth of the borehole axis, and β pg is the deviation angle of the borehole axis.

Figure 1 : 5 μ 2
Figure 1: Coordinate system transformation.(a) The principal stress coordinate system transforms to the geographic coordinate system.(b) The geographic coordinate system transforms to the borehole coordinate system.

Figure 4 :
Figure 4: The relationship between tangential stress and circumferential angle.

Figure 5 :
Figure 5: Transformation curve of the anisotropic tensile strength with the circumferential angle.

Figure 6 :
Figure 6: The relationship curves between tangential stress and circumferential angle of different σ v /σ H .

Figure 8 :Figure 9 :
Figure 8: The relationship curves between tangential stress and circumferential angle of different E v /E H .

Figure 10 :
Figure 10: The relationship curves between tangential stress and circumferential angle of different υ v /υ H .

Figure 13 :
Figure 13: The relationship curve between initiation pressure and m/s.

Figure 14 :
Figure 14: The relationship curves between tensile strength and circumferential angle of different A/B.
The far-field stress distribution tensor under the borehole coordinate system (MPa) σ xx,b , σ yy,b , σ zz,b , τ xy,b , τ xz,b , and τ yz,b : Young's modulus in different orientations of shale formation (GPa) υ xy , υ yz , υ xz , υ h , and υ v : Poisson's ratio in different planes of shale formation (dimensionless) G xy , G yz , G xz , G h , and G v : Shear modulus in different planes of shale formation (dimensionless) a:

Figure 15 :
Figure 15: The relationship curve between initiation pressure and A/B.
and E v are Young's modulus in different orientations of the shale formation, υ xy , υ yz , υ xz , υ h , and υ v are Poisson's ratio in different planes of the shale formation, and G xy , G yz , G xz , G h , and G v are shear modulus in different planes of the shale formation.In (3), λ i can be calculated by a 11 a 12 a 13 a 14 a 15 a 16 a 21 a 22 a 23 a 24 a 25 a 26 a 31 a 32 a 33 a 34 a 35 a 36 a 41 a 42 a 43 a 44 a 45 a 46 a 51 a 52 a 53 a 54 a 55 a 56 a 61 a 62 a 63 a 64 a 65 a 66 Parametric Sensitivity 4.2.1.Sensitivity Analysis of σ v /σ H .The theoretical values are calculated under the condition of Table