Numerical and Experimental Investigations of the Interactions between Hydraulic and Natural Fractures in Shale Formations

Natural fractures (NFs) have been recognized as the dominant factors that increase hydraulic fracture complexity and reservoir productivity. However, the interactions between hydraulic and natural fractures are far from being fully understood. In this study, a two-dimensional numerical model based on the displacement discontinuity method (DDM) has been developed and used to investigate the interaction between hydraulic and pre-existing natural fractures. The inelastic deformation, e.g., stick, slip and separation, of the geologic discontinuities is captured by a special friction joint element called Mohr-Coulomb joint element. The dynamic stress transfer mechanisms between the two fracture systems and the possible location of secondary tensile fracture that reinitiates along the opposite sides of the NF are discussed. Furthermore, the model results are validated by a series of large tri-axial hydraulic fracture (HF) tests. Both experimental and numerical results showed that the displacements and stresses along the NFs are all in highly dynamic changes. When the HF is approaching the NF, the HF tip can exert remote compressional and shear stresses on the NF interface, which results in the debonding of the NF. The location and value of the evoked stress is a function of the far-field horizontal differential stress, inclination angle of the NF, and the net pressure used in fracturing. For a small approaching angle, the stress peak is located farther away from the intersection point, so an offset fracture is more likely to be generated. The cemented strength of the NF also has an important influence on the interaction mechanism. Weakly bonded NF surfaces increase the occurrence of a shear slippage, but for a moderate strength NF, the hybrid failure model with both tensile and shear failures, and conversion may appear.


