Mathematical Problems in Engineering
Volume 2020 (2020), Article ID 1763713, 11 pages
https://doi.org/10.1155/2020/1763713

Research Article

Numerical Investigation of Fracture Network Formation under Multiple Wells

Xiaoxi Men¹ and Jiren Li²

¹College of Civil Engineering, Dalian Jiaotong University, Dalian 116028, China
²College of Civil Engineering, University of Science and Technology Liaoning, Anshan 114051, China

Correspondence should be addressed to Xiaoxi Men; menxiaoxi@126.com and Jiren Li; lijiren@ustl.edu.cn

Received 13 March 2020; Revised 18 June 2020; Accepted 15 July 2020; Published 20 August 2020

Academic Editor: Luca Chiapponi

Copyright © 2020 Xiaoxi Men and Jiren Li. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

A two-step fracturing method is proposed to investigate the hydraulic fracture evolution behavior and the process of complex fracture network formation under multiple 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. The 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. The 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. The constant pressure in wellbore H₂ 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. This study provides valuable guidance to unconventional reservoir reconstruction designs.

1. 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 low-porosity areas can have high permeability; by contrast, in high-porosity zones, permeability can be low due to poor connectivity and lack of fractures [4]. Therefore, 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–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.

The 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. The 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. Thousands 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. The fractures will propagate toward the regions of high local pore pressure [13]. This 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. They 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. They found that fracture path is affected not only by fracture spacing but also by the boundary conditions on the previously created fractures.

The 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)2D2.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.

2. Brief Introduction to the NumericalRFPA-Flow Code

Different from metal and glass materials, rock is a heterogeneous natural geological material [16]. Therefore, 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–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–29]. Finite element method is used in the RFPA-Flow code. In addition, the plane four-node isoparametric element is the basic element.

The 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. The 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. These mesoscopic elements are statistically distributed (e.g., normal, Poisson, and Weibull distributions) to describe the mechanical properties of the nonlinear macroscopic model (Figure 1).

Figure 1: Mechanical properties of RFPA model using a mesoscale linear heterogeneous distribution to describe nonlinearity at the macroscale.

The failure strength in rock can vary remarkably because of grain-scale heterogeneity. The 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].

Equation (1) represents the statistical distribution of the mechanical parameters of the mesoscopic elements in the RFPA-Flow code. The heterogeneity index [31, 32] is a parameter that defines the shape of the distribution function and the degree of material heterogeneity. A large value indicates the presence of highly homogeneous materials, whereas a small 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 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 . 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). The element mechanical properties change with the occurrence of the primitive phase transition.

Figure 2: 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. The 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. Then, 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:

In the above equations, denotes density; is the effective stress; is the total stress; is the total strain; is the volume strain; is the coefficient of the pore pressure; is the displacement; is the pore pressure; is the Lame coefficient; is the Kronecker delta; is the modulus of shear deformation; is Biot’s constant; is the permeability coefficient; is the initial permeability coefficient; is the coupling parameter, which reflects the influence of stress on the coefficient of permeability; and is the mutation coefficient of permeability, which accounts for the increase in permeability when the element reaches the damage state. The 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, and are assumed for the mesoscopic element in the elastic stage (as indicated by in equation (7)). Once damage occurs, the permeability of the mesoscopic element increases remarkably (equations (9) and (10)). Furthermore, is assumed for the damaged element (denoted by ) and ; is assumed for the fully damaged element (denoted by ).

In the elastic state, the relationship of stress and permeability coefficient is described by

where is the effective stress and 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. The 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)).

(a)

(b)

Figure 3: Elastic-brittle damage constitutive laws of mesoscopic elements subjected to uniaxial tensile stress (a) and uniaxial compression (b).

In equations (4) and (5), is the tensile failure strength, is the friction angle, and 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 represents the damage variable; and and are the elastic moduli of the damaged and undamaged elements, respectively. In the simulation, when the tensile-shear stress in an element reaches the failure threshold (equations (4) and (5)), the damage variable can be expressed as follows:

where is the residual tensile strength; is the maximum tensile strain at the elastic limit; is the ultimate tensile strain of the element, at which the element is completely damaged; is the residual compressive strength; and is the maximum compressive strain.

As the damage variable of the element changes in the damage process, the permeability coefficient of the mesoscopic element is expressed as

3. 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]. The 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 is the net pressure and is a coefficient that describes the relationship among , , and . The induced stresses substantially change with the distance to the fracture (equations (11)–(13)). The 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.

Figure 4: Schematic of the induced stresses around a preexisting fracture.

Table 1: Model parameters.

Parameter
Unit
Value

