A 2D explicit numerical scheme–based pore pressure cohesive zone model for simulating hydraulic fracture propagation in naturally fractured formation

Abstract

In this work, a 2D pore pressure cohesive zone model is presented to simulate the hydraulic fracture propagation in naturally fractured formation, in which the fracturing process is governed by a bilinear cohesive zone model equipped with Coulomb's friction law and the fluid flow within the hydraulic fracture is described by the lubrication equation. In this model, fluid leak‐off into rock matrix is ignored and a fully explicit temporal integration scheme is adopted to overcome the convergence issue of conventional implicit scheme (eg, Newton‐Raphson method). The advantage of the proposed model is that it does not require any special crossing criterion to determine the interaction behavior when the hydraulic fracture hits the natural fracture. Implementation of the model was described in detail, and then, the model was verified with well‐known analytical solution of KGD problem and criteria of hydraulic fracture crossing natural fracture. The capability of the proposed model to capture the interaction behavior between hydraulic fracture and natural fracture was demonstrated by the good agreement of the modeling result and analytical solution. Several numerical cases were performed to investigate the impact of key factors on fracture network evolution during hydraulic fracturing treatment. Results show that fracture network is not only governed by the spatial distribution of natural fracture, but also greatly affected by the initial hydraulic aperture of natural fracture and in situ stress contrast.

INTRODUCTION

Hydraulic fracturing is a widely used technique to create a highly permeable zone by the injection of a pressurized fluid into a formation and serves as the primary approach to enhance the recovery of hydrocarbon resources from low‐permeability reservoirs. A number of methods, such as laboratory experiment and microseismic mapping technology, have been developed to optimize the hydraulic fracturing treatment design. However, laboratory experiment is difficult to provide direct guide to field‐scale hydraulic fracturing, although it is a valuable way to understand the mechanism of hydraulic fracture initiation and propagation. In addition, because of the economic constraints, not all hydraulic fracturing treatments are monitored with microseismic mapping technology. Mathematical model is an alternative and, perhaps, more effective approach for optimizing fracturing treatment design and thus enhancing outcomes and reducing costs. Actually, mathematical models have been playing an increasingly crucial role in understanding and predicting the hydraulic fracture propagation for decades.

The history and recent trends of the development of hydraulic fracturing model can be found in Adachi et al and Lecampion et al. The earliest analytical 2D and 3D models are KGD, PKN, and penny‐shaped, which have now been largely replaced by pseudo‐3D models. In addition, some advanced models have been implemented in various numerical frameworks such as finite element method (FEM), boundary element method (BEM), and discrete element method (DEM). In the framework of FEM, adaptive remeshing technique, cohesive element, and the partition of unity methods (eg, XFEM and GFEM) are the most widely used tools for simulating hydraulic fracturing. Displacement discontinuity method (DDM), a special boundary element method which was first proposed by Crouch et al, is particularly suitable for hydraulic fracture propagation modeling due to its relatively computational efficiency as well as simplicity of meshing and remeshing. For example, Weng et al developed unconventional fracture model based on DDM and investigated the complex‐fracture‐network propagation in a formation with preexisting natural fractures. Considering the effect of perforation erosion, Long and Xu formulated a DDM approach to model the simultaneous propagation of multiple hydraulic fractures. Though the above models work quite well on the simulation of fracture network evolution during hydraulic fracturing treatment, the intersection behavior between HF and NF needs to depend on a predefined crossing criterion in the models. However, the crossing criterion may begin to lose effectiveness under complex fracture network condition, because the local stress field varies due to the opening or shear slip of surrounding fractures. DEM can naturally incorporate fractures into its framework, therefore being a potential way to study the interaction between HF and NF. For example, Zhang et al combined PFC and FLAC to study the behavior of a hydraulic fracture crossing natural fractures. Fatahi et al investigated the interaction mechanisms of HF and NF by using PFC2D. DEM may be useful for understanding the mechanisms of intersection behavior. However, it is difficult to link the microscopic properties of particles to macroscopic properties of rock. The literature review indicates that although extensive studies exist on the propagation of hydraulic fracture, it is still not easy to capture the behavior of hydraulic fracture propagation in naturally fractured reservoir.

