A Numerical Study on the Diversion Mechanisms of Fracture Networks in Tight Reservoirs with Frictional Natural Fractures

An opened natural fracture (NF) intercepted by a pressurized hydro-fracture (HF) will be diverted in a new direction at the tips of the original NF and subsequently form a complex fracture network. However, a clear understanding of the diversion behavior of fracture networks in tight reservoirs with frictional NFs is lacking. By means of the extended finite element method(XFEM), this study investigates the diversion mechanisms of an opened NF intersected by an HF in naturally fractured reservoirs. The factors affecting the diversion behavior are intensively analyzed, such as the location of the NF, the horizontal principal stress difference, the intersection angle between HF and NF, and the viscosity of the fracturing fluid. The results show that for a constant length of NF (7 m): (1) the upper length of the diverted fracture (DF) decreases by about 2 m with a 2 m increment of the upper length of NF ( L u p p e r ), while the length of DF increases 9.06 m with the fluid viscosity increased by 99 mPa · s; (2) the deflection angle in the upper parts increases by 30.8° with the stress difference increased by 5 MPa, while the deflection angle increases by 61.2° with the intersection angle decreased by 30°. It is easier for the opened NF in lower parts than that in upper parts to be diverted away from its original direction. It finally diverts back to the preferred fracture plane (PFP) direction. The diversion mechanisms of the fracture network are the results of the combined action of all factors. This will provide new insight into the mechanisms of fracture network generation in tight reservoirs with NFs.


Introduction
With the technological progress in petroleum industries, petroleum engineers are increasingly concerned with the exploration and development of tight reservoirs in recent years. Due to the ultra-low matrix permeability, hydraulic fracturing is a key technology for enhancing the recovery of tight hydrocarbon reservoirs [1][2][3][4][5][6][7][8][9][10][11]. Activation of preexisting natural fractures (NFs) during fracturing

Governing Equations of Hydraulic Fracturing Problems
As shown in Figure 1, the domain Ω denotes a tight reservoir, which includes an HF and an NF. The injection point is located on the middle left of the domain Ω, and the corresponding pump rate is denoted by Q 0 . As is known, the hydraulic fracturing problem is essentially a fluid-solid interaction process, so its governing equations consist of two parts: the stress equilibrium equation for rock skeleton and fluid pressure equation in the hydraulically driven fracture [52][53][54]. (1) The Stress Equilibrium Equation: According to the theory of elasticity, the stress equilibrium equation is expressed: where σ denotes the stress tensor; b denotes the body force vector in the rock skeleton. As shown in Figure 1, the boundary conditions are composed of a displacement boundary condition (Γ u ) and a force boundary condition (Γ t ). They are expressed as u =ū onΓ u , σ · n = t, onΓ t , σ · n HF = pn HF , onΓ HF , σ · n NF = t NF , onΓ NF , where p denotes the fluid pressure on the artificial fracture; t NF denotes the contact traction vector on the NF surface Γ NF ;ū denotes the displacement imposed on the boundary Γ u ; t denotes the traction vector imposed on the boundary Γ t .
For brittle rocks, there is a linear relationship between stress tensor and strain tensor under a small deformation assumption, so the corresponding constitutive equation is expressed as where D denotes the fourth elasticity tensor; ε denotes the strain tensor; the symbol ":" denotes the double dot product of the two tensors.
Under the assumption of small deformation, the relationship between displacement vector u and strain tensor ε are as follows: (2) Fluid Pressure in an HF: Under lubrication theory assumptions, the velocity profile of the fluid in the HF is that of a planar Poiseuille flow between two parallel plates. Therefore, the fluid pressure in the HF can be expressed as [15,33,34,55,56] ∂w ∂t where w denotes the fracture opening of HF; s denotes the crack propagation direction; t denotes the injection time; k denotes the fracture transmissivity. According to the cubic law, the fracture transmissivity can be expressed as where µ denotes the viscosity of fracturing fluid. The corresponding initial and boundary conditions can be expressed as where s tip denotes the tips of HFs; q denotes the injection rate of fracturing fluid at the crack point s and time t. It is can be seen that Equation (7) satisfies mathematically the Neumann boundary condition. In order to get a unique solution for the fluid pressure equation, the constraint condition should be additionally imposed. The necessary condition, i.e., the conservation of global mass in the HF can be written as