Heterogeneity coefficient (m)

2

Young’s modulus (E)
GPa
20

Poisson’s ratio (ν)

0.25

Friction angle (φ)
°
30

Compressive strength (fc)
MPa
80

Tensile strength (ft)
MPa
8

Permeability coefficient (K)
md⁻¹
0.000864

Table 2: Model setting parameters.

Model
Size (mm)
Mesh
Wellbores
Radius of wellbores (mm)
Induced crack
Water pressure (MPa)

Step 1 model
152 × 152
400 × 400
H₁, H₂, H₃
4
0.6 mm wide and 6 mm long
H₁: 1 initial pressure + 0.05 increment
H₂: 0.6, 1.0, 1.4, or 1.7

Step 2 model
152 × 192
400 × 600
H₁, H₂, H₃, H₄, H₅
4
0.6 mm wide and 6 mm long
H₃: 1 initial pressure + 0.05 increment

4. 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. The model is divided to form a 400 × 400 mesh. The radius of wellbores H₁, H₂, and H₃ is 4 mm. The distance between H₁ and H₂ is 50 mm, and that between H₂ and H₃ is 40 mm. A vertically induced crack (0.6 mm wide and 6 mm long) is found at the bottom of H₁ to provide a starting direction for the hydraulic fracture. In the experiment, the injection rate of wellbore H₁ varies from 0.5 cm³/min to 3.0 cm³/min, and the influence is insignificant. Therefore, 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₁. 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₂ in the four models. No pressure is applied in H₃. In the experiment, the thickness of the slab is 25 mm. Thus, the plane-strain calculation is used in the simulation. In addition, the boundary displacement is limited. Thus, a displacement load with 0 is applied to the model boundary. Step 2. It is based on Step 1. As shown in Figure 5(b), the model dimensions are 152 × 192 mm. The model is divided into a 400 × 600 mesh. The geometrical configuration of wells H₁, H₂, H₃, H₄, and H₅ and the physical parameters are the same as those in Step 1. Two artificial structural planes (Fractures 1 and 2) are created between H₁ and H₂ and between H₄ and H₅. These planes are based on the results of Step 1. The two induced cracks on both sides of H₃ 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₃. Table 1 shows the rock parameters used in the models [13, 36], and Table 2 shows the model settings of the two steps.

(a)

(b)

Figure 5: Schematic of the model steps: Step 1 (a) and Step 2 (b).

5. 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₁ and H₂. As the pressure at H₂ increases from 0.6 MPa to 1.7 MPa, the asymmetry of the pore pressure gradient is enhanced (Step 1, Figure 6(a)). As the pressure at H₁ 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.

(a)

(b)

Figure 6: Numerical results of pore pressure field and minimum principal stress distribution. (a) Pore pressure field. (b) Minimum principal stress distribution.

The 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₂ increases, the tensile stress region around H₁ and H₂ (the light area around the well) expands (Step 1, Figure 6(b)). Furthermore, the pressure in well H₂ is higher, and the fracture is closer to wellbore H₂ (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₂ is 1.7 MPa, the maximum tensile stress is perpendicular to the line between H₁ and H₂, and a large deflection can be observed in the fracture propagation direction. Therefore, 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]. The 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₂ (Figure 8). The variability in the breakdown pressure indicates that the pore pressure gradient is the main factor that controls model stability.

(a)

(b)

Figure 7: Comparison of numerical and experimental results. (a) H₂ = 1.4 MPa. (b) H₂ = 1.7 MPa.

Figure 8: Breakdown pressure of specimens at different water pressures in H₂.

6. 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₂ (H₄) as well as increasing water pressures at H₁ (H₅). As the water pressure increases, the fracture initiates and propagates at the induced crack tip of H₁ (H₅) and finally connects to H₂ (H₄) (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₃, three high-stress regions appear around well H₃ (Figure 9(a)). Two regions are located at the induced crack tips, and the other region is at the bottom of H₃. As the pressure in H₃ 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. Then, 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°. These regions form because of the hydrostatic pressure on the round wellbore and the heterogeneity of the rock materials.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 9: Numerically obtained evolution of minimum principal stress ((a)–(d)) and pore pressure field ((e)–(h)).

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₃ 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₃ (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)). Thus, 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)). Then, 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.

7. 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]. The notion of artificial structural plane is proposed in the present work. Different from a natural structural plane, an artificial structural plane can be set in any favorable position artificially. In Step 2, the simulation results show that the hydraulic fracture evolution is influenced by the preexisting fractures; the induced stresses are larger in the direction perpendicular to the fracture; and in the same direction, the induced stresses are larger when the distance to the preexisting fracture is smaller. The simulation results are supported by mechanics analysis. The hydraulic fractures propagate along the large induced stresses and perpendicularly contact the artificial structural planes, forming the two-step fracturing fracture network. This 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. This 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.

