Numerical Investigation of Fracture Network Formation under Multiple Wells

A two-step fracturing method is proposed to investigate the hydraulic fracture evolution behavior and the process of complex fracture network formation undermultiple wells. Simulations are conducted with Rock Failure Process Analysis code. Heterogeneity and permeability of the rocks are considered in this study. In Step 1, the influence of an asymmetric pressure gradient on the fracture evolution is simulated, and an artificial structural plane is formed. *e simulation results reflect the macroscopic fracture evolution induced by mesoscopic failure; these results agree well with the characteristics of the experiments. Step 2, which is based on the first step, investigates the influence of preexisting fractures (i.e., artificial structural planes) on the subsequent fracturing behavior. *e simulation results are supported by mechanics analysis. Results indicated that the fracture evolution is influenced by pressure magnitude on a local scale around the fracture tip and by the orientation and distribution of pore pressure on a global scale. *e constant pressure in wellboreH2 can affect fracture propagation by changing the water flow direction, and the hydraulic fractures will propagate to the direction of higher local pore pressure. Furthermore, the artificial structural planes influence the stress distribution surrounding the wellbores and the hydraulic fracture evolution by altering the induced stresses around the preexisting fractures. Finally, fracture network is formed among the artificial structural planes and hydraulic fractures when multiple wells are fractured successively. *is study provides valuable guidance to unconventional reservoir reconstruction designs.


Introduction
In 1947, the first experimental hydraulic fracturing treatment in the United States occurred in the Hugoton gas field in Grant County, Kansas [1]. Since then, this technology has been used worldwide, and it has become an important technical means for oil and gas exploitation from fractures. However, in recent years, conventional oil and gas have gradually been exhausted due to exploitation [2]. Moreover, unconventional oil and gas resources, which are abundant in widely distributed reserves with low permeability, have become a global energy focus [3]. Given that unconventional reservoirs (e.g., shale reservoirs) are characterized by wide distributions, large specific surfaces, low porosity, low permeability, and poor connectivity, such reservoirs cannot be efficiently and economically developed through a single fracture. Studies have shown that fracture development has a substantial influence on the permeability of reservoirs. As fractures develop in a rock mass, even lowporosity areas can have high permeability; by contrast, in highporosity zones, permeability can be low due to poor connectivity and lack of fractures [4]. erefore, methods of developing hydraulic fractures and forming complex fracture networks to enhance stimulated reservoir volume (SRV) are among the first considerations in fracture design [5][6][7]. Reconstructing unconventional reservoirs is necessary to improve their permeability, weaken the adsorption effect of organic carbon on oil [8], increase the drainage area, and improve oil and gas production. e exploitation of unconventional oil and gas resources is a current topic of interest. In recent years, many scholars have studied the hydraulic fracture evolution and the complex fracture network formation and obtained certain results that are positive to unconventional reservoir reconstruction. Yan et al. [9] studied the effects of crustal stress field, confining pressure, and natural fractures on the fracture initiation and propagation by laboratory tests. e results demonstrated that stress concentration around the hole would significantly increase the fracture pressure of the rock, and natural fractures in the borehole wall would eliminate stress concentration. He et al. [10] investigated the different hydraulic fracture extension patterns of shale through hydraulic fracturing experiments, and they believed that the typical bedding plane well developed in the shale formation plays an important role in the propagation of hydraulic fractures. Hou et al. [11] proposed stimulated rock area (SRA) as an evaluation index for hydraulic fracturing results, and they found that lower in situ stress difference and shorter distance between hydraulic fracture and bedding plane in brittle shale formation lead to a larger SRA and more complex fracture morphology.
A well is a prerequisite for oil and natural gas exploitation. ousands or even tens of thousands of wells are found in one oil field. Simultaneous fracturing in two or more offset wells can effectively improve production [12]. Multiwell fracturing is important for generating complex fracture networks in unconventional reservoirs. A pore pressure gradient is formed in a certain range by controlling the water pressure that is pumped into different wells to control the initiation, coalescence, and propagation of hydraulic fractures. Fractures are influenced by the pore pressure value on a local scale around the crack tip and by the orientation and distribution of pore pressure gradients on a global scale. e fractures will propagate toward the regions of high local pore pressure [13].
is pore pressure gradient also aids in the formation of a fracture network in a low permeability reservoir. Yao et al. [14] investigated the propagation regularity of hydraulic fractures in the mode of multiwell pads. ey found that when multiple wells are simultaneously fractured, adjacent fractures will propagate toward one another, and multiwell fracturing can induce larger areas of stress reversal compared to single-well fracturing. Sesetty and Ghassemi [15] developed a 2D coupled displacement discontinuity numerical model to simulate fracture propagation in simultaneous and sequential hydraulic fracture operations for single and multiple parallel wells. ey found that fracture path is affected not only by fracture spacing but also by the boundary conditions on the previously created fractures. e present study aims to numerically investigate how the pore pressure field can affect fracture propagation and determine how the fracture network forms when multiple wells are successively fractured in rock materials using the Rock Failure Process Analysis (RFPA) 2D 2.0-Flow code. Two steps in hydraulic fracturing are modeled as examples to illustrate the pore pressure distribution that affects fracture initiation and the propagation and preexisting fractures that influence the fracture behavior and stress distribution by altering the induced stresses.