Crack Propagation Criterion
According to the theory of fracture mechanics, the maximum circumferential stress criterion is adopted to determine the propagation direction of a hydraulically driven fracture at every time t. The artificial fracture will propagate along a direction perpendicular to the maximum circumferential stress. If the stress intensify factor K is no less than the fracture toughness of rock skeleton K IC , the crack will propagate along a certain direction. The interaction integral method in domain form is utilized to calculate the stress intensity factors K I and K II . The following equation of the interaction integral can be written as [57] where I (1,2) denotes the interaction integral; W (1,2) denotes the interaction strain energy as follows; q w (x) denotes the smooth weighting function, which takes a value from 0 to 1; δ denotes the Kronecker symbol; the superscripts (1) and (2), respectively, denote the current state and the auxiliary state for the stress and strain field. The corresponding calculation procedures can be described in detail in [40,57,58] W (1,2) = σ (1) ij ε ij ε (1) ij .
The direction of crack propagation θ can be computed in a local tip coordinate system: where the symbol "arctan" denotes the arc-tangent function.

The Cross Criterion between HF and Frictional NF
As is known, when an HF encounters a frictional NF, there are three possible scenarios: arrested, direction-crossing, or a crossing with an offset [16,23,41]. Here, the extended Renshaw and Pollard rule is adopted to determine the interaction behavior between HF and NF. As shown in Figure 2, if a new fracture initiates on the opposite side of the NF, the maximum principle stress σ 1 will reach the rock tensile stress. Meanwhile, a no-slip condition should be satisfied along the NF surface. Otherwise, the HF will cross directly or branch into NF with an offset. The expressions of the combined shear stress and normal stress are shown in Equation (12), which is described in detail in other references [15,23].
where β denotes the intersection angle between HF and NF; τ β denotes the combined shear stress on the NF surface under the action of remote stress and the local crack tip stress; σ βy denotes the combined normal stress; T 0 denotes the rock tensile strength; S 0 denotes the cohesion force of the frictional NF; µ f denotes the frictional coefficient of the NF surface.

XFEM and Discretization of the Governing Equations of the Hydraulic Fracturing Problem
The XFEM (extended finite element method) is utilized to approximate the displacement discontinuity on both sides of the HF. In order to represent the multiple cracks, a Junction enrichment function is introduced, as shown in Figure 3. The enriched displacement field can be written as [59].
where N, N dis , N tip and N jun , respectively, denote the set of standard nodes, Heaviside enrichment nodes, crack-tip nodes, and junction enrichment nodes; u denotes the standard nodal d.  The Heaviside enrichment function is expressed as [59] H The crack-tip enrichment function is defined as where (r, θ) denotes the local crack-tip coordinate in the polar coordinate system. The junction enrichment function J H (x) is defined as where ϕ m (x) and ϕ s (x), respectively, denote the signed distance function of the main crack and the secondary crack. It can be seen that J H (x) is equal to 1, −1, or 0 on different sub-domains divided by the secondary cracks. According to the finite element method, the pressure field p and the fracturing opening displacement vector w can be respectively approximated as where U and P, respectively, denote the global nodal displacement vector and the nodal pressure vector; N p (s) and N w (s), respectively, denote the matrix of the shape function of the fracture opening and pressure.
By substituting the above XFEM formulation, displacement and pressure approximations into the weak form of stress equilibrium equation and lubrication equation, the corresponding discretization forms are written as where K denotes the global stiffness matrix; Q denotes the coupling matrix; F ext denotes the external loading vector; H denotes the flow matrix; ∆t denotes the time step; and S denotes the source term. They are, respectively, defined as follows [33,34]: In the above Equation (21), D cont denotes the contact stiffness matrix of fracture interfaces: As shown in Figure 4, the flow rate in the main crack and secondary crack satisfies the law of conservation, i.e., Q 0 = Q 1 + Q 2 , where Q 1 and Q 2 denote the flow rate in Branches 1 and 2 of the secondary crack, respectively. The nonlinear fluid-solid coupling system of equations of hydraulic fracturing problems can be numerically solved by the Newton-Raphson method. More details are described in [33,34].

Verification of the XFEM Model
For the model verification, the results from our models is summarized in Appendix A [33,34]. The verification model of this XFEM code is described in detail in [33,34], so the related process is not repeated in this article. It is shown that the numerical results have good agreement with the experimental results of true tri-axial hydraulic fracturing by TerraTek, Inc. (Salt Lake City, UT, USA).
For further details of numerical and experimental procedures and the corresponding results, we refer the reader to References [33,34,60,61].

