Investigation of Processes of Interaction between Hydraulic and Natural Fractures by PFC Modeling Comparing against Laboratory Experiments and Analytical Models

Abstract

Hydraulic fracturing technology is usually used to stimulate tight gas reservoirs for increasing gas production. The stimulated volume depends in part on the pre-existing natural fractures in a reservoir. The mechanisms influencing the interaction between hydraulic fractures and natural fractures have to be well understood in order to achieve a successful application of hydraulic fracturing. In this paper, hydraulic fracturing simulations were performed based on a two-dimensional Particle Flow Code with an embedded Smooth Joint Model to investigate the interactions between hydraulic fractures and natural fractures and compare these against laboratory experimental results and analytical models. Firstly, the ability of the Smooth Joint Model to mimic the natural rock joints was validated. Secondly, the interactions between generated hydraulic fractures and natural fractures were simulated. Lastly, the influence of angle of approach, in situ differential stress, and the permeability of natural fractures was studied. It is found that the model is capable of simulating the variety of interactions between hydraulic fractures and natural fractures such as Crossed type, Arrested type and Dilated type, and the modeling examples agree well with the experimental results. Under high approach angles and high differential stresses, the hydraulic fractures tend to cross pre-existing natural fractures. Under contrary conditions, a hydraulic fracture is more likely to propagate along the natural fracture and re-initiate at a weak point or the tip of the natural fracture. Moreover, these numerical results are in good agreement compared with Blanton’s criterion. The variety of permeability of natural fractures has a great effect on their interactions, which should not be overlooked in hydraulic fracturing studies.

1. Introduction

Since shale gas reservoirs often have very low permeability [1], the technique of hydraulic fracturing is commonly used to enhance shale gas production. The successful application of hydraulic fracturing in low permeability reservoirs has resulted in a rapid increase of natural gas production in recent years. Studies of Barnett shale have demonstrated that natural fractures are common, narrow, sealed by calcite, and present in en-echelon [2]. Although the sealed narrow fractures cannot contribute to the reservoir storage, they act as weakness planes and reactivate during hydraulic fracturing treatments, which can produce additional flow paths from formation to borehole to enhance the reservoir permeability. Field microseismic monitoring observations during hydraulic fracturing show that the reactivation of the natural fracture network improves the efficiency of the stimulation [3,4,5]. It is generally believed that the presence of natural fractures will result in higher complexity of the hydraulic fracture network because of the interactions between hydraulic fractures and natural fractures. Due to the variety of the mechanical properties of rock the orientation of natural fractures, in-situ stress states, and the hydraulic fracturing treatment parameters, the interaction is undetermined, and thus the propagation path of the hydraulic fracture is undecided. It is essential to understand more details about the hydraulic fracturing processes, including the interaction between hydraulic fractures and natural fractures, and fracture patterns to achieve a complex fracture network.

In the past few years, the interaction between hydraulic fracture and natural fractures has been investigated with of laboratory experiments, analytical studies, and numerical simulations [6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27]. The main results are summarized in the following.

Experimental work is a direct way to study the interaction between hydraulic fractures and pre-existing fractures. Lamont and Jessen [6] conducted a series of experiments to study the fracture crossing phenomenon and found that all existing fractures with an angle of inclination from 0° to 45° can be crossed by natural fractures. Daneshy [7] stated that existing flaws, the dimension of which are small compared to the driven hydraulic fractures, have no effect on hydraulic fracture propagation, while the interactions between larger-scaled hydraulic fractures and natural fractures are complex. In addition, he stated that the weakest planes encountered in oil reservoirs are of the closed type, which can be crossed by hydraulic fractures. In the laboratory experiments on fractured Devonian shale and hydrostone, Blanton [9] observed that the induced fracture either crossed the natural fracture or opened it, partly depending on the approach angle and differential stress. These experiments suggest that, with high approach angles and high differential stress, crossing is favored and otherwise opening. Warpinski and Teufel [10] studied the effect of single artificial fractures on extending hydraulic fractures and suggested that hydraulic fractures crossed the artificial fractures at angles of 60° or higher under the high differential stress of 10 MPa. Zhou et al. [12] performed a number of hydraulic fracturing experiments in which three different surface roughness papers emulated the pre-existing fractures with different frictional coefficients. Three scenarios such as arrest, dilatation, and crossing were observed, depending on the approach angles, in situ differential stress, and the shear strength of natural fractures. Liu et al. [14] concluded based on experiments that natural fractures propagate following the most preferential propagation, the least resistance, and the shortest propagation path.