RFPA-Flow Code
Different from metal and glass materials, rock is a heterogeneous natural geological material [16]. erefore, the traditional calculation method based on material homogeneity and isotropy is not suitable for rocks. Rock failure is an accumulation process from mesoscale weakening to macroscale damage [17][18][19], which can be described as nonlinear but not elastic-plastic. RFPA-Flow is a code developed for heterogeneous and permeable geomaterials (e.g., rocks). It can simulate the fracture propagation and failure accumulation process of quasi-brittle heterogeneous materials [20][21][22][23][24][25][26][27][28][29]. Finite element method is used in the RFPA-Flow code. In addition, the plane four-node isoparametric element is the basic element.
e geomaterials (rock) modeled in RFPA-Flow code are assumed to be composed of many mesoscopic elements, which are used to represent the heterogeneity of materials and the random distribution of defects in the geomaterials. e mesoscopic elements in the numerical model are assumed to be isotropic and homogeneous, and their mechanical properties (e.g., constitutive relation, AE rate, and loading rate relation) are assumed to linearly vary. ese mesoscopic elements are statistically distributed (e.g., normal, Poisson, and Weibull distributions) to describe the mechanical properties of the nonlinear macroscopic model (Figure 1). e failure strength in rock can vary remarkably because of grain-scale heterogeneity. e material properties of different mesoscopic elements are randomly distributed throughout the domain of analysis based on a Weibull distribution to include the statistical variability of the bulk failure strength in the RFPA model [30]. (1) Equation (1) represents the statistical distribution f(μ) of the mechanical parameters of the mesoscopic elements in the RFPA-Flow code. e heterogeneity index m [31,32] is a parameter that defines the shape of the distribution function and the degree of material heterogeneity. A large m value indicates the presence of highly homogeneous materials, whereas a small m value denotes the existence of highly inhomogeneous materials. In addition, μ represents the mechanical parameters (i.e., Young's modulus, Poisson's ratio, tensional strength, and compressive strength) of the elements, and μ 0 denotes the scale parameters related to the average values of the mechanical parameters. A high homogeneity index value indicates that the values of most elements are concentrated closer to μ 0 . In the simulation, the mesoscopic elements are defined as damaged when the strength criterion is met. Numerous elements fail as the stress increases; and as these elements become connected, the heterogeneous materials can fail.
In the RFPA model, matrix, air, and contact primitives are used to describe the deformed accumulation and failure processes of heterogeneous materials ( Figure 2). Primitive is another expression of the RFPA mesoscopic element. Primitive phase transition occurs when the mesoelement is in tension or compression. If the element deformation is in elastic and residual deformation (I and II in Figure 2), then the mesoelement is the matrix primitive. If the element deformation exceeds the residual deformation (III in Figure 2), then the matrix primitive changes to air primitive (apart phase) or contact primitive (contact phase). e element mechanical properties change with the occurrence of the primitive phase transition.
In RFPA-Flow code, the element fails when the matrix primitive transfers to air or contact primitive (tension or compression failure). In the simulation, the parameters of the element are multiplied by an extremely small coefficient to achieve element failure, and the failure element is presented as black in the model. e material properties of different mesoscopic elements in the RFPA model are randomly distributed; thus, the elements fail successively and connect gradually to form fractures as the load increases. en, the fracture propagation and failure accumulation of quasi-brittle heterogeneous materials are obtained.
In the RFPA-Flow code, the geological medium (rock) is assumed to be fully saturated with fluid flow governed by Darcy's law. In addition, the coupled process of stress and seepage in the deforming rock mass is governed by Biot's theory of consolidation. In view of the effects of stress and damage on permeability, the basic formulations used in the analysis are as follows: (2) In the above equations, ρ denotes density; σ ij ′ is the effective stress; σ ij is the total stress; ε ij is the total strain; ε v is the volume strain; α is the coefficient of the pore pressure; U i is the displacement; p is the pore pressure; λ is the Lame coefficient; δ ij is the Kronecker delta; G is the modulus of shear deformation; Q is Biot's constant; K is the permeability coefficient; K 0 is the initial permeability coefficient; β is the coupling parameter, which reflects the influence of stress on the coefficient of permeability; and ξ(ξ ≥ 1) is the mutation coefficient of permeability, which accounts for the increase in permeability when the element reaches the damage state. e values of coefficients ξ, α, and β are determined experimentally, and they vary with the stress state and damage evolution of the rock. In the RFPA-Flow code, 0 ≤ α < 1 and ξ � 1 are assumed for the mesoscopic element in the elastic stage (as indicated by D � 0 in equation (7)). Once damage occurs, the permeability of the mesoscopic element increases remarkably (equations (9) and (10)). Furthermore, ξ � 5 is assumed for the damaged element (denoted by 0 < D < 1) and α � 1; ξ � 100 is assumed for the fully damaged element (denoted by D � 1).
In the elastic state, the relationship of stress and permeability coefficient is described by where σ ′ is the effective stress and b is the coupling parameter.
Continuum damage mechanics are used to describe the constitutive laws of the mesoscopic elements in RFPA. Initially, an element is considered as elastic, and the associated elastic properties can be defined by Young's modulus and Poisson's ratio. e stress-strain relationship of each element is considered as linearly elastic until the given damage threshold is reached, and tensile and shear failures are considered in the analysis (Figure 3). An element is considered to have failed in the tensile mode when the maximum tensile stress criterion is satisfied (equation (4)), and failure in shear mode occurs when the shear stress satisfies the Mohr-Coulomb failure criterion (equation (5)).
In equations (4) and (5), f t is the tensile failure strength, ϕ is the friction angle, and f c is the compression failure strength.
For an individual element, when the stress of the element satisfies the selected strength criterion, the element begins to undergo damage. On the basis of isotropic elastic damage theory, the elastic modulus of an element may gradually degrade as damage progresses, and the elastic modulus of the damaged material can be defined as follows: where D represents the damage variable; and E and E 0 are the elastic moduli of the damaged and undamaged elements, respectively. In the simulation, when the tensileshear stress in an element reaches the failure threshold (equations (4) and (5)), the damage variable can be expressed as follows: where f tr is the residual tensile strength; ε t0 is the maximum tensile strain at the elastic limit; ε tu is the ultimate tensile strain of the element, at which the element is completely damaged; f cr is the residual compressive strength; and ε c0 is the maximum compressive strain.
As the damage variable D of the element changes in the damage process, the permeability coefficient of the mesoscopic element is expressed as

