**1. Introduction**

With the rapid growth of the global economy, the conventional oil and gas resources with decreasing production cannot meet the need for energy. Unconventional energy resources such as shale gas, shale oil, gas hydrate, and dry hot rock geothermal energy are quite abundant, and their economically and technically feasible exploitation is an e ffective solution to the problem of future

energy shortages [1–4]. However, extremely tight pore structure and ultra-low permeability for the unconventional reservoir rocks restrict the flow of energy and e fficient production [5,6]. Hydraulic fracturing is an important technical method for reservoir reconstruction in the process of energy resources exploitation. By injecting high-pressure fluid into the reservoir, the rock is fractured, thereby increasing the conductivity [7–9]. A ffected by the heterogeneity and anisotropy of the reservoir rock, the initiation and propagation of fractures, the leak-o ff of the fracturing fluid, the interaction between hydraulic fractures and natural fractures, and many other physical and mechanical conditions, hydraulic fracturing is generally regarded as the multi-field coupling problem under complex in-situ stress [10–13]. Mastering the mechanism of hydraulic fracturing technology is crucial to building a complex fracture network system and forming a high-yield unconventional energy reservoir.

Laboratory testing is a direct and e ffective method to discover the mechanism of hydraulic fracturing. Various observation and analysis techniques have been utilized by domestic and foreign geologists to study the fracturing behaviors of reservoir rocks under the action of pressure fluids [14–17]. Guo et al. [18] used the true triaxial test system and the high-energy computed tomography (CT) scanning based on linear accelerator to capture the e ffects of flow rate, in-situ stress and fluid viscosity on fracturing horizontal well in shale. They indicated that the abovementioned multi-aspect factors dominate the initiating and propagating rules of hydraulic fractures, and the interaction between hydraulic fractures and sedimentary bedding. Zhang et al. [19] observed the propagation patterns of supercritical carbon dioxide (SC-CO2) induced hydraulic fractures in shale through CT scanning and digital radiography (DR) scanning. Comparing with hydraulic fracturing, they believed that the percolation e ffect of SC-CO2 reduces the initiation pressure needed to fractures by more than 50%, and induces more secondary fractures to form complex fracture networks. He et al. [20] investigated the surface characteristics of hydraulic fractures using the scanning electron microscope (SEM), and found that numerous micro cracks propagating along the boundaries of hard minerals and soft organic matters connected the pores and macro hydraulic fractures, thus increasing the permeability of shale. Zhuang et al. [21] performed a set of hydraulic fracturing tests on granites under true triaxial loading with di fferent injection schemes, and analyzed the mineral fracturing behaviors and induced seismicity according to thin section microscopic observations and acoustic emission technology. The results show that hydraulic fractures under six fluid injection schemes are all composed of the intragranular cracks splitting microcline, orthoclase and quartz, but the cyclic pulse pressurization scheme improves the injectivity and decreases induced seismicity.

In addition to these laboratory tests, low-cost and high-e fficiency numerical simulation becomes one of the most common ways to solve the hydraulic fracturing problems as the computer science develops [22]. Various numerical simulation methods including the finite-element method (FEM), the finite-di fference method (FDM), the discrete element method (DEM), and the discontinuous displacement method (DDM) are successful in reproducing the hydraulic fracturing process of reservoir rocks and revealing the mechanism of hydraulic fracture propagation [23,24]. Li et al. [25] simulated the multi-cluster hydraulic fracturing based on the FEM software ABAQUS, and demonstrated that stress interference from multiple-clusters causes the suppression and transfer of the fracture network. Liao et al. [26] stated that horizontal stress anisotropy plays a key role in the interaction between hydraulic fractures and natural fractures in term of the three-dimensional fracturing simulation by the FDM software FLAC 3D. Zangeneh et al. [27] established a conceptual reservoir model with faults in the DEM software UDEC to confirm the influence of fluid di ffusion on triggering fault slip and to evaluate the maximum magnitude of seismic events. Janiszewski et al. [28] discussed the interaction of fractures in granite under di fferent approach angle conditions with the help of the DDM software FRACOD, and concluded that the existence of low dip-angle fractures allows more complex fracture networks.