Cohesive zone model (CZM), which originated from Dugdale and Barenblatt, provides an alternative method for characterizing the fracturing process. A fracture process zone characterized by a traction‐separation constitution law is assumed at the fracture tip, which avoids the stress singularity at the tip in linear elastic fracture mechanics (LEFM). CZMs have been widely used to analyze the fracture problem in various solids since it was first proposed and a variety of CZM constitutive laws are summarized in Park and Paulino. Although CZMs have been playing an increasingly influential role in the simulation of hydraulic fracture propagation since the pore pressure cohesive zone (PPCZ) element was provided in ABAQUS, an obvious limitation is that the PPCZ element can only solve hydraulic fracturing problems with very simple geometry. In addition, almost all of the existing PPCZ elements use implicit integration scheme, which could encounter convergence difficulties when solving problems with high nonlinearity and/or complex contact interactions. Unfortunately, the process of hydraulic fracture propagation in naturally fractured formation is a typical problem characterized by these two features.

The goal of this paper is to present a fully coupled CZM based on explicit integration scheme to capture the evolution of fracture network during hydraulic fracturing treatment in naturally fractured reservoir. The structure of this paper is as follows. Section 2 describes the numerical simulation methodology, including governing equations and finite element implementation. Next, the capability of the proposed model is verified in Section 3. The impact of spatial distribution and initial hydraulic aperture of natural fractures as well as the stress contrast on the interaction behavior is studied in Section 4. A summary and some implications based on numerical results are offered in Section 5.

SIMULATION METHODOLOGY

Modeling hydraulic fracture propagation in naturally fractured formation is of challenge, in which at least three components should be coupled simultaneously: (a) fracture initiation and propagation; (b) fluid flow within the fracture; and (c) rock deformation induced by the fluid pressure on the fracture surfaces. In addition, another important aspect of hydraulic fracture propagation in naturally fractured formation is the interaction between propagating fractures and preexisting geological discontinuities, such as natural fractures, joints, and bedding planes. Although a large number of numerical models have been developed for modeling the process of hydraulic fracturing in recent years, it is still a tough task to capture the behavior of fracture network evolution in naturally fractured reservoir. To this end, an explicit temporal integration–based pore pressure cohesive zone (Explicit‐PPCZ) model is developed. As illustrated in Figure A, the rock domain is discretized by triangular elements, and zero‐thickness 6‐node cohesive elements are inserted between triangular elements to characterize the fracture growth and fluid flow with fractures. At the intersection, the associated cohesive elements are connected by sharing the same pore pressure node, as shown in Figure B. By utilizing this meshing scheme, the difficulties to handle the intersection between hydraulic fracture and natural fracture can be avoided. In this section, the governing equations, weak forms, finite element discretization, and implementation of the Explicit‐PPCZ method will be provided.

Mesh discretization for fracture network (A) and the connectivity of CZ elements at intersection point (B)

Governing equations

Consider a two‐dimensional domain Ω containing a hydraulic fracture Γ_HF driven by the injection of Newtonian fluid with an injection rate of Q_inj, as shown in Figure . Domain Ω includes natural fractures Γ_NF. The displacements of boundary Γ_u are fixed, and stress t is imposed on the boundary Γ_t. In addition, assuming that the permeability of rock is extremely low and thus it is reasonable to neglect the poroelastic effect of rock matrix and the fluid leak‐off from the fracture surfaces. The details of the governing equations to describe the deformation of rock matrix and the fluid flow within fractures are provided in the following sections.

Problem of hydraulic fracture growth and boundary conditions

Deformation of rock

The rock is assumed as an isotropic, homogenous, and elastic medium, and thus, the momentum balance can be expressed as[Image Omitted. See PDF]where σ_ij is the stress tensor, b_i is the body force, ρ_s is rock density, and u_i is displacement field, and α is the damping coefficient.