Induced Stresses around a Pre-Existing Fracture
Many natural and manmade fractures are distributed in rock masses, and these fractures can influence the initiation and propagation of new fractures [33]. e interactions among fractures can change the local in situ stress and further affect the propagation direction of the fractures in the bedrock [34]. On the basis of the semi-infinite fracture model proposed by Green and Sneddon [35], the induced stresses around a preexisting fracture were analyzed. As shown in Figure 4, given the existence of a fracture, the induced stresses at any point around the fracture are as follows: where p net is the net pressure and η is a coefficient that describes the relationship among σ x , σ y , and σ z . e induced stresses substantially change with the distance to the fracture (equations (11)- (13)). e induced stresses in the direction perpendicular to the fracture are larger than those in the direction of the fracture. In the same direction, the smaller the distance is, the greater the induced stresses will be. Accordingly, as the propagating fracture approaches the preexisting fracture, such fracture propagates more easily.

Numerical Model Establishment
In this study, the numerical simulation of fracture network formation can be divided into two steps.
Step 1. It strictly follows the experiment by Bruno and Nakagawa [13], which can verify the numerical results and obtain the corresponding stress and pore pressure distribution that the test cannot achieve. As shown in Figure 5(a), the dimension of the model is 152 × 152 mm. e model is divided to form a 400 × 400 mesh. e radius of wellbores H 1 , H 2 , and H 3 is 4 mm. e distance between H 1 and H 2 is 50 mm, and that between H 2 and H 3 is 40 mm. A vertically induced crack (0.6 mm wide and 6 mm long) is found at the bottom of H 1 to provide a starting direction for the hydraulic fracture. In the experiment, the injection rate of wellbore H 1 varies from 0.5 cm 3 /min to 3.0 cm 3 /min, and the influence is insignificant. erefore, in the simulation, the initial water pressure is 1 MPa, and this pressure is increased at an increment of 0.05 MPa per step in wellbore H 1 . To compare the fracture mode with different pore pressure gradients, constant pressures of 0.6, 1.0, 1.4, and 1.7 MPa are sequentially applied in wellbore H 2 in the four models. No pressure is applied in H 3 . In the experiment, the thickness of the slab is 25 mm. us, the plane-strain calculation is used in the simulation. In addition, the boundary displacement is limited. us, a displacement load with 0 is applied to the model boundary.
Step 2. It is based on Step 1. As shown in Figure 5 ese planes are based on the results of Step 1. e two induced cracks on both sides of H 3 provide a starting direction for the hydraulic fractures, and the angles between the cracks and horizontal direction are 45°and 135°, respectively. An initial water pressure of 1 MPa and an increasing water pressure at an increment of 0.05 MPa per step are applied via wellbore H 3 . Table 1 shows the rock parameters used in the models [13,36], and Table 2 shows the model settings of the two steps.