Each kind of numerical simulation software developed on the basis of the continuous or discontinuous methods has its own unique advantages. As a typical representative of DEM software, the numerical model created in PFC is composed of thousands of basic circular particles bonded together by the contact bonds [29]. Furthermore, the failure of a contact bond corresponds to the generation of a microcrack in the model. Because of the simple topology and the fast updating mechanism of particle information, PFC shows obvious advantages in emulating the initiation, propagation and aggregation of microcracks and the ultimate macro failure for rock and soil mass under complex stress conditions [30–33].

The fluid-mechanical coupling algorithm in PFC was first proposed by Cundall [34], which not only considers the interaction process of pore pressure and fracture aperture, but also realizes the fluid flow in both fractures and rock matrix. Al-Busaidi et al. [35] recreated the laboratory test results of hydraulic fracturing for Lac du Bonnet granite by simulation, which proved the e ffectiveness of this algorithm. Yoon et al. [36–38] probed into the seismic activity induced by fluid injection during the development of deep granite reservoir. The detailed seismic parameters involving temporal and spatial distribution, moment magnitudes and radiated seismic energy of the induced seismic events were obtained through the reservoir fracturing simulations, which provided constructive guidance for the optimization of field fracturing schemes. Tomac et al. [39,40] extended this algorithm and explored the influences of mineral grains, fluid pressure and temperature on permeable hot granite during injecting pressurized cold fluid. They presented that the thermal damage area of wellbore caused by rock cooling and cold fluid infiltration increases as the fluid dynamic viscosity decreases. In our previous research, the calculation method of fluid pressure in the domain after contact bond failure was modified to make the update of the simulation results more accurate when the numerical model contains a large-volume injection well [41,42]. Following the improved algorithm, we comprehensively studied the hydraulic fracturing mechanism in isotropic and laminated reservoirs under di fferent in-situ stress state, injection rate and fluid viscosity, and validated the ability in modeling the variety of the interaction scenarios between hydraulic fractures and natural fractures such as direct crossing and crossing with an o ffset [43–47].

Although the PFC has been employed in many hydraulic fracturing studies, most of them only consider the simple case of single injection well and ignore the interference of the multiple injection wells. However, in order to obtain a better fracturing e ffect in the field of production, it is necessary to arrange multiple injection wells in a suitable space and select a reasonable fracturing sequence. In this study, taking the Alxa porphyritic granite as an example, a two-well fracturing model based on our modified algorithm was established to discuss the e ffects of injection sequence and well spacing on hydraulic fracturing behaviors. The rest of this paper is organized as follows. First, the characteristics of contact bond and the modified fluid-mechanical coupling algorithm in PFC were introduced. Second, the micromechanical parameters were carefully calibrated and validated according to the experimental results of uniaxial loading and the analytical solutions of the breakdown pressure, respectively. Third, the boundary conditions of the numerical model were transformed with reference to diverse schemes so as to analyze the change in breakdown pressure, fracture propagation pressure, fluid pressure distribution and fracture propagation pattern. Finally, the stress variation along the height direction of hydraulic fracture was extracted to assess the stress shadow e ffect of hydraulic fractures.

## **2. Modeling Methodology**

#### *2.1. Introduction of Particle Flow Code*

PFC is a commercial DEM software that applies the rigid, unbreakable circular particles as the basic elements of the model. The circular particles in contact are available for connection as a whole by giving a contact bond model. During the process of solution, the displacements, contact forces and moments of particles are continuously updated with the calculation time by Newton's second law and force-displacement law until the unbalanced force and unbalanced moment in the model reach the program termination conditions [29].

Plentiful models of contact bonds built into PFC including the Parallel-Bond Model (PBM), the Smooth-Joint Model (SJM), the Flat-Joint Model (FJM) and the Hertz Contact Model (HCM) determine the macro mechanical properties of various solid materials [48]. In particular, envisioned