Based on the aforementioned experimental results, various analytical criteria have been developed to predict the propagation direction of induced fractures [9,10,11,13]. Sarmadivaleh and Rasouli [28] presented a summarization of the analytical criteria for interaction between hydraulic and natural fractures, including the criteria of Blanton, Warpinski and Teufel, and Renshaw and Pollard, as well as the modified criteria of Renshaw and Pollard. In Blanton’s criterion [9], under various differential stresses and angles of approach, if the pressure at the intersection point exceeds the normal stress, the natural fracture will open, while the hydraulic fractures cross the natural fractures and re-initiate on the other side if the required pressure for re-initiation is less than the opening pressure. In Warpinski and Teufel’s criterion [10], the pore pressure distribution is considered to predict the dilatation or shear slippage of natural fractures. Renshaw and Pollard [11] developed a criterion for predicting whether a hydraulic fracture will orthogonally cross a frictional interface based on the linear-elastic fracture mechanics solution for the stresses near the fracture tip. Renshaw and Pollard’s criterion assumes that crossing-type interaction will not occur if a slip takes place, which is in contrast to some experimental results [12,14]. Gu and Weng [13] extended Renshaw and Pollard’s criterion to determine whether a fracture crosses a frictional interface at nonorthogonal approach angles. There are some basic assumptions among these analytical models such as the flat geometry of a single natural fracture, its uniform roughness, and the impermeable sides of the natural fracture. However, the analytical models are useful to obtain an initial knowledge of the interaction mechanism.

With the rapid development of computer science, numerical modeling becomes a very important method in research. Lots of numerical models have been developed to simulate fluid-driven complex fracture network formation in the subsurface. To deal with more simple problems, the Perkins-Kern-Nordgren model (PKN) [29,30] and the Khristianovic-Geertsma-de-Klerk model (KGD) [31,32] have been developed. McClure [33] investigated the fracture propagation in a pre-existing discrete fracture network based on the two-dimensional displacement discontinuous method (2D-DDM). Wu and Olson [34], using the modified 2D-DDM approach, studied the influence of stress shadow on the fracture extension pattern when considering several perforation clusters generated in a horizontal well. Chen [35] used the commercial finite element code ABAQUS, wherein the stress intensity model is replaced by a cohesive-zone fracture tip model, to simulate fracture propagation driven by fluid, and their modeling results are in agreement with the 2D PKN and KGD solutions. Wang [36] also used an ABAQUS embedded extended finite element method and a cohesive zone method to model fracture initiation and propagation in different brittle rocks, and the results show that the fluid pressure field and fracture geometry are significantly affected by the inelastic deformations of rock. However, this kind of simulation is still a challenging problem because the complex coupling mechanisms between fluid and solid, the heterogeneity of reservoirs, and the changing boundary conditions have to be taken into account.

In recent years, the discrete element method (DEM) [37] has been developed to investigate hydraulic fracture growth in naturally fractured reservoirs. Choi [38] used the Universal Distinct Element Code (UDEC) to compare several methods for determining the shut-in pressure during hydraulic fracturing because, owing to the relationship between the behavior of hydraulic fractures and remote in situ stress, the shut-in pressure is not considered certain. Zangeneh et al. [39,40,41] further demonstrated the potential of UDEC for simulating hydraulic fracturing, compared different injection schedules from adjacent wellbores, modeled fault slip, and simulated the interaction of hydraulic fractures and natural fractures. Hamidi and Mortazavi [42,43] used 3-Dimension Distinct Element Code (3DEC) for small models to study the effects of different fracture fluid properties and fluid rates, in situ stress states, and the mass properties of rock on hydraulic fracture propagation.