Fracture Propagation Trends under an Asymmetric Pressure Gradient
As shown in Figure 6(a), the asymmetric pressure gradient is established after the pressure is applied via H 1 and H 2 . As the pressure at H 2 increases from 0.6 MPa to 1.7 MPa, the -f tr asymmetry of the pore pressure gradient is enhanced (Step 1, Figure 6(a)). As the pressure at H 1 increases, the fracture initiates at the induced crack tip, and crack propagation is influenced by the pore pressure. Conversely, the fracture propagation affects the pore pressure evolution. As the fracture propagates, the fluid flow in the fracture causes the asymmetry of the pore pressure gradient to be further enhanced. When the model loses stability, the calculation terminates. e fracture is generally caused by tensile failure, and the fracture propagation direction is perpendicular to the tensile stress direction. Figure 6(b) shows that the stress distribution of the specimen is controlled by the pore pressure. As the pressure in wellbore H 2 increases, the tensile stress region around H 1 and H 2 (the light area around the well) expands (Step 1, Figure 6(b)). Furthermore, the pressure in well H 2 is higher, and the fracture is closer to wellbore H 2 (final step in Figure 6(b)). All the fractures in the four models are jagged fractures, and such fractures exhibit typical characteristics of tensional failure and rock heterogeneity because a fracture will deterministically select the path of least energy resistance through a rock. When the pressure in well H 2 is 1.7 MPa, the maximum tensile stress is perpendicular to the line between H 1 and H 2 , and a large deflection can be observed in the fracture propagation direction. erefore, the pore pressure value and the asymmetry of the pore pressure gradient directly influence the tensile fracture initiation and propagation in the rock. Figure 7 shows that the numerical results are consistent with the experimental results [13]. e numerical results of fracture evolution are generally consistent with the experimental results at the macroscale for different asymmetric pore pressure gradients. In addition, for different pore pressure fields, the breakdown pressure increases linearly with increasing pressure at H 2 (Figure 8). e variability in the breakdown pressure indicates that the pore pressure gradient is the main factor that controls model stability.