Effect of the Location of Natural Fractures on the Diversion of Fracture Network Propagation
In this section, the effect of the location position of an NF on HF propagation paths is determined by numerical simulation using XFEM. The input parameter values of this model are as shown in Table 1. Under isotropic stress state conditions, the intersection angle between HF and NF is equal to 90 • . As shown in Figure 1, the domain is a 25 m × 25 m square, where the injection point is located at the midpoint of the left edge. In this domain, the HF is 2.6 m in length, and the length of NF is equal to 7 m. Based on the above input parameters, the mechanical NF-HF interaction processes in hydraulic fracturing are numerically simulated at different lengths of NF in lower and upper parts, i.e., corresponding to L lower and L upper in Figure 1, respectively. The corresponding crack propagation paths are as shown in Figure 5. It is obvious that fracture diversion occurs near the tips of the NF in all cases. In Figure 5a, when the HF intersects with the NF, the fracturing fluid flows into the opened NF. In the lower parts of the NF, the opened NF firstly propagates along a vertically downward path for a certain length, and it is then diverted along a new direction; in the upper parts of the NF, the opened NF is directly diverted at the upper tip of the original NF. In Figure 5b, both the lower and upper parts of the NF firstly extend vertically downward and upward for a short distance, respectively, and are then diverted to the right-hand side of the graph. However, the length of the lower parts of the diverted fracture (DF) is longer than that of the upper parts of the DF, corresponding to the red line in Figure 1. In Figure 5c, in the upper parts of the NF, the opened fracture can only propagate vertically upward, and cannot be diverted near the tip of the NF; in the lower parts of the NF, the opened fracture can be diverted away from the original NF. By making use of the data in Figure 5, the length of the DF propagation in the upper parts is calculated from the upper tip of the original NF. When the upper length of NF, i.e., L upper , is equal to 4, 5, and 6 m, the corresponding length is 3.14, 2.15, and 1.14 m, respectively. If L upper is increased by 2 m, the upper length of the DF will decrease by 2 m. Therefore, it is shown that the longer the upper parts of the original NF are, the more difficult it is for the opened NF to be diverted away from the upper tip of the NF under the conditions of an isotropic stress state, while the lower parts of the original NF is more easily diverted to the right-hand side than the upper parts of the original NF under this circumstance.  The Von-Mises stress distributions are shown in Figure 6, where a blue color represent a relative stress value, while a yellow or red one represent a higher stress value. For all cases of the model, it is shown that there is a small region of stress concentration near the two tips of the original NF, which indicates that a higher pressure is required to divert the opened NF away from the original NF. The fracture aperture and net pressure curves of the diverted fracture are shown in Figures 7 and 8, respectively. With the decrease of L lower , both curves have revealed an asymmetrical characteristic, which peak at the diverted point in the lower parts. In addition, at the intersection point between HF and NF, their corresponding values take second place. The fracture aperture and net pressure in the lower parts of the DF are much greater than those in the upper parts. By comparison, in the case of L lower = 3 m and L lower = 4 m, their curves are nearly symmetrical. This indicates that, under the combined action of remote stress and local crack-tip stress, the variation tendency of the fracture aperture and the net pressure is quite different from that in the case of only a single HF.  The flow rate in the diverted fracture is shown in Figure 9. When fluid flows into the intersection point between HF and NF, it will flow upward and downward, respectively. If the value of the flow rate is negative, fluid will flow upward; otherwise, it will flow downward. It is shown that the flow rate in the lower parts of the DF is much greater than that in the upper parts. Therefore, the fracture aperture and net pressure in the lower parts of the DF is greater than that in the upper parts under the condition of the same fluid viscosity. With the decrease of L lower , fluid flows downward more easily. Thus, the lower parts of the original NF is more easily diverted than the upper parts of the original NF. This may explain the results in Figure 5.