The Particle Flow Code (PFC) is another typical DEM based on the bonded particle method (BPM) [44]. Compared with UDEC, PFC modeling of hydraulic fracturing allows not only fluid flow through fractures, but also fluid leak-off into the rock matrix. In recent years, PFC has been widely applied in geomechanics to model the failure processes of geomaterials and hydraulic fracturing [18,23,24,45,46,47,48,49,50]. Zhao and Young [18] validated the 2-Dimension Particle Flow Code (PFC^2D) for modeling hydraulic fracturing. Shimizu et al. [47] conducted a series of hydraulic fracturing simulations in competent rock by using flow-mechanically coupled PFC^2D code and investigated the influence of the fluid viscosity and particle size distribution. Their results show that tensile cracks are dominantly generated in hydraulic fracturing process, while the energy from shear type acoustic emission is larger than from tensile types. After fracture creation, the velocity of the fluid infiltrating the fracture highly depends on the viscosity of the fluid. Eshiet et al. [48] employed PFC^2D to model the pressure development and the subsequent fracturing and/or cavity propagation of bulk rock and sand, respectively, considering fluid flow based on the Navier-Stokes equation of an incompressible fluid with contact density. Yoon et al. [49] used PFC^2D to study seismicity during the development of an Enhanced Geothermal System by fluid injection deep underground and simulated multistage hydraulic fracturing processes in a fractured reservoir. They then investigated the optimization of a hydraulic fracture network according the stress shadow [50]. Zhou et al. [23,24,25] studied hydraulic fracturing processes in isotropic and laminated reservoirs based on a modified fluid-mechanically coupled mechanism, and this model has been validated by comparing numerical values of breakdown pressures against analytical solutions. In addition, some influence factors on the driven fracture patterns and geometries such as the fluid viscosity, fluid injection rate, and the in situ stress state were analyzed.

Although the PFC has been used in hydraulic fracturing studies as mentioned above, the research into the details of interactions between hydraulic fractures and natural fractures is less reported. In this paper, using the laboratory experimental results of Zhou et al. [12] as a comparison, the details of interactions between hydraulic fractures and natural fractures are studied by PFC^2D based on our fluid-mechanically coupled code [23,24]. Firstly, the simulation method is introduced and summarized, including the Particle Flow Code, the Smooth Joint Model, and the fluid-mechanical coupling mechanism in bonded granular matrix. Secondly, the numerical model is set up and the microscale parameters are calibrated according to the mechanical properties of the laboratory sample (cement plaster). Thirdly, a series of numerical simulations is carried out for comparison with the experimental results and to analyze the influence of approach angle and differential stress on the interaction between hydraulic and natural fractures. Finally, the influence of the permeability of natural fractures on the interaction is analyzed.

2. Simulation Methodology

2.1. Particle Flow Code

PFC^2D is a distinct element modeling method in which the solid materials are represented as an assembly of circular particles. Although PFC^2D is based on the discontinuum method, with the help of bond models at the contacts between round particles, it can also be used to model the deformation behavior of the continuum. With the properties of normal and shear stiffness, as well as shear and tensile strength, the bonds can simulate deformation and micro-crack development based on the relationship presented by Potyondy and Cundall [44]. In rock mechanics research, the parallel bond model is among the most frequently used models. The corresponding microscale properties and deformation and failure behaviors are presented in Figure 1. Normally, the Young’s modulus E of an emulated rock sample is related to the specified contact micro-stiffness. Poisson’s ratio ν is affected by the ratio of normal to shear stiffness. The microscale parameters in this method are different from macroscale parameters such as E and ν, which can be directly measured in the laboratory. These parameters have to be calibrated according to the laboratory results from the confined/unconfined compressive test and the direct tensile test. Generally, by adjusting the micro-stiffness and micro-strength of particles, a realistic rock can be reproduced. Under the applied load, when the maximum shear stress or tensile stress acting on the bond exceeds the shear or tensile strength, the bond will break in shear or tensile mode, resulting in shear or tensile micro-cracks, respectively. With the ongoing generation of micro-cracks, a macro-fracture can be formed by the linking of these individual micro-cracks.

2.2. Smooth Joint Model (SJM)

The general method to emulate natural joints in PFC is the bond removal method, in which the particle contacts lying on a joint track are left unbonded. This approach has been used to study the shear behavior of rock joints in a number of studies [52,53,54]. According to Bahaaddini et al. [51], the ability of the bond removal method to reproduce the shear behavior of rock joints is limited because the circular shape, unequal size of the particles, and non-uniform distribution lead to an unrealistic shear behavior of the rock joints. In order to overcome the shortcomings of this approach, Pierce et al. [50] introduced the Smooth Joint Model (SJM) into PFC.

The SJM emulates the behavior of joints by micro-scale slip surfaces at the contacts between particles that lie opposite to each other at the joint interface, as shown in Figure 2a. At these contacts, the parallel bonds are removed and the new bonding model is assigned with pre-defined orientation, whereas contacting particles can overlap and pass through each other [55]. The details of the fundamental algorithm of the SJM can be acquired from the software manuals [56], while a summary of this model in PFC^2D code is given here.

In SJM, the force F and relative displacement U at a contact point can be given as follows [56]:

$$F=F_n{\widehat{n}}_j+F_s$$

(1)

