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 ( σ_{s}^{i})^{k−1}_{total} and shear displacement discontinuity ( D_{s}^{i})^{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 (σ_{s}^{i})_{0} is the initial shear stress; the sign selection of the yield stress is consistent with the total shear stress:

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 (σ_{s}^{i})_{0} and (σ_{n}^{i})_{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.

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

## 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).

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

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.

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

## 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.

**Table 1.** The default value of the input parameters.

**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.

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.

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

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].

**Figure 9.** Normalized displacement distribution along the natural fracture (NF). (a) Fractures gap ρ = 0.1 m; (b) Fractures gap ρ = 0 m.

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 Figure 12 and Figure 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.

*Figure 10.** Normalized stresses along the NF. (a) Fractures gap ρ = 0.1 m; (b) Fractures gap ρ = 0.0 m.*

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.

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

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.

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

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 13. **Offset of the maximum principal stress peak at different inclinations.*

### 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.

*Figure 14.** Offsets of the maximum principal stresses as a function of hydraulic fracture (HF) net pressure.*

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.

*Figure 15. **Mohr diagrams for the entire sets of elements along the NF. ( a) Moderate strength NF (c = 2.2 MPa, *

*ф*

*= 26.6º); (*

**b**) Weak strength NF (c = 0.0 MPa,*ф*

*= 5.7º).*

## 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.

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

## 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.

*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.*

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.

*Figure 18.** Diagram of the testing specimen.*

* *

*Table 2.** Experimental scheme and parameters.*

* *