The rock deformation is described by Hooke's law as[Image Omitted. See PDF]where D_ijst is the elasticity material property matrix of the triangular element and ε_ij is the infinitesimal strain tensor defined as ${\epsilon }_{\mathit{ij}}=\frac{1}{2}\left(u_{i,j}+u_{j,i}\right)$.

Fluid flow within fracture

Fluid flow within a fracture can be described by the fluid mass conservation equation[Image Omitted. See PDF]where w is the local fracture width, ρ_f is the fluid density, q is the fluid flux, t is the elapsed time, and s is the longitudinal coordinate along the fracture. For slightly compressible fluid, the relationship between pressure and density can be written as[Image Omitted. See PDF]where p₀ is the datum pressure, which is assumed to be zero here; ρ₀ is the fluid density at datum pressure; and K_f is the bulk modulus of the fluid.

Fluid flow within smooth parallel fracture yields[Image Omitted. See PDF]where μ is the viscosity of fluid.

The lubrication equation can be obtained by introducing Equations (4) and (5) into (3) as[Image Omitted. See PDF]

Fracture aperture and flow rate naturally yield the following boundary conditions[Image Omitted. See PDF][Image Omitted. See PDF]

in which l is the length of fracture.

Cohesive zone model

Cohesive zone model (CZM) is a phenomenological model introduced by Barenblatt and Dugdale to overcome the limitations of classical linear fracture mechanics such as the stress singularity at the fracture tip. The CZM assumes that there is a fracture process zone ahead of the fracture tip in which the mechanical behavior is governed by a traction‐separation law. Further details on cohesive zone models can be found in Refs. For rocklike materials, friction between fracture surfaces plays a significant role in the shear activation of natural fractures, and therefore, a CZM equipped with frictional contact capability should be utilized to simulate the fracture network in naturally fractured formation. In this work, we consider that the static and dynamic friction coefficients are equal. As shown in Figure , we utilize the following CZM that combined the bilinear traction‐separation law and Coulomb friction law to calculate the normal and shear tractions in the interface element[Image Omitted. See PDF]

Cohesive zone model equipped with Coulomb's friction law: (A) tension; and (B) shear

in which[Image Omitted. See PDF]where $\overline{\sigma }$ and $\overline{\tau }$ are cohesive tensile and cohesion strength of the interface element. ${\delta }_n$ and ${\delta }_t$ are the normal and the shear separation of the interface element at any instant of time. ${\delta }_n^0$ and ${\delta }_t^0$ are the normal and the shear separation at corresponding damage initiation. ${\delta }_n^1$ and ${\delta }_t^1$ are the normal and the shear relative displacements at the complete failure. ${\tau }_f$ is the interfacial friction. In Figure A, the area below the traction‐separation curve is denoted by $G_{\mathit{IC}}$, which represents the fracture energy per unit area of the fracture surface.

A quadratic criterion is implemented to determine the mixed‐mode failure, as follows[Image Omitted. See PDF]where $G_{\mathrm{I}C}$ and $G_{\text{II}C}$ are tensile and shear energy release rates, respectively; $G_{\mathrm{I}}$ and $G_{\text{II}}$ are the energies absorbed in normal and shear directions, respectively; and D is a damage index ranging from 0 (undamaged) to 1 (completely failed). Once D reaches 1, all cohesive tractions (not including the friction component) are zero.

For a linear elastic fracture under plane strain condition, the relation between energy release rate $G$ and fracture toughness $K$ is given as [Image Omitted. See PDF]in which $K_{\mathrm{I}C}$ and $K_{\text{II}C}$ are mode I and mode II fracture toughness, respectively; E is Young's modulus; and $v$ is the Poisson ratio of rock.

Finite element implementation

Weak forms