$$U=U_n{\widehat{n}}_j+U_s$$

(2)

where ${\widehat{n}}_j$ is the normal unit vector at contact, as shown in Figure 2a. The positive values of F_n and U_n represent the compression force and the overlapping of particles, respectively, while F_s and U_s represent the shear force vector and the shear displacement vector, respectively. The force-displacement relationship of each smooth joint contact point follows the Coulomb sliding model with dilation, as shown in Figure 2b. Micro-scale parameters such as the SJM normal stiffness ${\overline{k}}_{nj}$, SJM shear stiffness ${\overline{k}}_{sj}$, SJM friction coefficient ${\mu }_j$, and SJM dilation angle ${\phi }_j$ reflect the mechanical behavior of a smooth joint contact with the area of cross Section A.

The increments of joint force are calculated by the elastic components of the displacement increment, multiplying the SJM normal and shear stiffness. The normal force and shear force are updated as follows [56]:

$$F_n:=F_n+{\overline{k}}_{nj}A\mathrm{\Delta }{U}_n$$

(3)

$${\mathbf{F}}_{\mathrm{s}}^{\mathbf{\prime }}:=F_{\mathrm{s}}-{\overline{k}}_{sj}A\mathrm{\Delta }{U}_{\mathrm{s}}$$

(4)

If $\left|{\mathbf{F}}_{\mathrm{s}}^{\mathbf{\prime }}\right|\le (F_s^{\ast }={\mu }_j{F}_n)$, then $\left|F_{\mathrm{s}}\right|=\left|{\mathbf{F}}_{\mathrm{s}}^{\mathbf{\prime }}\right|$. If $\left|{\mathbf{F}}_{\mathrm{s}}^{\mathbf{\prime }}\right|>F_s^{\ast }$; then sliding occurs at the smooth joint contact, and the shear force and normal force are updated as follows [56]:

$$\left|F_{\mathrm{s}}\right|=F_s^{\ast }$$

(5)

$$F_n:=F_n+[\mathrm{\Delta }{U}_s^{\ast }\tan {\phi }_j]{\overline{k}}_{nj}A=F_n+(\frac{\left|{\mathbf{F}}_{\mathrm{s}}^{\mathbf{\prime }}\right|-F_s^{\ast }}{{\overline{k}}_{sj}}){\overline{k}}_{nj}\tan {\phi }_j$$

(6)

In order to investigate the ability of the smooth joint model to reproduce the shear behavior of rock joints, a series of direct shear tests was simulated. Figure 3 shows the model in which the length and width are 0.4 m and 0.4 m, respectively. A joint modeled by SJM was generated in the model center, which cuts the model into an upper and a lower part. The lower part was fixed by three immovable walls, and constant normal stress was applied on the top wall of the upper part, while a constant horizontal velocity was applied to the two side walls of the upper part. The shear displacement was recorded at the center of the left wall of the upper part, and the shear force was recorded according to the unbalanced force on the right wall of the lower part. The coefficient of friction of the joint was set to 0.38, and its cohesion was 0. The direct shear test results are shown in Figure 4. Under four different normal stresses, the shear stress increases up to a stable value, with shear displacement increasing. The stable shear stresses are 1.14 MPa, 1.90 MPa, 2.66 MPa, and 3.80 MPa, corresponding to the normal stress of 3 MPa, 5 MPa, 7 MPa, and 10 MPa, respectively (Figure 4a). The linear regulation relationship between shear stress and normal stress shows that the friction coefficient is 0.38 (Figure 4b), which is the same as the input. According to the test results, it can be concluded that the SJM can reproduce the shear behaviors of rock joints well.

2.3. Fluid Flow in a Bonded Granular Matrix

The original concept of the fluid flow algorithm in a bonded granular matrix is based on the network flow model as first proposed by Cundall [45]. In Cundall’s fluid flow algorithm there are two assumptions. One is that the fluid domain is surrounded by contacted round particles, and the other is that these fluid domains are linked by flow channels at the contacts. This algorithm has been used in a number of studies [18,23,24,25,45,46,47,48,49,50], where more details can be found. The summarized equations will be given in the following.

Under differential pressure, the fluid flows from one element domain to another, which is assumed to be a process in which the fluid passes two parallel plates in laminar flow and is modeled by the cubic law. The volumetric flow rate q is given as follows:

$$q=\frac{e^3}{12\mu }\frac{\mathrm{\Delta }P}{L}$$

(7)