Fracture Network Formation Based on Multiple Wells
On the basis of the numerical simulation results of Step 1, a constant pressure of 1.7 MPa is applied at H 2 (H 4 ) as well as increasing water pressures at H 1 (H 5 ). As the water pressure increases, the fracture initiates and propagates at the induced crack tip of H 1 (H 5 ) and finally connects to H 2 (H 4 ) ( Figure 5(b)). As the artificial structural planes of Fracture 1 (Fracture 2) are fractured, all injection stops. Subsequently, Step 2 of the simulation starts. As an increasing pressure is applied via well H 3 , three high-stress regions appear around well H 3 (Figure 9(a)). Two regions are located at the induced crack tips, and the other region is at the bottom of H 3 . As the pressure in H 3 increases, a fracture initiates at the induced crack tip (under the induced stresses) and deviates from the original direction, and a nonuniform pore pressure field is established (Figure 9(f )). Moreover, the two high-stress regions move with the crack tips. en, three high-stress regions form an equilateral triangle (Figure 9(b)), and the angle between each high-stress point and the well center is approximately 120°. ese regions form because of the hydrostatic pressure on the round wellbore and the heterogeneity of the rock materials.
Given the existence of Fractures 1 and 2, the induced stresses around the two preexisting artificial structural planes alter the stress distribution surrounding well H 3 and further influence fracture propagation. As the hydraulic pressure increases, the hydraulic fractures initiate and propagate at the tips of the two induced cracks (high-stress Regions 1′ and 2′) that are closer to preexisting Fractures 1 and 2 than the bottom of well H 3 (high-stress Region 3′; Figures 9(b) and 9(c)). In addition, during the propagation process of the two hydraulic fractures, the high stress in Region 3′ is gradually transferred to Regions 1′ and 2′ (Figure 9(c)). us, fracture propagation occurs appropriately when new fractures are close to preexisting fractures, because the induced stresses are large when the distance to the nearest preexisting fracture is short (Figures 9(f ) and 9(g)). en, two hydraulic fractures propagate at the crack tips as the pore pressure increases and perpendicularly contacts Fractures 1 and 2 because of the effects of the largely induced stresses in the direction perpendicular to Fractures 1 and 2 (Figures 9(d) and 9(h)). A fracture network can form if multiple wells are sequentially fractured to create effective channels for oil and gas flow, enhance the permeability of the rock mass, and improve reservoir production.

Discussion
In this study, a two-step fracturing method is proposed to investigate the hydraulic fracture initiation and propagation behavior using the RFPA-Flow code, and a fracture network is obtained. In Step 1, the simulation results show that the asymmetric pressure gradient controls the fracture propagation in gradient scale, and the fractures will propagate toward the regions of higher local pore pressure; these results are consistent with those of Bruno and Nakagawa [13]. e notion of artificial structural plane is proposed in the  Figure 5: Schematic of the model steps: Step 1 (a) and Step 2 (b). Step 1 Final step  is study is based on the existing wells that do not require drilling. In comparison with other multiwell fracturing methods, the proposed approach can save fracturing cost and utilize existing wells.
is study can provide   valuable guidance to the unconventional reservoir reconstruction designs. To obtain a more complex fracture network, third and fourth steps can be added in future research.

Conclusions
(i) e pore pressure value and the nonuniform pore pressure gradient directly influence the tensile Mathematical Problems in Engineering 9 fracture initiation and propagation in rock. e fracture will propagate in the direction of higher local pore pressure. (ii) Asymmetric fracture propagation can increase the asymmetry of the pore pressure gradient because of the fluid flow in the propagating fracture. e asymmetrical pressure gradient can also influence the fracture propagation direction. e pore pressure gradient is the main controlling factor of model stability. e breakdown pressure of the specimens increases as the asymmetry of the pore pressure gradient is enhanced. (iii) By altering the locally induced stress around the preexisting fractures, the preexisting fractures influence the initiation position and propagation direction of the hydraulic fractures under hydrostatic pressure. (iv) A fracture network can be formed when multiple wells are successively fractured. e influence of the pore pressure on tensile fracture propagation must be understood to design efficient fracturing scenarios and production operations in an oil field.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.