The weak form of the rock equilibrium equation is derived from the principle of virtual work. By ignoring the body force, the summation of virtual strain energy and virtual kinetic energy of rock domain is equal to the summation of virtual work done by cohesive traction and fluid pressure on the hydrofracture surfaces (Γ_HF) as well as the external traction on boundary (Γ_t)[Image Omitted. See PDF]where $\delta \mathbf{\epsilon }$, $\delta \mathbf{u}$, and $\delta \mathbf{\Delta }$ are virtual strain, virtual displacement, and virtual fracture separation, respectively.

Assume that the datum pressure $p_0$ is zero, and the weak form of the fluid flow within fracture can be obtained by multiplying lubrication Equation (6) by a test function $\delta p$ from left and then integrating the product over the domain Γ_HF[Image Omitted. See PDF]

Using integration by parts and divergence theorem, the right‐hand side of Equation (14) can be recast as[Image Omitted. See PDF]

Substituting Equation (15) into Equation (14), the weak form of the equation governing the fluid flow within fracture can be rewritten as[Image Omitted. See PDF]where $\partial \mathrm{\Gamma }$ denotes the tip of hydraulic fracture.

Discretization in time and space

Due to the convergence issue of implicit algorithms for simulating complex fracture network evolution during stimulation treatment, explicit integration scheme was employed in this work. By introducing triangle elements and PPCZ elements to discretize $\mathrm{\Omega }$ and ${\mathrm{\Gamma }}_{\mathit{HF}}$, respectively, the discrete form of coupled governing Equations (13) and (16) can be given by[Image Omitted. See PDF][Image Omitted. See PDF]where $\mathbf{M}$ is the lumped mass matrix; $\mathbf{C}$ is the lumped damping matrix; $\mathbf{K}$ is the lumped stiffness matrix; ${\mathbf{R}}_{\text{ext}}$ is the applied load vector; $\mathbf{H}$ is the lumped capacitance matrix; and ${\mathbf{I}}_{\text{int}}$ and ${\mathbf{I}}_{\text{ext}}$ are the internal flux vector and applied nodal source vector, respectively. In this work, the “Rayleigh damping” provided in ABAQUS was used ($\mathbf{C}={\beta }_1\mathbf{M}+{\beta }_2\mathbf{K}$). In the following numerical cases, two parameters, ${\beta }_1$ and ${\beta }_2$, are 4000 and 0, respectively.

The equation of motion for the domain (17) is integrated using the explicit central‐difference integration rule[Image Omitted. See PDF]

The equation of fluid flow within fracture (18) is integrated using the explicit forward‐difference time integration rule[Image Omitted. See PDF]

The above model is implemented as VUEL and VEXTERNALDB subroutines in ABAQUS. Since the explicit central‐difference and forward‐difference temporal integration schemes are conditionally stable, the upper limit of the time increment for stability consideration should be provided. The maximum stable time increment $\mathrm{\Delta }t$ of the aforementioned explicit temporal integration algorithm can be estimated[Image Omitted. See PDF]in which $L_c$ is the characteristic length of rock element or cohesive element; $c_d$ is the stress wave velocity in rock; $\overline{\rho }$ and $K_n$ are the linear density and stiffness of cohesive element, respectively; and $w_{\text{avg}}$ is the average aperture of the cohesive element. It should be noted that the available time increment of explicit temporal integration scheme may be very small, which results in expensive computation. However, it has advantage in convergence and stability.

VERIFICATION

The purpose of this section is to verify the capability of the Explicit‐PPCZ model on modeling the growth of hydraulic fracture in naturally fractured reservoir. The simulation results of the KGD problem and interaction behavior of hydraulic fracture crossing natural fracture obtained from the proposed model are compared with the analytical solution and well‐accepted criteria, respectively.

Model verification by the KGD problem

The KGD model is a classical hydraulic fracture propagation model which characterizes the growth of a single, biwing, plane strain fracture driven by injecting fluid into the wellbore at a constant flow rate of $Q_0$, as shown in Figure .

Sketch of KGD model