Effect of Horizontal Stress Differences on the Diversion of Fracture Network Propagation
Based on input parameters in Table 1, the effect of the remote stress difference on the diversion of fracture network propagation is numerically simulated at different levels of minimum horizontal stresses σ h : 5, 4, and 0 MPa; the maximum horizontal stress σ H is kept constant (5 MPa) in all cases of this model. At the same time, L lower = 3 m, L upper = 4 m, and the NF-HF intersection angle is equal to 90 • .
The corresponding crack propagation paths are shown in Figure 10, in which the deflection angle is defined the angle between NF and DF at the tips of the original NF. According to the data in Figure 10, the deflection angle in the upper parts is calculated. When the horizontal stress difference is, respectively, equal to 0 and 5 MPa, the corresponding deflection angle is 59.2 • and 90 • , respectively. If the stress difference is increased by 5 MPa, the deflection angle in the upper parts is increased by 30.8 • . The higher the horizontal stress difference is, the greater the deflection angle is. In Figure 10c, i.e., ∆σ = 5 MPa, the opened NF firstly diverts and propagates along the direction of minimum horizontal stress and finally tends to extend along the preferred fracture plane (PFP) direction in petroleum engineering. By contrast, when the stress state is approximately isotropic, crack propagation in both the lower and upper parts will extend along the minimum horizontal stress direction for some length. This indicates that crack propagation of the opened NF is a complex mechanical process under the combined action of the local crack-tip and the remote stress state.
As shown in Figure 11c, there is a lower Von-Mises stress region (corresponding to blue area on the contour) located on the right of DF for 5 MPa stress difference. This propagates the diverted fracture along the PFP. In Figure 11a,b, both cases correspond to a lower stress difference, and the lower stress region is mainly near the two crack tips. The local stress distribution is an explanation to interpret crack paths in Figure 10.
The fracture aperture and net pressure curves of DF are shown in Figures 12 and 13, respectively. The fracture aperture curve reveals an asymmetrical characteristic, while the net pressure curve reveals a nearly symmetrical characteristic. The fracture aperture peaks at a global maximum value at the inflection point, where the opened fracture in the upper parts diverts to the right side in Figure 10; at the NF-HF intersection point, the fracture aperture takes the second place; at the inflection point in the lower parts, it takes the third place. The net pressure reaches a maximum at the intersection point, and there are two inflection points on this curve. This means that the opened fracture diverts to the right side. The higher the stress difference is, the greater the net pressure it requires, which leads to a greater fracture aperture.  The flow rate in the DF is shown in Figure 14. It is obvious that the flow rate in the upper parts is much greater than that in the lower parts. Therefore, the fracture aperture and net pressure in the upper parts are greater than those in the lower parts for the same fluid viscosity. This might explain the results of fracture aperture and net pressure in Figures 12 and 13.

Effect of the NF-HF Intersection Angle on the Diversion of Fracture Network Propagation
Based on the input parameters in Table 1, the effect of the NF-HF intersection angle on the diversion of fracture network propagation is numerically simulated at different levels of intersection angles β: 75 • , 60 • , and 45 • ; both the maximum and minimum horizontal principle stresses are equal to 5 MPa in all cases of this model. Meanwhile, L lower and L upper are, respectively, 3 m and 4 m.
The corresponding crack propagation paths are as shown in Figure 15. Under the condition of isotropic stress state, the NF-HF intersection angle will have a significant impact on the propagation direction of the primary HF, i.e., the black dotted line in Figure 15, when HF is approaching NF. With the decrease of the intersection angle, the primary HF deflects from the horizontal line. When the NF-HF intersection angle is greater than 60 • (in Figure 15a,b), the opened NF in the upper parts is more easily diverted away from the original NF than that in the lower parts under the combined action of remote stress and the crack-tip stress field. However, the intersection angle is less than 60 • (in Figure 15c). The opened NF in the lower parts is more easily diverted away from the original NF than that in the upper parts under the combined action of remote stress and the crack-tip stress field. By making use of the data in Figure 15, the deflection angle for the primary HF was calculated. When the intersection angle is decreased from 75 • to 45 • , the corresponding deflection angle is increased from 0 • to 61.2 • . This indicates that the NF-HF intersection angle will have a significant impact on the diversion propagation of the primary HF and the secondary opened NFs.
The Von-Mises stress distributions at different levels of NF-HF intersection angle are shown in Figure 16. In Figure 16a, there is a stress concentration region near the diversion point in the upper parts of the NF, which indicates that it will require a high net pressure to divert the opened fracture upward. In Figure 16b, the Von-Mises stress in the upper parts is greater than that in the lower parts, which causes the subsequent fracture to easily divert upward. In Figure 16c, the Von-Mises stress on the left of the NF is greater than that on the right, so it is more easily diverted downward.
The fracture aperture and net pressure curves of the DF are shown in Figures 17 and 18, respectively. In the case of β = 75 • , the fracture aperture and net pressure in the upper parts is greater than that in the lower parts. In particular, the fracture aperture and net pressure near the diversion point in the lower parts is close to zero, and this is consistent with the results of Von-Mises stress at this point. In the other two cases, the maximum values of the fracture aperture and the net pressure are at the diversion point in the lower parts. The smaller the intersection angle is, the greater the net pressure it requires to divert the fracture.
The flow rate in DF is shown in Figure 19. It is obvious that, in the case of β = 45 • , the flow rate in the upper parts is much greater than that in the lower parts. This indicates that it is easy for fracturing fluid to flow upward when the intersection angle is small. This is a possible reason for explaining the results in Figures 17 and 18.