where e is hydraulic aperture, $\mathrm{\Delta }P$ is the fluid pressure difference between the two neighboring reservoir domains, L is the length of the flow channel, and μ is the fluid dynamic viscosity. The out-of-plane thickness is assumed to be of unit length.

During the fluid flow calculation, the increment of fluid pressure ($\mathrm{\Delta }P$) in a reservoir is computed from the bulk modulus of fluid ($K_f$), the volume of the domain ($V_d$), the sum of the flow volume for a time step ($\mathrm{\Delta }t$), and the volume change of the domain due to mechanical loading, which is neglected in some studies [31,32] but not in this study. Thus the equation used is given by

$$\mathrm{\Delta }P=\frac{K_f}{V_d}\left({\displaystyle \sum q\mathrm{\Delta }t}-\mathrm{\Delta }{V}_d\right)$$

(8)

Reservoir deformations are caused by the fluid pressure exerted on the surrounding particles (see Figure 2b in Zhou et al. [24]). This force (f) is a product of fluid pressure (P_f), the length (l) exposed to the fluid in a domain, and the unit thickness (1 m) in an out-of-plane direction. The resulting force is then applied to the particles in an outward direction normal to the boundary exposed to the fluid.

In the present model, the hydraulic aperture e is not the real aperture of the fluid flow channel but an empirical value based on the rock permeability test under different confining pressures. As shown in Equation (9) [49,50], the hydraulic aperture e of each flow channel is given as function of normal stress ${\sigma }_n$ (in MPa) at the contact. $e_0$ is assumed to describe the residual aperture of the flow channel under 0 MPa normal stress, while $e_{\inf }$ represents the aperture value under infinite pressure. The value of $e_0$ and $e_{\inf }$ can be recalculated according to Equation (10) based on the laboratory test results of the permeability $k$ of the rock specimen, the length L of the flow channel, and the total volume V of the reservoir rock.

$$e=e_{\inf }+(e_0-e_{\inf })\exp (-0.15{\sigma }_n)$$

(9)

$$k=\frac{1}{12V}{\displaystyle \sum _{\mathrm{pipes}}L{e}^3}$$

(10)

Moreover, in our fluid-mechanically coupled code, the assumption is made that the hydraulic aperture $e_0$ of the flow channel connecting two fluid domains will increase fivefold if the tensile or shear failure takes place between two fluid domains.

3. Model Setup and Micro-Parameter Calibration

3.1. Model Setup

Zhou et al. [12] performed a series of experiments to investigate the influence of different horizontal stresses, angles of approach, and shear strengths of pre-fractures on the interaction between driven fractures and pre-existing natural fractures. In Zhou et al. [12], the mixture of Chinese cement no. 325 and fine sand was used to prepare the test samples, and three types of paper (Type I, II, and III) with different frictional coefficients presented the pre-existing natural fractures. Young’s modulus, Poisson’s ratio, and the uniaxial compressive strength of the cement mortar sample were 8.4 GPa, 0.23 MPa, and 28.34 MPa, respectively. The frictional coefficients of the three types of paper were 0.38 (Type I), 0.89 (Type II), and 1.21 (Type III), respectively, and the cohesion was always 3.2 MPa. In this section, we generate the bonded particle models according to the mechanical properties of test samples in the study of Zhou et al. [12]. The SJM is used to mimic the pre-existing natural fractures with different shear strengths.

Figure 5 presents the cement mortar models with two parallel natural fractures at dip angles of 30°, 60°, and 90°. Each model is assembled by 22,600 circular particles. The width and height are 8000 mm and 6000 mm, respectively. The bonded circular particle radii range from 4 mm to 6 mm following a normal distribution. A borehole with a diameter of 200 mm is created at the center of the model for hydraulic fracturing fluid injection. The parallel natural fractures are emulated by the SJM, for which the length is 2000 mm. According to the servo mechanism, the four walls surrounding the model can move to apply constant confining pressure to the rock specimen, σ₁ in the horizontal direction, and σ₃ in vertical direction.

3.2. Micro-Parameters Calibration

The parameters used in PFC for simulating cement mortar samples are microscale parameters that cannot be directly acquired from laboratory experimental tests. However, they can be optimized based on a series of uniaxial compressive tests (UCT) and direct tensile tests (DTT). According to our experience, the elastic modulus of the model is related to the stiffness of the contact, and Poisson’s ratio is decided by the ratio of normal stiffness to shear stiffness at the contact. The microscale tensile strength, cohesion, and friction angle influence the strength of the sample. By trial and error, the microscale parameters as listed in