To compare with the analytical solution, we performed a numerical simulation of half domain due to the symmetry of the problem. The dimension of the domain is 50 × 200 m. Assume that the hydraulic fracture propagates along a predefined path (30 m length) coaxial with x‐axis (see Figure A). The size of the wellbore is relatively small compared to the fracture size. Therefore, the wellbore is modeled as an injection point located at the left side of the predefined path. Figure B shows the generated mesh of the model in which the predefined path and rock matrix are discretized with cohesive elements and triangular elements, respectively. The element sizes increase gradually from 0.5 m in the vicinity of the predefined path to 10 m far away from the fracture. The first two cohesive elements near the injection point are assigned as starter fracture (Figure C) with initial hydraulic opening of 0.02 mm. σ_H and σ_h are imposed, respectively, along the x and y directions on the triangular elements as initial stresses. Symmetric boundary condition is applied to the left edge, and normal displacements at the other boundaries are fixed. Rock properties and other input parameters of the numerical simulation and the analytical solution are listed in Table .

(A) Model geometry; (B) mesh scheme; and (C) zoomed‐in view of the finite element mesh near the injection point

Rock properties and other input parameters

Parameters	Values
Young's modulus (GPa)	30
Poisson's ratio	0.2
Rock density (kg/m3)	2700
In situ stresses (MPa)	1/2/5
Normal stiffness (GPa/m)	200
Tension strength (MPa)	0.1 (viscosity‐dominated)
10 (toughness‐dominated)
Tension fracture energy (J/m2)	20 (viscosity‐dominated)
2000 (toughness‐dominated)
Shear stiffness (GPa/m)	200
Cohesion strength (MPa)	0.1
Shear fracture energy (J/m2)	20
Friction coefficient	0.6
Injection rate (m3/s)	0.001
Viscosity (Pa·s)	0.001
Bulk modulus of fluid (GPa)	2.2
Fluid density (kg/m3)	1000

The simulation results under viscosity‐ and toughness‐dominated conditions are compared with the analytical solution provided by Detournay, respectively. Temporal evolution of the net pressure, fracture aperture at injection point, and fracture length is shown in Figure A–C, respectively. Figure A,B shows the net pressure and fracture aperture along the fracture length at the end of injection. It is obvious that the numerical solutions agree well with the analytical solutions, which indicates that the developed model can be effectively used to simulate the propagation of a single hydraulic fracture.

Time evolutions of (A) net pressure, (B) fracture aperture at injection point, and (C) length of fracture

(A) Net pressure and (B) fracture aperture along the fracture length at the end of injection

Model verification by the interaction between hydraulic fracture and natural fracture

Given that the interaction between hydraulic fracture and natural fracture plays a key role in the growth of hydraulic fracture network in naturally fractured formation, accurately predicting the fracture crossing behavior is essential in the numerical simulation. Therefore, the capability of the proposed Explicit‐PPCZ model to simulate the intersection behavior should be verified before utilizing it to model the evolution of fracture network in naturally fractured formation. Therefore, we performed a series of numerical cases and compared the simulation results with criteria of hydraulic fracture crossing natural fracture proposed by Gu et al and Cheng et al. In this work, a residual hydraulic aperture (ie, 1.0E−5 m) is assigned to the broken PPCZ elements if their mechanical aperture is zero (eg, shear failure).

As illustrated in Figure , a 2D plane strain hydraulic fracture is propagated against the least principle stress (ie, propagating vertically) by injecting a viscous fluid from the injection point, which may intersect the preexisting natural fractures at an approaching angle β. The size of the domain is 10 × 10 m, and the injection point is located at the center of the domain. The coordinates of the expected intersection points between HF and natural fractures are (0, 1.5) and (0, −1.5), respectively. Approaching angles of 90° and 75° are considered in this section, and the corresponding finite element meshes are shown in Figure . The natural fracture surfaces are treated as cohesive elements with relatively small cohesive strength and negligible fracture energy. Input parameters are provided in Table .

Schematic of an HF propagating toward natural fractures at an approaching angle β

Finite element meshes. (A) β = 90° and (B) β = 75°

Rock properties and other input parameters