Effect of Fluid Viscosity on the Diversion of Fracture Network Propagation
Based on the input parameters in Table 1, the effect of the viscosity of fracturing fluid on the diversion of fracture network propagation is numerically simulated at different levels of viscosity µ : 100 mPa·s, 10 mPa·s, and 1 mPa·s [62]. In this model, the maximum and minimum horizontal stresses are kept constant (5 MPa) for all cases. Meanwhile, the NF-HF intersection angle is equal to 90 • .
The corresponding crack propagation paths are shown in Figure 20. By making use of the data in Figure 20, the length of DF is calculated. When the fluid viscosity is increased from 100 mPa·s to 1 mPa·s, the length of DF is increased from 14.29 to 5.23 m. It is clear that the smaller the viscosity is, the more difficult the opened fracture diverts into a new direction. When the viscosity is equal to 1 mPa·s, artificial fracture will propagate along the NF direction under given conditions. This indicates that more energy is required to divert the opened NF upward and downward. d i s t a n c e a l o n g t h e f r a c t u r e l e n g t h ( m ) a n g l e = 7 5 o a n g l e = 6 0 o a n g l e = 4 5 o Figure 17. The fracture aperture curves of the DF along the fracture length direction at different levels of intersection angle between HF and NF. The distance in the x-axis is along the direction from the lower part to the upper part of the DF. d i s t a n c e a l o n g t h e f r a c t u r e l e n g t h ( m ) a n g l e = 7 5 o a n g l e = 6 0 o a n g l e = 4 5 o Figure 18. The net pressure curves of the DF along the fracture length direction at different levels of intersection angle between HF and NF. The distance in the x-axis is along the direction from the lower part to the upper part of the DF.  The Von-Mises stress distributions at different levels of fluid viscosity are shown in Figure 21. In Figure 21c, there is a lower Von-Mises stress area on the right of the NF, which makes the opened NF propagate along the original NF direction. This is consistent with the results in Figure 20.
The fracture aperture and net pressure curves of the DF are shown in Figures 22 and 23, respectively. They are close to symmetrical about the axis of the original HF under the condition of an isotropic stress state. In the cases of µ = 100 mPa ·s and µ = 10 mPa ·s, there are two inflection points on the curves, which correspond to the diversion point in the lower and upper parts of the NF. The greater the fluid viscosity is, the greater the fracture aperture and net pressure are.
The flow rate in the DF is as shown in Figure 24. Obviously, the greater the viscosity is, the greater the flow rate is, and thus the easier it is for the secondary fracture to divert. This might explain the results in Figures 20 and 21 [63].

Conclusions
This paper investigates the diversion mechanisms of a fracture network in tight formations with frictional NFs by means of the XFEM technique. The effects of some key factors such as the location of the NF, the intersection angle between the NF and HF, the horizontal stress difference, and the fluid viscosity on the mechanical diversion behavior of the HF were analyzed in detail. The following main conclusions can be drawn: (1) Fracture diversion propagation will occur near the two tips of the opened NF after an HF is intersecting with an NF. The numerical results show that some key factors such as the NF position, the NF-HF intersection angle, the horizontal stress differences, and the fluid viscosity have a significant impact on the diversion propagation in the upper and lower parts of the opened NF.

Conflicts of Interest:
The authors declare no conflicts of interest.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A
The result of XEFM is here compared with results of analytical solutions. As is known, depending on the dimensionless fracture toughness, the analytical solutions of the Kristianovich-Geertsma-de Klerk (KGD) model have different expressions. The dimensionless fracture toughness can be written as where K m denotes the dimensional fracture toughness; K IC denotes the rock fracture toughness; µ denotes the viscosity of fracturing fluid; Q 0 denotes the injection rate; E denotes the rock Young's modulus; ν denotes the Poisson's ratio of the rock matrix. If K m is greater than 4, the fracture propagation regime is toughness dominated; if K m is less than 1, the fracture propagation regime is viscosity dominated, which is much more common in most hydraulic fracturing treatments. The input parameters of the verification model are listed in Table A1. According to Equation (A1), the dimensionless fracture toughness is equal to 0.313, which indicates that the fracture propagation regime is viscosity-dominated. In this model, the HF is located at the center of a symmetrical model with a length of 100 m and 180 m along the xand y-direction directions, respectively. The domain is divided into 3080 bilinear quadrilateral elements. The initial half-length of the HF is equal to 1.25 m, and it is assumed that a constant fluid pressure acting on the fracture wall is equal to 3.9 MPa. The curves of fluid pressure at the injection point and the fracture width at 30 s are shown in Figure A1a,b, respectively, which is compared with the corresponding analytical solutions. There is very good agreement between the numerical results and analytical solutions, which indicates that the XFEM model can obtain reliable results.