Table 2. Input microscale parameters for the modeling of cement mortar and natural fractures.

Input Micro-Parameters
Lower bound of particle radius (mm)	Rmin	4.0
Ratio of particle radius	Rmin/Rmax	2.6
Particle density (kg/m3)	ρ	2350.0
Young’s modulus of the particle (GPa)	Ec	6.7
Ratio of normal to shear stiffness of the particle	kn/ks	2.6
Friction coefficient of particle	μ	0.50
Young’s modulus of the parallel bond (GPa)	${\overline{E}}_c$	6.7
Ratio of normal to shear stiffness of the parallel bond	${\overline{k}}_n/{\overline{k}}_s$	2.6
Tensile strength of the parallel bond (MPa)	${\overline{\sigma }}_{\mathrm{t}}$	8.0
Shear strength of the parallel bond (MPa)	$\overline{\tau }$	18.0
Radius multiplier	$\overline{\lambda }$	1.0
Moment contribution factor	$\beta$	0.2
Normal stiffness of SJM (GPa/m)	kn	1300.0
Shear stiffness of SJM (GPa/m)	ks	500.0
Tensile strength of SJM (MPa)	${\overline{\sigma }}_{\mathrm{t}}{}^{\prime }$	3.2
Friction angle of SJM (°)	${\varphi }^{\prime }$	20.81/41.67/50.42
Cohesion of SJM (MPa)	${\overline{c}}^{\prime }$	3.2

were arrived at, which can be used to reproduce the macroscale mechanical behavior of the test samples. Figure 6 shows the micro-crack distribution and stress-strain curves of the numerical models in the uniaxial compressive test and direct tensile test from which the macroscale mechanical properties were obtained, as listed in

Table 1. Macroscale properties of cement mortar.

Cement Mortar Parameters	Laboratory Values [12]	Modeling Values
Uniaxial compressive strength (MPa)	28.34	28.70
Tensile strength (MPa)	-	3.73
Young’s modulus (GPa)	8.40	8.34
Poisson’s ratio	0.23	0.23
Permeability (mD)	0.1	0.1

. The values of the macroscale properties of the cement mortar agree well with the laboratory values. The parameters of the SJM are also listed in Table 2, representing the three types of paper, which mimic pre-existing natural fractures with different shear strengths.

4. Modeling Results of Interactions between Hydraulic Fractures and Natural Fractures

4.1. Comparison with Laboratory Experimental Results

In the experimental works of Zhou et al. [12], three interaction types have been observed under tri-axial stress states, including Crossed, Dilated, and Arrested type. The characteristics of the three interaction types in laboratory experiments are summarized as follows.

Crossed (Experiment): A natural fracture is crossed by a hydraulic fracture.

Dilated (Experiment): A natural fracture is arrested by an opening and dilating natural fracture.

Arrested (Experiment): A natural fracture is arrested by a shear slippage natural fracture with no dilation and by fluid flow in the natural fracture.

A series of hydraulic fracturing models were performed in two-dimensional stress states to investigate the interactions between hydraulic fractures and natural fractures, considering that the viscosity of the fracturing fluid is 1.0 × 10⁻³ Pa·s and the rate of injection is 2.0 × 10⁻⁴ m³/s/m. Also, in our modeling (

Table 3. Summary of numerical and experimental conditions and the results for the hydraulic fracturing test using the Type I pre-existing interface.

Approaching Angle	σ1 (MPa)	σ3 (MPa)	Experimental Results (Types)	Modeling Results
				Types	100,000 Steps	300,000 Steps
90°	8	3	Crossed	Crossed
90°	8	5	Crossed	Crossed
60°	8	3	Dilated	Dilated-Arrested
30°	8	3	Dilated	Arrested

,

Table 4. Summary of numerical and experimental conditions and the results for hydraulic fracturing test using the Type II pre-existing interface.

Approaching Angle	σ1 (MPa)	σ3 (MPa)	Experimental Results (Types)	Modeling Results
				Types	100,000 Steps	300,000 Steps
90°	10	3	Crossed	Crossed
60°	8	5	Dilated	Dilated
30°	10	5	Dilated	Dilated-Arrested
30°	13	3	Arrested	Arrested

and

Table 5. The summary of numerical and experimental conditions and the results for hydraulic fracturing test using Type III pre-existing interface.

Approaching Angle	σ1 (MPa)	σ3 (MPa)	Experimental Results (Types)	Modeling Results
				Types	100,000 Steps	300,000 Steps
90°	8	3	Dilated	Crossed
60°	13	3	Dilated	Crossed
60°	10	3	Dilated	Crossed
30°	13	5	Dilated	Dilated

), there are three basic scenarios of interaction between hydraulic fractures and natural fractures, which are similar to those in laboratory experiments. The details of the interaction processes are described in the following.