Introduction
Since the beginning of this century, shale gas has become increasingly important source of natural gas supplies.In 2000, shale gas provided only 1% of the natural gas production in the United States, but by 2011 it had become more than 34% [1,2].The US Energy Information Administration (EIA) predicts that by 2035, this proportion will be expanded to 46% [3].Unlike conventional gas reservoirs, shale gas reservoirs are characterized by nanoscale porosity and ultralow permeability, they are not usually associated with natural productivity and have to be economically developed by means of horizontal drilling and hydraulic fracturing [4].Moreover, shale gas reservoirs often contain natural fractures.Gale, et al. [5] has reviewed the types and characteristics of common fractures based on core and outcrop studies from several different shale plays, including Mississippian Barnett shales from the Fort Worth Basin and the Devonian Woodford shale from the Permian Basin.It is found that natural shale fractures are heterogeneous (those at a high angle to bedding, early compacted fractures, bed-parallel fractures, faults and calcite filling cracks).Previously, most hydraulic fractures were considered to be simple and planar features, but recent micro-seismic monitoring in the Barnett shale suggested a more complex fracture network [6][7][8].Therefore, understanding the hydraulic behavior in naturally fractured scenarios is critical to the stimulation design.
As an early attempt, Lamont and Jessen [9] conducted a series of experiments using 106 models made from cement blocks and natural rocks.In their experiments, hydraulic fractures were found to cross closed pre-existing natural fracture at all of the approaching angles and combinations of stresses.Daneshy [10] performed three hydraulic fracturing experiments in granite.He found that when a hydraulic fracture approaching a large open plane of weakness, the hydraulic fracture reoriented itself parallel to the plane of weakness.Anderson [11] studied the frictional properties of unbonded interfaces on the growth of hydraulically driven fractures.Test results showed that the hydraulic fractures will be arrested at the interface when the normal stress below a threshold.Blanton [12] systematically investigated the interactions between hydraulic fracture (HF) and existing natural fractures (NFs) in large scale laboratory experiments.It was found that HF tend to cross the NF only under higher differential stresses and approaching angle.Warpinski and Teufel [13] first discussed the photography and mapping schemes from mine-back hydraulic fracturing experiments.They observed that geologic discontinuities can significant affect the geometry and effectiveness of hydraulic fractures.Beugelsdijk, et al. [14] found that fluid tends to divert more into the NF under a lower horizontal stress difference.Zhou, et al. [15] have focused their attention on the shear strength of the NF.In their experiments, three types of paper with the different coefficient of friction were used to simulate the NFs.Olson, et al. [16] conducted laboratory tests on gypsum cement blocks containing sandstone, glass, and gypsum, to represent the cemented natural fractures.Based on these experiments, complex interactions were observed, e.g., interfacial separation, fracture bypass, diversion, and continued fracture propagation away from natural fracture tip lines.Lee, et al. [17] utilized a semicircular bend test to study the veins on the crack propagation paths.They found that the inclination and thickness of the veins could exert a strong influence on the HF propagation pattern.Alabbad and Olson [18] indicated that the outcome of hydraulic and natural fractures interactions was largely dictated by the mineral cementation inside the natural fractures because the cement fill type determined the cohesion and frictional coefficient along the interface.
Additionally, various analytical criteria have been proposed to predict the behavior of a HF when it intersects a pre-existing NF [12,13,[19][20][21][22][23][24][25][26][27].However, these criteria, which generally consider a constant in-situ stress along the NF and assumed the HF already intersected the NF are oversimplified.In fact, with the extension of the hydraulic fracture, the state of the stress within the rock is in a dynamic changing process, and the new fracture may appear along the posterior side of the NF before the HF intersects the NF.
In this study, the details of the interactions between a natural fracture and a hydraulic fracture were quantitatively studied.A fluid uncoupled model based on the displacement discontinuity method is employed to simulate the HF propagation, initiation, and intersection.The nonlinear deformation, e.g., bonded, sliding, and opening, of the natural fracture is accurately described by the Mohr-Coulomb joint element.The dynamic stress transfer and sliding response of the NF are investigated using a rigorous sensitivity.Furthermore, the probable location of the secondary tensile fracture along the opposite of the NF is determined by the maximum principal stress.Lastly, the model results are validated by a laboratory, true tri-axial hydraulic fracture test.

Displacement Discontinuity Method
The displacement discontinuity method is an efficient and economical numerical method which has been extensively used in fracture mechanics.The two-dimensional displacement discontinuity method was introduced by Crouch [70].It is based on the analytical solution to the problem of a constant dislocation over a finite line segment in an elastic body.
Through discretizing an arbitrarily curved fracture with N elements, as shown in Figure 1, the displacements and stresses along the fracture can be easily formulated: where σ i s and σ i n are the shear and normal components of stress at the midpoint of the i-th element; u i s and u i n are the shear and normal displacements of the i-th element; D j s and D j n are the components of discontinuity in the s and n directions of the j-th element, as shown in Figure 1b, each element has two surfaces, one is on the positive side of n (n = 0 + ) and another is on the negative side, denoted n = 0 − .When crossing from one side of the line element to the other, the displacements undergo a constant change in value.Therefore, the displacement discontinuities D j s and D j n can be written as: A ss (i, j), etc., are the boundary influence coefficients for the stresses, while B ns (i, j), etc., are the boundary influence coefficients of the displacements.Given the prescribed stress or displacement boundary conditions, Equation (1) is a system of equations with 2N unknown displacement discontinuity components.These equations can be solved by standard methods, either direct or indirect.In this study, an iterative method called the "method of successive over-relaxation" is used.Details are provided in Appendix A.

Fracture Initiation and Propagation
The maximum circumferential stress criterion is used to predict the fracture initiation angle and propagation path [71].This failure criterion states that a fracture will start its propagation when the stress intensities satisfy the following inequality: where K I and K II are the stress intensity factors; K IC is fracture toughness; θ is the fracture initiation angle and satisfy the following equation: With DDM, the stress intensity factors K I and K II can be directly calculated by using the fracture tip element displacement discontinuities (DDs): where E is the modulus of elasticity; v is the Poisson's ratio; ∆a is the length of the tip element; D n and D s are the displacement discontinuities of the tip element; C is an empirical correction coefficient, C = 0.806 proposed by Olson [72] is used in this paper.
The fracture propagation path is obtained by adding a new element to the end of the fracture.The length of the newly initiated element is equal to the original fracture tip element.The calculation accuracy can be improved by shortening the element length.

Mohr-Coulomb Joint Element
A special element, referred to as the joint element is used to model the effects of geologic discontinuities, i.e., faults and natural fractures in the rock mass.Unlike the constant displacement discontinuity element, a joint element can be visualized as an elemental displacement discontinuity whose opposite surface are connected by a spring, as shown in Figure 2. A joint element with two degrees of freedom, the stress-strain relations can be easily written as: where K n and K s are the normal and stress stiffnesses of the springs, which are chosen to representative of the properties of the fracture filling material.
If a closed natural fracture is divided by N joint element, the equations between tractions and the DDs can be written as: where D i n and D i s are the total deformations of the i-th joint element, which usually be represented as the sum of the initial displacements ((D i n ) 0 and (D i s ) 0 ) and the induced displacements (D i n and D i s ).In this study, the initial deformations of the natural fracture were not considered, so the total deformations of the joint element is equal to the induced displacements.
In reality, there must be a constraint on the maximum shear stress that can be developed across the joint.Such a constraint is given by the Mohr-Coulomb condition, as shown in Figure 3: where c i and ϕ i are the cohesive and angle of friction of the joint-filling material.When a joint element satisfies the constraint expressed by Equation ( 7), it will be called a Mohr-Coulomb element.The inelastic deformation of the natural fracture can be accurately simulated by using this special element.It is required that the element be allowed to undergo a certain amount of inelastic deformation or permanent slip, when the total shear stress on the joint element exceeds the total yield stress.Suppose that loading increment k − 1 has been completed, and no inelastic deformation occurred during this or any previous increment, the total shear stress σ i s k−1 total and shear displacement discontinuity D i s k−1 will lie at some point A (see Figure 4).At the next loading increment, the i-th element may load or unload elastically (points B and C), or it may fail in accordance with Equation (1) (points D and E).If the magnitude of the shear stress exceed the value of the yield stress σ s,yield from Equation ( 7), the total shear stress must equal the yield stress, and so the corresponding equations can be modified as follows: where σ i s 0 is the initial shear stress; the sign selection of the yield stress is consistent with the total shear stress: sign σ i s,yield = sign σ i s total (9) In the above discussion, we have tacitly assumed that the normal stress across a joint element is compressive.According the Mohr-Coulomb condition (Equation ( 7)), if the stress greater than c cot ϕ as shown in Figure 3, another failure model is possible, namely joint separation or crack open.In this study, we assume that natural fractures cannot bear tensile forces.Therefore, if the normal stress on a joint element changes from compressive stress to tensile stress, an opened Mohr-Coulomb element is formed.For such an element, the boundary conditions are that the total normal and shear stresses are zero.In this case, a cracked Mohr-Coulomb element is equivalent to a constant displacement discontinuity element, so the following equations must be satisfied: where σ i s 0 and σ i n 0 are the initial shear and normal stresses.

Numerical Procedure
The occurrence of slip along a joint is a nonlinear, path-dependent phenomenon, which must be modeled by an incremental process.In this paper, the final state is reached after K increments.The detailed iterative procedure is shown in Figure 5.

Model Verifications
The accuracy of the numerical model is validated by calculating the stress distribution along a bonded horizontal interface just ahead of a vertical opening fracture.The geometry of the model is shown in Figure 6.For simplicity, a uniform remote tensile (5 MPa) is used.The distance between the interface and the fracture tip is s; the vertical fracture length is 2a, which is 1-m-long; Young's modulus is 20 GPa; Poisson's ratio is 0.25.The analytical solutions can be easily derived from the Westergaard function (Appendix B).
As shown in Figure 7, a good agreement between the numerical results and analytical solutions is displayed.The small deviation may arise owing to the numerical discretization of the model and the use of the constant displacement discontinuity elements.

Case Studies
In this section, a series of numerical experiments have been performed to investigate the interaction between a hydraulic fracture and a pre-existing natural fracture.The geometry of the model is illustrated in Figure 8, a HF propagates along the direction of maximum horizontal stress and ultimately intersects an inclined NF.For simplicity, a fluid-uncoupled model is used herein owing to the desire specify a prescribed uniform pressure distribution along the HF.The parameters selected for these analyses are listed in Table 1.In the process of numerical calculation, the HF was discretized with 160 constant displacement discontinuity elements, each 1.25 cm wide.In addition, the NF was divided into 400 Mohr-Coulomb elements with an equal length of 0.5 cm.
To simplify the expression, some dimensionless parameters are introduced: where x and y are the dimensionless coordinates; ρ are the normalized gap between the two fractures; D s,n are the normalized displacement discontinuities; ∏ is the normalized net pressure; ∆ is the normalized horizontal stress difference; σ m and σ d are the average and difference of the horizontal stresses, respectively.

Displacements and Stresses along the Natural Fracture (NF)
The normalized shear and normal displacements along the NF with different fractures gap are shown in Figure 9.The results clearly indicated that the NF has been activated before intersecting with the HF, and the magnitudes of the displacements rapidly increase as the fractures gap decrease.In this case, the peak of normalized normal displacement is −0.27 when the fractures coalesce, which is three times bigger than that of the default fractures gap.The normalized stresses along the NF are shown in Figure 10.From this figure, we can see that near the center of the NF, there is an opening zone where the shear and normal stresses are equal to zero, and the maximum principal stresses are located at the end of this zone.Meanwhile, the maximum principal stress has two different peak points, and the bigger one is located at the positive side of the NF, so secondary tensile cracks are more likely appear at this location.This result can be confirmed by the field test of Jeffrey [73].
Then, we analyse the distribution of the normalized maximum principal stress along the NF with different fractures gaps.As shown in Figure 11, we can see that with a decreasing fracture gap, the maximum principal stress increases dramatically, and the location of the peak is farther from the junction.Note that the induced maximum principal stress along the opposite side of the NF may exceed the tensile strength of the rock mass before the fracture coalescence.

Parametric Study of Stress Distribution and Unstable Tensile Zone along the NF
Firstly, the natural fracture inclination is considered.As plotted in Figures 12 and 13, we can see that NF slope angle has a significant impact on the stress distribution along the NF.With the increasing of the NF inclination, the magnitude of the maximum principal stress rapidly increases, but the peak value location is farther from the central of the NF.It means that the HF is prone to directly cross the high angle natural fracture.Secondly, we focused our attention to HF net pressure Π and horizontal stress difference ∆.As shown in Figure 14, we can see that the offset of the maximum principal stress peak increases proportionally to the HF net pressure.It means that we can increase the fracture complexity by increasing the HF net pressure.However, the horizontal stress difference may have a negative effect on the fracture complexity.In this case, the offset of the maximum principal stress peak in a nonuniform in-situ stress field is smaller than that under a uniform stress.

Effects of Natural Fracture on the Hydraulic Fracture (HF) Propagation Path
The accurate prediction and monitoring of hydraulic fracturing propagation pathways is very important for reservoir stimulation practice.Therefore, the flowing numerical simulations are conducted to investigate the influence of NF on the HF trajectories.The main parameters required for the model are derived from [44].One issue that needs to be especially pointed out is that the initial deformation of the natural fracture is also not considered in this case.The HF trajectories for different NF inclination angles are depicted in Figure 16.From this figure, we can see that the NF inclination angle has an important effect on the shape of the HF propagation path.When the NF inclination is 90 • , the HF extends parallel to the x-axis until it crosses the NF.However, as the slope angle decreases, the left side of the NF (i.e., y < 0) is closer to the HF surface, and the direction of the maximum principal stress in front of the NF will change.Thus, the HF is more prone to deviate from its original propagation path and reorientation near the NF.

Experimental Investigation
For better understanding of the mechanics of how the NFs affect the HF initiation and propagation and validation of the numerical model results, a series of laboratory stimulated experiments were carried out on the Longmaxi shale outcrops collected from Sichuan Base, China.

Experimental Setup and Sample Preparation
The apparatus used in this study was a large-scale, true tri-axial hydraulic fracturing test system (as shown in Figure 17), capable of subjecting an 800 × 800 × 800 mm block to magnitudes of confining stress up to 30 MPa.The hydraulic fluid injection rate and pressure were precisely controlled by an individual hydraulic pump pressure servo control system.The maximum hydraulic fluid injection pressure could reach values up to 100 MPa.Apart from these variables, an acoustic emission monitoring system (that contained eight probes, operated within a frequency range of 15-75 kHz, with resonant frequencies of 40 kHz, and a minimum noise threshold of 18 dB), was used to detect the acoustic signal during fracture initiation and propagation.Two cube samples of 300 mm × 300 mm × 300 mm were prepared for these experiments.Each block had a simulated wellbore (diameter: 24 mm; length: 170 mm) in its center, with the axis of the borehole parallel to the minimum horizontal stress.Two slotted steel pipes (diameter: 20 mm; pipe thickness: 2.5 mm) were then cemented with epoxy onto the simulated wellbore, leaving a 5 mm open-hole section at the lower side of the wellbore.The direction of wellbore and in-situ stresses are shown in Figure 18.For each experiment, water with red fluorescent dye was injected into the wellbore.The experimental scheme is shown in Table 2.

Experimental Results and Analyses
Photographs of the cube after the fracturing tests are shown in Figures 19 and 20.It is found that the geometry of the shale hydraulic fracture is more complex, which is different from the traditionally bi-winged, planar fracture.The morphology of the HF is strongly influenced by the opened NFs and sedimentary beddings.The crack initiation and evolution can also be analyzed by the acoustic emission (AE) detection results during experiments.One red dot is defined as an AE event.For specimen Y-7-1 (shown in Figure 21), the HF firstly began to propagate along the overburden stress, lately extended to an inclined pre-existing natural fracture (fracture #1).Shortly, some acoustic emission events were gradually appeared along the natural fracture (fracture #2).These AE events were most probably generated by shear slippage along the NF surface since fluid that leaked off increased the pressure and reduced the cemented strength.As time elapsed, the fluid pressure was sufficient to create a new secondary fracture at the weak location on the surface with the natural fracture.Thus, a step-like fracture with an obvious lateral offset was formed in the vertical direction (fracture #3).At the same time, some additional secondary fractures near the crack junction were also formed by the higher fluid pressure (fracture #4).Finally, the hydraulic fracture was arrested by a partially opened bedding fracture at the upper side and extended to the outer surface of the specimen at the lower left corner.
For another specimen, labeled Y-7-3 (shown in Figure 22), the HF firstly propagated along the direction of the overburden stress (fracture #1), and thus turned into an opened bedding fracture (fracture #2).Soon afterwards, the partially opened bedding fracture near the wellbore was activated by the higher fluid pressure (fracture #3).As time elapsed, another pre-existing bedding fracture (fracture #4) was also activated by the main hydraulic fracture (fracture #1).Finally, a discrete fracture network was formed, but the complexity of the final fracture network was lower than that of specimen Y-7-1.The variation curves of the pump pressure and total injected volume during the tests are shown in Figure 23.It can be observed that there is a significant difference between the pump pressure curves of the two model blocks.In test a typical pump pressure curve was recorded, and an obvious fracture breakdown was indicated by a sharp decline of the pump pressure.However, for specimen Y-7-3, it was difficult to see obvious fracture breakdown characteristics from the pump pressure curve, and the pump pressure continued to rise until a stable hydraulic flow channel had been fully formed.This was because the hydraulic fractures were mainly propagating along the opened bedding fracture, and a higher fluid pressure was needed to overcome the overburden stress and stress shadow.In addition, the larger total injected volume and less AE events indicate that the specimen Y-7-3 was dominated by tensile failure and the sample Y-7-1 was mainly present in the form of shear failure.Based on these experimental results, we found that the morphology of the shale hydraulic fracture was strongly influenced by the characteristics of the pre-existing natural fractures and sedimentary beddings.The density, orientation, and cemented strength of the pre-existing natural fractures are the three factors dominating the formation of complex fracture networks.The hydraulic fractures are either diverted or arrested by the NF.The experiment results are in good agreement with the numerical results.

Conclusions
Based on these numerical and experimental analyses, the following conclusions can be proposed: (1) The hydraulic fracture morphology of shale was strongly influenced by the characteristics of its natural fractures.The NF density, orientation, and cemented strength, are three main factors dominating the formation of complex fracture networks.(2) The cemented strength of the NF caused a noticeable nonlinear behavior for this problem.For an NF with a low strength, only shear failure occurred, and the HF was more likely to terminate.However, for an NF with a moderate strength, the hybrid failure model (tensile failure and shear failure co-existence, and conversion) may occur, and the HF is more inclined to step over at the contact.(3) The displacements and stresses along the NF all changed in a highly dynamic manner.During the stage of approach of the HF to an NF, the HF tip could exert a remote compressional and shear stress on the NF interface, which could lead to the debonding of the natural fracture.Meanwhile, a maximum principal stress peak is generated at the end of the opening zone, where a new tensile crack is more likely to occur.(4) The location and value of the stress is the function of NF inclination angle, far-field differential stress as well as HF net pressure.For small approaching angles, the stress peak is located farther from the intersection point, so a step over (offset) fracture is more likely to occur.The same effect can be found by a higher HF net pressure, because the stress perturbation ahead the NF is proportional to the fluid pressure in the HF.(5) Meanwhile, the HF was found to be more prone to deviate from its original propagation path and reoriented near the NF because the deformation of the NF.
If we denote the k-th approximations of D i s and D i n by D i s k and D i n k , the following iteration formulas can be written as: In these equations, ω is the over-relaxation factor (1 ≤ ω < 2).Omitting the intermediate steps of the computations, ∆D i s k+1 and ∆D i n k+1 may be written with the aid of Equation (A2) as:

Appendix B
One of the most typical crack problems in fracture mechanics is an infinite plane with a line crack of length 2a subjected to biaxial stress σ 0 at infinity, as shown in Figure A1.The solution to the preceding boundary value problem was firstly given by Westergaard [75].To find the explicit expressions of the stress field, one usually uses the polar coordinates.Therefore, the stress components can be obtained: In the vicinity of the crack tip, we have: Along the crack extended line (θ = θ 1 = θ 2 = 0), these near-tip stresses are:

Figure 1 .
Figure 1.A curved fracture discretized by N elemental displacement discontinuities.

Figure 2 .
Figure 2. Representation of a joint element.

Figure 3 .
Figure 3. Mohr diagram of a joint element under different stress conditions.

Figure 4 .
Figure 4.The curve of shear stress and shear deformation.

Figure 5 .
Figure 5. Flowchart of the computing procedure.

Figure 6 .
Figure 6.Model used to simulate a vertical opening fracture approaching a horizontal bonded interface.

Figure 7 .
Figure 7.Comparison of the analytical solution of the stresses with numerical results.(a) Shear stress when gap is 0.5 cm; (b) Maximum tension under different gap.

Figure 8 .
Figure 8. Geometry of hydraulic and inclined natural fractures.

Figure 11 .
Figure 11.Maximum principal stresses along the NF with different fracture gaps.

Figure 12 .
Figure 12.Maximum principal stresses along the NF at different inclinations.

Figure 13 .
Figure 13.Offset of the maximum principal stress peak at different inclinations.

Figure 14 .
Figure 14.Offsets of the maximum principal stresses as a function of hydraulic fracture (HF) net pressure.Finally, the cemented strength of NF is being considered.The Mohr diagrams for the entire set of elements along the NF are depicted in Figure 15.The two inclined straight lines in each diagram represent the yield condition and are usually called Mohr-Coulomb envelope lines.The bracketed numbers refer to the element numbers.From this figure, it can be observed that for a moderate strength NF, the shear and tensile damages occur simultaneously.In this case, the elements (No. 201-252) in the range of 0.0 m < y < 0.26 m open, while the elements (No. 253-262) in the range of 0.26 m < y < 0.31 m undergo sliding.However, for a weak strength NF, only shear failure occurs.The elements (No. 202-365) in the range of 0.0 m < y < 0.82 m undergo sliding, which means that the HF is more likely to turn toward the NF direction.

Figure 16 .
Figure 16.HF propagation path at different NF inclination angles.

Figure 17 .
Figure 17.The main devices of the true tri-axial fracturing system.(a) Tri-axial loading system; (b) Force transmission frame; (c) Reserved holes of the acoustic emission (AE) sensor.

Figure 18 .
Figure 18.Diagram of the testing specimen.
) and (A4) constitute the iterative solution of Equation (A1).The convergence of the iteration process is said when the absolute values of ∆D i s k+1 and ∆D i n k+1 are less than or equal to a specified small number ε for all of the elements.

1 σ
of the j-th element (L) D j n normal displacement discontinuity of the j-th element (L) D s,n normalized displacement discontinuities E modulus of elasticity (ML −1 T −2 ) K I stress intensity factor of kind I (ML −1/2 T −2 ) K II stress intensity factor of kind II (ML −1/2 T −2 ) K IC fracture toughness (ML −1/2 T −2 ) K s shear stress stiffness (ML −2 T −2 ) K n normal stress stiffness (ML −2 T −2 ) u i s shear displacements of the i-th element (L) u i n normal displacements of the i-th element (L) ν Poisson's ratio x dimensionless coordinate in the x direction y dimensionless coordinate in the y direction σ m average of the horizontal stresses (ML −1 T −2 ) σ d difference of the horizontal stresses (ML −1 T −2 ) σ i s shear components of stress at the midpoint of the i-th element (ML −1 T −2 ) σ i n normal components of stress at the midpoint of the i-th element (ML −1 T −2 ) yield shear stress (ML −1 T −2 ) θ fracture initiation angle ρ normalized gap between the two fractures ∏ normalized net pressure ∆ normalized horizontal stress difference

Table 1 .
The default value of the input parameters.

Table 2 .
Experimental scheme and parameters.