Parameters	Values
Rock
Young's modulus (GPa)	30
Poisson's ratio	0.2
Density (kg/m3)	2700
Rock fracturing parameters
Normal stiffness (GPa/m)	300
Tensile strength (MPa)	3
Type I fracture energy (J/m2)	10
Shear stiffness (GPa/m)	3000
Cohesion strength (MPa)	30
Type II fracture energy (J/m2)	100
Friction coefficient	0.6
Natural fracture parameters
Normal stiffness (GPa/m)	300
Tensile strength (MPa)	0.3
Type I fracture energy (J/m2)	1
Shear stiffness (GPa/m)	3000
Cohesion strength (MPa)	0.5
Type II fracture energy (J/m2)	10
Fracturing fluid
Injection rate (m3/s)	0.0001
Viscosity (Pa·s)	0.001
Bulk modulus (GPa)	2.2
Density (kg/m3)	1000

Figure shows the fracture profile and displacement field at different injection time with 90° approaching angle, 0.2 friction coefficient, and zero differential stress (σ_H: σ_h = 10:10). The deformation is magnified 300 times for better illustration. It can be observed that the natural fractures are activated when the HF hits them, and then, the HF is arrested by the NF and then propagates along the natural fractures.

Fracture profiles and displacement field at different injection time for the case of zero differential stress. The deformation is magnified 300 times

Figure shows the HF‐NF interaction behavior of a case with a larger differential stress (σ_H:σ _h = 25:10). Three zoomed‐in images at different time are given to illustrate the detail process of the intersection. At t = 26.0 seconds, the HF arrives the intersection points and the HF tips are blunted due to the shear slip of the NFs. The HF crosses the NFs at t = 26.1 seconds as the interfacial forces increase to a critical value which is sufficient to split the rock on the other side of NFs. The NFs remain closed during the injection process, which indicates that HF directly crosses the NFs under this condition.

Fracture profiles and displacement field at different injection time for the case of 15 MPa differential stress. The deformation is magnified 300 times

In addition to the mentioned two cases, more cases with different NF friction coefficient, approaching angle, and differential stress are performed. The fracture crossing results are summarized and compared to the criteria of HF crossing NF proposed by Gu et al and Cheng et al, as shown in Figure . The area in the upper‐right side of the criterion line indicates conditions that a HF will cross a NF when they encounter, while the area in the down‐left side of the criterion line indicates conditions that a HF will turn into a NF. The figure shows great agreement between numerical results and criteria output, which demonstrates the capability of the proposed method to simulate the growth of HF in naturally fractured formations.

Numerical results of interaction between HF and NFs: (A) Approaching angle is 90°; and (B) approaching angle is 75°

HF GROWTH IN NATURALLY FRACTURED FORMATION

Model setup

Naturally fractured formation generally contains several NF sets, each of which comprises fractures that have a similar orientation but different sizes. Within a specific set, fracture intensity shows a size‐dependent power‐law distribution. NF can be usually characterized in terms of its location, orientation, length, mechanical and hydraulic property, etc. The main purpose of this section is to reveal how the existing discontinuities affect the growth of hydraulic fracture. Therefore, the natural fracture networks in this research are generated stochastically, and fractures in the same set are assumed to be the same orientation for simplicity. Fracture number and lengths were derived from a power‐law distribution function [Image Omitted. See PDF]

in which $n\left(l\right)\mathrm{d}l$ is the number of NFs having a length in the range $\left[l,l+\mathrm{d}l\right]$, $\beta$ is a coefficient of proportionality, and $a$ is an exponent varying between 1 and 3.

With the above definition, the number of NFs having a length in the range $\left[l,l+\mathrm{d}l\right]$ in a unit thickness, square domain with side length $L$ may be written as[Image Omitted. See PDF]