Crossed (Modeling): A hydraulic fracture initiates and propagates along the direction of maximum principle stress, crosses the pre-existing natural fracture at the interaction with no or little fracturing fluid leaking off into it, and continues to propagate with no diversion from the original propagation direction.

Arrested (Modeling): A hydraulic fracture initiates and propagates along the direction of maximum principle stress; as the hydraulic fracture approaching the natural fracture, the natural fracture slips due to the deformation caused by the near-tip stress field of the hydraulic fracture, and a lack of fluid flow in the damaged natural fracture; when the hydraulic fracture encounters the natural fracture, the fracturing fluid leaks off into the damaged natural fracture quickly; and the hydraulic fracture is arrested by the damaged natural fracture, which becomes a part of the hydraulic fracture. Lastly, the hydraulic fracture restarts at one tip of the natural interface, with its orientation rerouted from the dip direction of the natural fracture to the direction of maximum principle stress.

Dilated (Modeling): Initially, the hydraulic fracture propagates along the direction of maximum principle stress; before the hydraulic fracture encounters the natural fracture, no or limited local shear slippage takes place in the natural fracture; while the hydraulic fracture encounters the natural fracture, the fracturing fluid penetrates into the natural fracture and dilates it, which results in shear slippage of the natural interface by decreasing the effective stress while the pore pressure increases at the intersection; and the hydraulic fracture is arrested by the damaged natural fracture, and it will reinitiate at the tip of the natural fracture or a weak point where the tensile stress is larger than the bond tensile strength.

Table 3, Table 4 and Table 5 present the modeling examples of the interactions between hydraulic fractures and natural fractures compared with the experimental results under different coefficients of friction of the natural fracture. Moreover, the figures, including the fracture patterns and pore-pressure distributions at two different times, are shown in Table 3, Table 4 and Table 5. In the figures, the original parallel natural fractures are shown by green lines. The red short line represents a tensile crack in the cement material, and the black short lines are the damaged part of the natural fractures under shear slippage. The pore-pressure in each fluid domain is presented by a solid blue circle, the size of which is related to the magnitude of the pore-pressure.

Based on the evolution of fracture patterns and pore-pressure distributions, the interaction types in our modeling were also summarized and listed in Table 3, Table 4 and Table 5, and they agreed well with the experimental results, except for the Type III pre-existing interfaces. This exception can be attributed to the fact that the thickness of Type III pre-existing interfaces is significantly larger than the thickness of the two other types, which results in a permeability increase. In our modeling, the pre-existing natural fracture was crossed by the hydraulic fracture at high differential stress and angles of approach of 60° or greater, while the pre-existing natural fracture with low shear strength was arrested at angles of approach of 30°. In addition, an interesting phenomenon has been found, which is that combined forms of the three basic interaction types (Crossed, Arrested, and Dilated type) exist in several simulations.

Thus, according to the modeling results, we can arrive at a basic conclusion that the approach angle, in situ differential stress, and the permeability of natural fractures affect the interactions between hydraulic fractures and natural fractures. The influence will be discussed in following sections.

4.2. The Influence of Approach Angle and In-Situ Differential Stress

From previous studies, it is known that the interaction between hydraulic and natural fractures is a complex process, which is related to lots of influence factors such as the formation rock mechanical properties, the strength of the natural fractures, the approach angle, in situ stress, and the permeability of natural fractures. According to the scenarios described in the preceding section, the influence of the approach angle and in situ stress are further studied in this section, with the intention to provide a basic conclusion for selecting stimulation methods in the hydraulic fracturing process.

A series of hydraulic fracturing models were performed considering in situ differential stress changing from 2 MPa to 8 MPa and approach angles varying from 30° to 90°, while the maximum principle stress was kept at 10 MPa in the horizontal direction. In addition, the cohesion and tensile strength of pre-existing fractures were both set to a negligible value of 0.05 MPa. The stimulation parameters were the same as in Section 4.1.

