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.
Xin Chang1,2, Yintong Guo1,2, Jun Zhou1,2, Xuehang Song2,3 and Chunhe Yang1,2
1State Key Laboratory of Geomechanics and Geotechnical Engineering, Institute of Rock and Soil Mechanics, Chinese Academy of Science, Wuhan 430071, China. 2University of Chinese Academy of Sciences, Beijing 100049, China. 3Shanghai Advanced Research Institute, Chinese Academy of Sciences, Shanghai 201210, China
Received: 8 August 2018 – Accepted: 18 September 2018 – Published: 23 September 2018
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.
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% . 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 .
Moreover, shale gas reservoirs often contain natural fractures. Gale, et al.  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  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  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  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  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  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.  found that fluid tends to divert more into the NF under a lower horizontal stress difference. Zhou, et al.  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.  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.  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  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 recent years, many numerical methods, including finite element methods (FEM) [28,29,30,31,32,33], cohesive zone methods (CZM) [34,35,36,37,38,39,40,41], displacement discontinuity methods (DDM) [21,42,43,44,45,46,47,48,49,50,51,52,53,54,55], discrete element methods (DEM) [56,57], combined finite-discrete element methods (FDEM) [58,59], rock failure process analysis (RFPA) [60,61], and extended finite element methods (XFEM) [62,63,64,65,66] have been proposed to solve the interaction problem between hydraulic and natural fractures.
Additionally, the phase-field approach  and the peri-dynamics method [68,69] were also employed to investigate the interactions between an HF and the NF. Although these advanced models enable us to better understand the hydraulic fracture propagation with complex interactions, very few models have investigated the evolution of stress transfer mechanisms along the natural fracture during mechanical interaction between hydraulic and natural fractures.
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 . It is based on the analytical solution to the problem of a constant dislocation over a finite line segment in an elastic body.
Figure 1. A curved fracture discretized by N elemental displacement discontinuities.
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 σsi and σni are the shear and normal components of stress at the midpoint of the i-th element; usi and uni are the shear and normal displacements of the i-th element; Dsj and Dnj 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 Dsj and Dnj can be written as:
Ass ( i , j ), etc., are the boundary influence coefficients for the stresses, while Bns ( 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 . This failure criterion states that a fracture will start its propagation when the stress intensities satisfy the following inequality:
where KI and KII are the stress intensity factors; KIC is fracture toughness; θ is the fracture initiation angle and satisfy the following equation:
With DDM, the stress intensity factors KI and KII 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; Dn and Ds are the displacement discontinuities of the tip element; C is an empirical correction coefficient, C = 0.806 proposed by Olson  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.
Figure 2. Representation of a joint element.
A joint element with two degrees of freedom, the stress-strain relations can be easily written as:
where Kn and Ks 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 Dni and Dsi are the total deformations of the i-th joint element, which usually be represented as the sum of the initial displacements ((Dni)0 and (Dsi)0) and the induced displacements ( Dn’i and Ds’i ). 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 ci and φi are the cohesive and angle of friction of the joint-filling material.
Figure 3. Mohr diagram of a joint element under different stress conditions.