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

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.

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 10 and Lecampion et al. 11 The earliest analytical 2D and 3D models are KGD, 12 PKN, 13,14 and penny-shaped, 15,16 which have now been largely replaced by pseudo-3D models. 10 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, [17][18][19][20] cohesive element, [21][22][23][24][25] and the partition of unity methods (eg, XFEM and GFEM) [26][27][28][29][30][31][32][33][34] 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, 35 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 36 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 37 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 38 combined PFC and FLAC to study the behavior of a hydraulic fracture crossing natural fractures. Fatahi et al 39 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 40 and Barenblatt,41 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. 42 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, 21,[43][44][45][46] 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 1A, 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 1B. 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.

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

| Deformation of rock
The rock is assumed as an isotropic, homogenous, and elastic medium, and thus, the momentum balance can be expressed as 47 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 48 where D ijst is the elasticity material property matrix of the triangular element and ε ij is the infinitesimal strain tensor defined as ij = 1 2 u i,j + u j,i .

| Fluid flow within fracture
Fluid flow within a fracture can be described by the fluid mass conservation equation 49 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 50 where p 0 is the datum pressure, which is assumed to be zero here; ρ 0 is the fluid density at datum pressure; and K f is the bulk modulus of the fluid.
Fluid flow within smooth parallel fracture yields 51 where μ is the viscosity of fluid.

| Cohesive zone model
Cohesive zone model (CZM) is a phenomenological model introduced by Barenblatt 41 and Dugdale 40 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 tractionseparation law. Further details on cohesive zone models can be found in Refs. 42,52 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 3, we utilize the following CZM that combined the bilinear traction-separation law 42 and Coulomb friction law 53   n and 1 t are the normal and the shear relative displacements at the complete failure. f is the interfacial friction. In Figure 3A, the area below the traction-separation curve is denoted by G IC , which represents the fracture energy per unit area of the fracture surface.
A quadratic criterion is implemented to determine the mixed-mode failure, 54 as follows where G IC and G IIC are tensile and shear energy release rates, respectively; G I and G 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 55 in which K IC and K IIC are mode I and mode II fracture toughness, respectively; E is Young's modulus; and v is the Poisson ratio of rock.

| 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  ) where ε, u, and 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 p from left and then integrating the product over the domain Γ HF Using integration by parts and divergence theorem, the right-hand side of Equation (14) can be recast as Substituting Equation (15) into Equation (14), the weak form of the equation governing the fluid flow within fracture can be rewritten as where Γ 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 Ω and Γ HF , respectively, the discrete form of coupled governing Equations (13) and (16) can be given by where M is the lumped mass matrix; C is the lumped damping matrix; K is the lumped stiffness matrix; R ext is the applied load vector; H is the lumped capacitance matrix; and I int and I ext are the internal flux vector and applied nodal source vector, respectively. In this work, the "Rayleigh damping" provided in ABAQUS was used (C = 1 M + 2 K). In the following numerical cases, two parameters, 1 and 2 , are 4000 and 0, respectively.
The equation of motion for the domain (17) is integrated using the explicit central-difference integration rule The equation of fluid flow within fracture (18) is integrated using the explicit forward-difference time integration rule 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 Δt of the aforementioned explicit temporal integration algorithm can be estimated 56 in which L c is the characteristic length of rock element or cohesive element; c d is the stress wave velocity in rock; ̄ and K n are the linear density and stiffness of cohesive element, respectively; and w 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.

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 4.
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  5A). 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 5B 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 5C) 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  Table 1.
The simulation results under viscosity-and toughnessdominated conditions are compared with the analytical solution provided by Detournay, 57 respectively. Temporal evolution of the net pressure, fracture aperture at injection point, and fracture length is shown in Figure 6A-C, respectively. Figure 7A,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.

| 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 58 and Cheng et al. 59 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 8, 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 9. The natural fracture surfaces are treated as cohesive elements with relatively small cohesive strength and negligible fracture energy. Input parameters are provided in Table 2. Figure 10 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. Figure 11 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.
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 58 and Cheng et al,59 as shown in Figure 12. 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.

| 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. 60 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 [60][61][62]  in which n (l) dl is the number of NFs having a length in the range [l,l + dl], 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 [l,l + dl] in a unit thickness, square domain with side length L may be written as Here, = 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 13A. 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 13B,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 13 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 14. 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 2 unless otherwise stated.
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.

| 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 13 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 −4 m 3 /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 15 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. 63

| 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

F I G U R E 1 4
The mesh of the naturally fractured formation 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 16 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 16A 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  16A,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 shearactivated 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.

| 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 17 in a similar manner as the previous section. Figure 17A,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 17C,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 17, 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. 9 Our simulation results indicate that high stress contrast is also favorable for the formation of complex fracture network.

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