According to the general mechanism of interaction between hydraulic fractures and natural fractures, the Dilated type and Arrested type are collectively called Arrested type. Figure 7 shows the interaction results related to different in situ differential stresses and approach angles. It can be found that the natural fracture is more favorable for opening and diverting fracturing fluid under low approach angles and low differential stresses, which results in the propagation of a hydraulic fracture along the natural fracture and its re-initiation at a weak point or the tip of the natural fracture. On the contrary, under high approach angles and high differential stresses, the hydraulic fracture tends to cross the pre-existing natural fracture. Moreover, the series of modeling results under different approach angles and in situ stress states indicated very good agreement compared to the analytical results based on Blanton’s criterion [9], except in a few cases (Figure 7).

4.3. The Influence of Natural Fracture Permeability

In previous numerical modeling of and analytical research into the interaction between hydraulic fractures and natural fractures, the permeability of natural fractures was not taken into consideration [9,10,11,12,14]. However, the experimental results in Zhou et al. [12] showed that the influence of the permeability of natural fractures on the interaction between hydraulic fractures and natural fractures cannot be omitted. For example, the thickness of the Type III pre-existing interfaces is significantly larger than that of the two other types, which results in the effect that fracturing fluid more easily leaks off into the natural fracture and causes dilation.

The permeability of Type III pre-existing interfaces was assumed to be 100 mD, which is three orders of magnitude higher than the original value.

Table 6. The summary of numerical and experimental conditions and the results for the hydraulic fracturing test assuming that the permeability of the Type III pre-existing interface is 100 mD.

Approaching Angle	σ1 (MPa)	σ3 (MPa)	Experimental Results (Types)	Modeling Results
				Types	100,000 Steps	300,000 Steps
90°	8	3	Dilated	Dilated-Crossed
60°	13	3	Dilated	Dilated-Crossed
60°	10	3	Dilated	Dilated-Crossed
30°	13	5	Dilated	Dilated

presents the summary of the modeled results of interaction between hydraulic fractures and natural fractures. When the driven hydraulic fracture encountered the natural fracture, the fracturing fluid was subsequently diverted into the natural fracture. Obviously, the effective normal stress acting on two walls of a natural fracture decreased owing to an increase in pore-pressure and a decrease in shear resistance force. If shear resistance is lower than shear force, a microscale shear crack will be initiated. Comparing the fracture patterns and pore-pressure distributions of the figures in Table 5 and Table 6, it can be observed that the interactions between hydraulic fractures and natural fractures also depend on the permeability of natural fractures. At this time, the interaction types of the modeling agree well with the related experimental results at different conditions.

5. Conclusions

A series of simulations for hydraulic fracturing in rock was performed in PFC^2D by our own fluid-mechanically coupled model to investigate the interactions between hydraulic fractures and natural fractures compared against a laboratory experimental study. The modeling results agree well with previous laboratory experimental works, and the main results are summarized as follows:

(1)

In this study, The Smooth Joint Model (SJM) shows an ability to emulate the elastic deformation, shear failure, and tensile failure of rock joints, which is firstly used to simulate the details of interactions between hydraulic fractures and natural fractures in PFC^2D.

(2)

The modeling results show that there are three basic scenarios of interaction between hydraulic fractures and natural fractures, namely, Crossed type, Arrested type, and Dilated type, which are similar to those in laboratory experiments by Zhou et al. Overall, the modeling examples agree well with the experimental results.

(3)

Moreover, the natural fracture is more favorable for opening and diverting fracturing fluid under low approach angles and low differential stresses, which results in the propagation of a hydraulic fracture along the natural fracture and re-initiation at a weak point or the tip of the natural fracture. On the contrary, under high approach angles and high differential stresses, the hydraulic fractures tend to cross pre-existing natural fractures. In addition, a series of modeling results under different approach angles and in situ stress states indicated very good agreement compared with analytical results based on Blanton’s criterion, except for a few cases.

(4)

In natural fractures with higher permeability, the fracturing fluid was diverted into natural fractures more easily. Undoubtedly, the effective normal stress acting on two walls of the natural fracture decreased owing to an increase in pore-pressure and a decrease in shear resistance force, which results in shear slippage when the shear resistance is lower than the shear force. The permeability of natural fractures will affect the interaction mechanism between hydraulic fractures and natural fractures, which results in a variety of the fracture patterns and pore-pressure distributions in naturally fractured reservoirs.

To summarize, the propagation of hydraulic fractures in naturally fractured reservoirs is a complex process. In this study, the SJM embedded in PFC^2D was successfully applied to simulate different scenarios of the interactions between hydraulic fractures and natural fractures compared with the experimental results. This method is a powerful tool for investigating the hydraulic fracturing behavior and optimizing fracture design in naturally fractured shale gas reservoirs.