as a collection of springs with constant strength and sti ffness (Figure 1a), the PBM is mostly exerted into the simulations of rock materials by assigning the particle-interaction laws. Moreover, the normal component *Fn i* (*t* + Δ*t*) and the tangential component *Fs i* (*t* + Δ*t*) for the contact force and the contact moment *M*(*t* + Δ*t*) at time *t* + Δ*t* depend on the corresponding values at the start of the timestep with the increments. The detailed equations used are given by [43]:

$$F\_i^n(t + \Delta t) = F\_i^n(t) + \Delta F\_i^n = F^n(t)n\_i + \left(-k^n A \Delta t I\_i^n\right)n\_i \tag{1}$$

$$F\_i^\mathfrak{s}(t + \Delta t) = F\_i^\mathfrak{s}(t) + \Delta F\_i^\mathfrak{s} = F\_i^\mathfrak{s}(t) + \left(-k^\mathfrak{s} A \Delta l I\_i^\mathfrak{s}\right) \tag{2}$$

$$M(t + \Delta t) = M(t) + \Delta M = M(t) + (-k^n I)(\omega^{[B]} - \omega^{[A]})\Delta t \tag{3}$$

where *kn* and *ks* are the normal and shear sti ffness of the parallel bond, respectively; Δ *Un* and Δ *Us* are the relative normal-displacement increment and the relative shear-displacement increment in one timestep Δ*t*, respectively; ω[*A*] and ω[*B*] are the relative rotation velocities of the two contacting particles, respectively; *ni* is the normal vector of each contact (*i* = 1 or 2); *A* is the cross sectional area of the bond, and *I* is the moment of inertia of the bond.

**Figure 1.** Schematic diagram of the Parallel-Bond contact in PFC2D (modified from Itasca Consulting Group, Inc.: Minneapolis, MN, USA [48]): (**a**) conceptual models; (**b**) the deformation and failure mechanisms; and (**c**) the failure mode.

Then the normal stress σ and shear stress τ distributed on the cross-section of the bond are updated by the following equations [48]:

$$
\sigma = \frac{-F^n}{A} + \beta \frac{|M|}{I} R \tag{4}
$$

$$
\pi = \frac{\left| F\_i^s \right|}{A} \tag{5}
$$

where β is the moment-contribution factor.

At any time of solution, if the absolute value of normal stress |σ| inside the bond exceeds the tensile strength of the bond σc, or the absolute value of shear stress |τ| inside the bond exceeds the shear strength of the bond τc, the tensile failure or shear failure occurs respectively in the PBM (Figure 1b). Correspondingly, a tensile microcrack or a shear microcrack is generated in the numerical model. After the failure of the parallel bond, the force, moment, and sti ffnesses related to the contact are removed, resulting in the rotation of the two particles around each other under external force (Figure 1c). Recent studies have pointed out that the rotation after the bond failure reduces the self-locking e ffect between the particles, causing an extremely low ratio of uniaxial compressive strength (UCS) to uniaxial tensile strength (UTS) [49,50]. Nonetheless, the PBM revised by potyondy [51] and Ding and Zhang [52] is still employed in this study, in which the moment-contribution factor β is added to the bond failure criterion to overcome the above defects.

#### *2.2. Modified Fluid-Mechanical Coupling Algorithm*

Cundall's fluid-mechanical coupling algorithm achieves the flow of viscous fluid in the parallel bond. As illustrated in Figure 2a, the fluid network topology is constructed after the circular particles are assembled by the PBM. More specifically, connecting the centers of the contacting particles by blue lines forms a series of enclosed domains (blue polygons) which are regarded as "reservoirs" to store the fluid pressure. The centers of the adjacent reservoir domains (blue circles) are linked by the virtual pipes (magenta lines) to ensure fluid flow. Any pipe is approximated as two parallel plates existing in the parallel bond. A di fferential pressure between the two connected reservoir domains drives the flow of fluid from the high-pressure region to the low-pressure region. The volumetric flow rate *q* is calculated by the cubic law [34]:

$$q = \frac{e^3}{12\mu} \frac{P\_2 - P\_1}{L} \tag{6}$$