Here, $\beta =0.1$, $a=0.2$, and $L=150$ were used to generate the number of fractures with different length, and the fracture length range from 5 to 10 m was employed in the model. In addition, the coordinates of NFs’ center are also needed to draw a NF set. The uniform distribution function was utilized to generate these coordinates in advance. Fractures with an angle of 15° from horizontal (ie, x‐axis) were embedded into a square domain of size 150 × 150 m as joint set 1. To avoid the effect of NFs near the domain's boundaries, a smaller square domain of size 100 × 100 m was cut from the center of the above domain, as shown in Figure A. To investigate the effect of the spatial distribution of NFs on the growth of fracture networks, a domain containing joint set 2 of angle 30° and a domain containing both set 1 and set 2 were generated using the same strategy, as shown in Figure B,C respectively. Some tactics have been introduced to reduce the complexity of mesh generation, the numerical difficulties, and the computational cost. First, all NFs in whole domain are not allowed to intersect at a very low angle (eg, 20°). Second, the distance of nearly parallel nonintersecting NFs must be greater than a user‐specified threshold value (eg, the value of 1 m was employed). Third, the length of the segments of the intersecting NFs was slightly adjusted to make sure that it is either zero or >1 m. The reasons of these tactics are to avoid elements that have high aspect ratio or are very small, which may lead to large computational error or result in extremely small stable time increment. Figure shows the natural fracture networks of the mentioned three cases. It should be pointed out that the finite element meshes of these cases are the same, as shown in Figure . The whole domain is discretized by using 21 384 triangular and 31 876 PPCZ elements having edge length ~1 m. Necessary simulation parameters are listed in Table unless otherwise stated.

Three cases of NF distribution pattern with different orientation and NF set

The mesh of the naturally fractured formation

The injection fracture, located in the center of the domain as black line, is defined having same compressive stiffness as the conventional PPCZ elements, but with initial hydraulic aperture 0.02 mm, null cohesion, and zero tensile strength.

Numerical results

Influence of NFs spatial distribution

To study the influence of spatial distribution of NFs on the evolution of fracture networks, three scenarios shown in Figure are considered. The normal displacements of outer boundaries of the model were fixed. Minimum and maximum horizontal principle stresses, σ_x: σ_y = 10:15 MPa, were imposed on the triangular elements as initial in situ stresses. The fluid injection rate of 2 × 10⁻⁴ m³/s was applied at the center of the domain and a total of 500 seconds of injection was simulated.

It is well known that hydraulic fractures are apt to propagate perpendicular to the direction of minimum horizontal principle stress in homogeneous formation, while in formations with small stress contrast (the difference in the horizontal principle stresses), the geometry of hydraulic fracture may be impacted significantly by the preexisting NFs. Figure shows the simulation results of stimulated fracture network at the end of simulation for the cases with different NFs’ spatial distribution. The fractures (hydraulic fractures and activated natural fractures) that are stimulated during the treatment are shown in blue color, and other inactivated ones are shown in gray. We can observe that the orientation of natural fractures has a significant influence on the overall propagation direction of hydraulic fractures. Furthermore, the copresence of NF sets 1 and 2 increases the tortuosity of hydraulic fracture path, which may significantly limit the proppant placement, consistent with the understanding in hydraulic fracturing community that using fracturing fluid with low sand concentration to prevent screenout of proppants during hydraulic fracturing treatment in naturally fractured reservoirs.

HF growth in naturally fractured formation with different NFs’ spatial distributions

Influence of initial hydraulic aperture of NFs

Not all NFs in formations are nonpermeable, and many of them actually allow fluid to flow through under in situ conditions. That indicates some NFs have certain initial hydraulic aperture, which can significantly affect the stress and pore pressure within NFs during stimulation treatment and consequently influence the shear activation of NFs and thus the evolution of fracture network. To this end, we varied the initial hydraulic aperture of NFs to investigate its effects on hydraulic fracture growth in naturally fractured reservoir. Three cases with different values of that, that is, 0, 0.01, and 0.05 mm, were performed in this study. The mesh, boundary, and initial conditions as well as other input parameters were similar to what the third case has in Section 4.2.1. Figure depicts the fracture path for NFs with different initial hydraulic apertures. The result of the third case in Section 4.2.1 with zero initial hydraulic aperture is also plotted in Figure A as reference. Different colors represent the Explicit‐PPCZ elements failed in different modes: red for tensile (mode I) and blue for shear (mode II), respectively. The simulation results show that the NFs with a very small initial hydraulic aperture have high potential of being activated during hydraulic fracturing treatment, which in turn indicates that it is more difficult to form complex fracture networks in reservoir with nonpermeable NFs under the same conditions. Although the mechanical aperture of these shear‐activated NFs in Figure A,C are too small to accommodate any proppant, their hydraulic aperture may increase due to the shear dilation of NFs’ rough walls. If hydraulic conductivity increase of these shear‐activated NFs remains for a relatively long period of time, the effective producing fracture area will increase significantly, which is critical for economical unconventional resources development.