8. Conclusions

(i)The pore pressure value and the nonuniform pore pressure gradient directly influence the tensile fracture initiation and propagation in rock. The 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. The asymmetrical pressure gradient can also influence the fracture propagation direction. The pore pressure gradient is the main controlling factor of model stability. The 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. The 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

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

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This study was supported by the Fund of Education Department of Liaoning Province (Grant no. 2017LNQN12).

References

D. R. Suchy and K. D. Newell, “Hydraulic fracturing of oil and gas wells in Kansas,” Kansas Geological Survey Public Information Circular, vol. 32, pp. 1–6, 2012.

M. F. M. Sena, L. P. Rosa, and A. Szklo, “Will venezuelan extra-heavy oil be a significant source of petroleum in the next decades?” Energy Policy, vol. 61, pp. 51–59, 2013.

C. Jia, “Breakthrough and significance of unconventional oil and gas to classical petroleum geology theory,” Petroleum Exploration and Development, vol. 44, no. 1, pp. 1–10, 2017.

Z. T. Li, J. Z. Su, and R. S. Zhang, Modern Fracturing Technology of Horizontal Well in Shale, Sinopec Press, Beijing, China, 2015.

C. C. Bose, B. Fairchild, T. Jones, A. Gul, and R. B. Ghahfarokhi, “Application of nanoproppants for fracture conductivity improvement by reducing fluid loss and packing of micro-fractures,” Journal of Natural Gas Science and Engineering, vol. 27, pp. 424–431, 2015.

Q. Wen, S. Wang, X. Duan et al., “Experimental investigation of proppant settling in complex hydraulic-natural fracture system in shale reservoirs,” Journal of Natural Gas Science and Engineering, vol. 33, pp. 70–80, 2016.

H. Hofmann, T. Babadagli, and G. Zimmermann, “Numerical simulation of complex fracture network development by hydraulic fracturing in naturally fractured ultratight formations,” Journal of Energy Resources Technology, vol. 136, 2014.

J. A. Breyer, “Shale reservoirs: giant resources for the 21st century,” American Association of Petroleum Geologists, pp. 1–20, 2012.

T. Yan, W. Li, and X. Bi, “An experimental study of fracture initiation mechanisms during hydraulic fracturing,” Petroleum Science, vol. 8, no. 1, pp. 87–92, 2011.

J. He, C. Lin, X. Li, and X. Wan, “Experimental investigation of crack extension patterns in hydraulic fracturing with shale, sandstone and granite cores,” Energies, vol. 9, no. 12, p. 1018, 2016.

B. Hou, M. Chen, Z. M. Li, Y. H. Wang, and C. Dao, “Propagation area evaluation of hydraulic fracture networks in shale gas reservoirs,” Petroleum Exploration and Development, vol. 41, no. 6, pp. 763–768, 2014.

Y. Tang, X. Tang, G. Y. Wang, and Q. Zhang, “Summary of hydraulic fracturing technology in shale gas development,” Geological Bulletin of China, vol. 30, no. 2, pp. 393–399, 2011.

M. S. Bruno and F. M. Nakagawa, “Pore pressure influence on tensile fracture propagation in sedimentary rock,” International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, vol. 28, no. 4, pp. 261–273, 1991.

J. Yao, Q. Zeng, Z. Huang, H. Sun, and L. Zhang, “Numerical modeling of simultaneous hydraulic fracturing in the mode of multi-well pads,” Science China Technological Sciences, vol. 60, no. 2, pp. 232–242, 2017.

V. Sesetty and A. Ghassemi, “A numerical study of sequential and simultaneous hydraulic fracturing in single and multi-lateral horizontal wells,” Journal of Petroleum Science and Engineering, vol. 132, pp. 65–76, 2015.

L. C. Li, Y. J. Xia, and C. A. Tang, “Investigation on the effect of pore pressure gradient on fracture propagation in rock materials,” Materials Research Innovations, vol. 19, pp. 237–240, 2015.

O. Pourhosseini and M. Shabanimashcool, “Development of an elasto-plastic constitutive model for intact rocks,” International Journal of Rock Mechanics and Mining Sciences, vol. 66, pp. 1–12, 2014.

X. Tan, H. Konietzky, and T. Fruhwirt, “Numerical simulation of triaxial compression test for brittle rock sample using a modified constitutive law considering degradation and dilation behavior,” Journal of Central South University, vol. 22, pp. 3091–3107, 2015.