where *e* is the hydraulic aperture; *P*2 − *P*1 is the fluid pressure di fference between the reservoir domains; *L* represents the length of the flow channel; μ is the fluid dynamic viscosity.

Under the action of confining pressure and fluid pressure, the hydraulic aperture *e* changes with the normal stress σ of the bond, which is described as [53,54]:

$$
\epsilon = c\_{\rm inf} + (c\_0 - c\_{\rm inf}) \exp(-0.15\sigma) \tag{7}
$$

when σ tends to infinity, the aperture decreases gradually to *e*inf; when σ = 0, the aperture *e*0 represents the residual aperture due to no force interaction between the particles on both sides of the flow channel. Therefore, this mechanism enables the fluid to di ffuse into the surrounding reservoirs whether fractures are in the model. According to the Equation (8), the values of *e*0 and *e*inf are determined by the permeability *k*, total volume *V* of the reservoirs and flow channel length *L* [35].

$$k = \frac{1}{12V} \sum\_{\text{pips}} Lx^3 \tag{8}$$

In the fluid calculation process of one timestep, the increment of fluid pressure Δ *P* for a reservoir domain is derived from the apparent volume change of this domain Δ *Vd* and the fluid flow of surrounding reservoir domains Σ*q*, as shown in Equation (9) [43]:

$$
\Delta P = \frac{K\_f}{V\_d} \left( \sum q \Delta t - \Delta V\_d \right) \tag{9}
$$

where *Kf* is bulk modulus of fluid; *Vd* is the apparent volume of the reservoir domain.

After the fluid calculation step, the fluid pressure acts on the particles in the form of the body force (Figure 2b), and the magnitude of the force is the product of the fluid pressure *Pf* and the length of the connecting line *li* at the contact point. As the parallel bond between adjacent domains is broken owing to the effect of confining pressure and fluid pressure, Al-Busaidi et al. [35] and Zhao and Young [55] used the average value of the fluid pressures *Pf*1 and *Pf*2 in the two domains before bond failure to express the fluid pressure *P f*in the new domain:

$$P\_f' = \frac{P\_{f1} + P\_{f2}}{2} \tag{10}$$

Obviously, the Equation (10) is not suitable for the cases of large-volume injection well in the model because of the inaccurate average fluid pressure in the domains around the well. An optimized computing method is as follows [41]:

$$P\_f' = \left[\frac{(V\_{f1} + V\_{f2})}{(V\_{r1} + V\_{r2})q} - 1\right] \mathbf{K}\_f \tag{11}$$

where *Vf*1 and *Vf*2 are the fluid volume in the two domains before bond failure, respectively; *Vr*1 and *Vr*2 are the volume of the two domains under the fluid pressure of 0 MPa; ϕ is the porosity of the model.

**Figure 2.** The sketch of the fluid-mechanical coupling model in PFC2D: (**a**) the generation of fluid network; and (**b**) fluid-particle interaction mechanism (yellow circles: solid particles, magenta lines: flow channels, blue polygons: reservoir domains, blue circles: the centers of reservoirs, red lines: Parallel-Bond contacts).

#### **3. Modeling Calibration and Validation**

The input micro-parameters involved in the fluid-mechanical coupling model in PFC2D are divided into three categories to characterize the physical and mechanical properties of basic particles, parallel bonds and fluid networks, most of which cannot be directly obtained from the laboratory tests. Debugging the parameters to match the calibration results with experimental results or analytical solutions is absolutely essential for the subsequent reliable numerical simulations. Fortunately, the relationships between the input parameters and macro mechanical characteristics of the models have been presented in several studies [56–59], which is helpful for the quick selection and calibration of the parameters. Combined with the uniaxial compression and tensile test results of Alxa porphyritic

granite, a group of input parameters which reproduced the mechanical behaviors of the real granite under uniaxial loading were determined in this section. To verify the accuracy of the parameters in simulating hydraulic fracturing, the granite model with a single injection well were created and contrasted with the analytical solutions of the breakdown pressure.