HF growth in naturally fractured formation with different NFs’ initial hydraulic apertures

Influence of stress contrast

In situ stress contrast, (σ_H – σ_h), has a crucial influence on the activation of natural fractures and consequently the growth of complex fracture network in naturally fractured reservoirs. In this study, four cases with different stress contrasts were performed to investigate its impact on the evolution of fracture network. The initial hydraulic aperture of NFs was assumed as 0.1 mm, and the simulation results are provided in Figure in a similar manner as the previous section. Figure A,B shows the stimulated fracture network for the cases of low stress contrast (ie, 0 and 2 MPa, respectively), in which the fracture network has preferred growth direction in the upper wing along joint set 2. This is probably because that the low stress contrast state has little impact, and consequently, the overall fracture network propagation direction is dominated by the preexisting natural fractures. For the cases with high stress contrast, the fracture network propagation direction is dominated by the in situ stress and is along the direction perpendicular to the direction of the minimum principle stress, as shown in Figure C,D. In addition, it is also worth mentioning that the activation mechanisms of natural fractures under high and low stress contrasts are different. As shown in Figure , it is clear that the activation mechanism of fracture network is dominated by tensile failure at low stress contrast, while it transitions to the shear failure dominated at high stress contrast. There is a widely accepted consensus in hydraulic fracturing community that low stress contrast facilitates the formation of complex fracture network. Our simulation results indicate that high stress contrast is also favorable for the formation of complex fracture network.

HF growth in naturally fractured formation with different stress contrasts

CONCLUSIONS

In this paper, a 2D, fully coupled hydromechanical model has been developed for a better understanding of the hydraulic fracture propagation in naturally fractured formation. In the proposed method, the fracture propagation and fluid flow within hydraulic fracture are governed by cohesive zone model equipped with Coulomb's friction law and cubic law, respectively. The fracturing fluid leak‐off into the rock matrix is ignored. The explicit temporal integration scheme for discretizing and solving the governing equations is adopted to improve the convergence during the simulation. The capability of the model was then verified with the analytical solution of KGD problem and well‐accepted criteria of hydraulic fracture crossing natural fracture. Finally, the model was used to investigate three important factors that affect the evolution of fracture network in naturally fractured formation. The following conclusions can be obtained from the simulation results.

1. The orientation of natural fractures has a significant influence on the overall propagation direction of hydraulic fractures. Furthermore, the copresence of NF sets increases the tortuosity of hydraulic fracture path, which may significantly limit the proppant placement. The numerical results support the understanding in hydraulic fracturing community that using fracturing fluid with low sand concentration to prevent screenout during hydraulic fracturing treatment in naturally fractured reservoirs.
2. Fracturing fluid leak‐off and associated pressure increasing in natural fractures play a key role in the activation of natural fractures. It is more difficult to obtain complex fracture networks if the natural fractures in reservoir are all nonpermeable.
3. Both low and high in situ stress contrasts facilitate the formation of complex fracture network with the presence of natural fractures. The difference between these two stress conditions is that the activation of fracture network is dominated by tensile failure at low stress contrast, while it transitions to the shear failure dominated at high stress contrast.

ACKNOWLEDGMENTS

This work was financially supported by National Natural Science Foundation of China (Nos. 51804268, 51874253, and 51604232) and Sichuan Science and Technology Program (2019YJ0521 and 2018FZ0069).