K. Ma, C. A. Tang, L. C. Li, P. G. Ranjith, M. Cai, and N. W. Xu, “3D modeling of stratified and irregularly jointed rock slope and its progressive failure,” Arabian Journal of Geosciences, vol. 6, no. 6, pp. 2147–2163, 2013.

N. W. Xu, C. A. Tang, H. Li et al., “Excavation-induced microseismicity: microseismic monitoring and numerical simulation,” Journal of Zhejiang University Science A, vol. 13, no. 6, pp. 445–460, 2012.

T. Xu, Q. Xu, M. Deng, and T. Ma, “A numerical analysis of rock creep-induced slide: a case study from jiweishan mountain, China,” Environmental Earth Sciences, vol. 72, no. 6, pp. 2111–2128, 2014.

L. C. Li, T. H. Yang, and Z. Z. Liang, “Numerical investigation of groundwater outbursts near faults in underground coal mines,” International Journal of Coal Geology, vol. 85, pp. 276–288, 2011.

S. Y. Wang, S. W. Sloan, C. A. Tang, and W. C. Zhu, “Numerical simulation of the failure mechanism of circular tunnels in transversely isotropic rock masses,” Tunnelling and Underground Space Technology, vol. 32, pp. 231–244, 2012.

X. Chen, J. Yu, C. A. Tang, H. Li, and S. Wang, “Experimental and numerical investigation of permeability evolution with damage of sandstone under triaxial compression,” Rock Mechanics and Rock Engineering, vol. 50, no. 6, pp. 1529–1549, 2017.

T. Xu, C. A. Tang, and J. Zhao, “Modeling of rheological deformation of inhomogeneous rock and associated time-dependent response of tunnels,” International Journal of Geomechanics, vol. 12, no. 2, pp. 147–159, 2012.

J. Liu and E. Y. Wang, “Effects of coal-rock strength on mechanical behavior and failure mechanism of samples composed of coal and rock,” Disaster Advances, vol. 6, pp. 57–68, 2013.

Q. He, F. T. Suorineni, T. Ma, and J. Oh, “Effect of discontinuity stress shadows on hydraulic fracture re-orientation,” International Journal of Rock Mechanics and Mining Sciences, vol. 91, pp. 179–194, 2017.

S. Q. Yang, T. Xu, L. He, H. W. Jing, S. Wen, and Q. L. Yu, “Numerical study on failure behavior of brittle rock specimen containing pre-existing combined flaws under different confining pressure,” Archives of Civil and Mechanical Engineering, vol. 15, no. 4, pp. 1085–1097, 2015.

X. X. Men, J. R. Li, and Z. H. Han, “Fracture propagation behavior of jointed rocks in hydraulic fracturing,” Advances in Materials Science and Engineering, vol. 2018, Article ID 9461284, 6 pages, 2018.

L. C. Li, Y. J. Xia, B. Huang, L. Y. Zhang, M. Li, and A. S. Li, “The behaviour of fracture growth in sedimentary rocks: a numerical study based on hydraulic fracturing processes,” Energies, vol. 9, 2016.

C. A. Tang, H. Liu, P. K. K. Lee, Y. Tsui, and L. G. Tham, “Numerical studies of the influence of microstructure on rock failure in uniaxial compression—part I: effect of heterogeneity,” International Journal of Rock Mechanics and Mining Sciences, vol. 37, no. 4, pp. 555–569, 2000.

T. F. Wong, R. H. C. Wong, K. T. Chau, and C. A. Tang, “Microcrack statistics, weibull distribution and micromechanical modeling of compressive failure in rock,” Mechanics of Materials, vol. 38, no. 7, pp. 664–681, 2006.

J. T. Shakib, A. Ghaderi, and A. A. Shahri, “Analysis of hydraulic fracturing in fractured reservoir: interaction between hydraulic fracture and natural fractures,” Life Science Journal-Acta Zhengzhou University Overseas Edition, vol. 9, pp. 1854–1862, 2012.

X. S. Shang, Y. H. Ding, L. F. Yang, and Y. J. Lu, “Network fracturing technique in shale reservoirs based on weak discontinuity and fracture interaction,” Natural Gas Geoscience, vol. 10, pp. 1883–1891, 2016.

A. E. Green and I. N. Sneddon, “The distribution of stress in the neighborhood of a crack in an elastic solid,” Proceedings of the Royal Society of London, vol. 187, no. 1009, pp. 229–260, 1946.

T. H. Yang, “Study on infiltrate character and coupling analysis of seepage and stress in rock failure process,” Northeastern University, Shenyang, China, 2001, Ph.D. thesis.