Numerical Investigation of Injection-Induced Fracture Propagation in Brittle Rocks with Two Injection Wells by a Modiﬁed Fluid-Mechanical Coupling Model

: Hydraulic fracturing is a key technical means for stimulating tight and low permeability reservoirs to improve the production, which is widely employed in the development of unconventional energy resources, including shale gas, shale oil, gas hydrate, and dry hot rock. Although signiﬁcant progress has been made in the simulation of fracturing a single well using two-dimensional Particle Flow Code (PFC2D), the understanding of the multi-well hydraulic fracturing characteristics is still limited. Exploring the mechanisms of ﬂuid-driven fracture initiation, propagation and interaction under multi-well fracturing conditions is of great theoretical signiﬁcance for creating complex fracture networks in the reservoir. In this study, a series of two-well fracturing simulations by a modiﬁed ﬂuid-mechanical coupling algorithm were conducted to systematically investigate the e ﬀ ects of injection sequence and well spacing on breakdown pressure, fracture propagation and stress shadow. The results show that both injection sequence and well spacing make little di ﬀ erence on breakdown pressure but have huge impacts on fracture propagation pressure. Especially under hydrostatic pressure conditions, simultaneous injection and small well spacing increase the pore pressure between two injection wells and reduce the e ﬀ ective stress of rock to achieve lower fracture propagation pressure. The injection sequence can change the propagation direction of hydraulic fractures. When the in-situ stress is hydrostatic pressure, simultaneous injection compels the fractures to deﬂect and tend to propagate horizontally, which promotes the formation of complex fracture networks between two injection wells. When the maximum in-situ stress is in the horizontal direction, asynchronous injection is more conducive to the parallel propagation of multiple hydraulic fractures. Nevertheless, excessively small or large well spacing reduces the number of fracture branches in fracture networks. In addition, the stress shadow e ﬀ ect is found to be sensitive to both injection sequence and well spacing.


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 effective solution to the problem of future

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 Energies 2020, 13, 4718 4 of 26 as a collection of springs with constant strength and stiffness (Figure 1a), the PBM is mostly exerted into the simulations of rock materials by assigning the particle-interaction laws. Moreover, the normal component F n i (t + ∆t) and the tangential component F s 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 n i (t + ∆t) = F n i (t) + ∆F n i = F n (t)n i + −k n A∆U n i n i (1) where k n and k s are the normal and shear stiffness of the parallel bond, respectively; ∆U n and ∆U s 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; n i 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. 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 stiffness (Figure 1a), the PBM is mostly exerted into the simulations of rock materials by assigning the particle-interaction laws. Moreover, the normal component F i n (t + ∆t) and the tangential component F i s 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]: where k n and k s are the normal and shear stiffness of the parallel bond, respectively; ΔU n and ΔU s 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.  Then the normal stress σ and shear stress τ distributed on the cross-section of the bond are updated by the following equations [48]: Energies 2020, 13, 4718 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 stiffnesses 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 effect 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.

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 differential 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]: where e is the hydraulic aperture; P 2 − P 1 is the fluid pressure difference 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]: 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 diffuse 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].
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 ∆V d and the fluid flow of surrounding reservoir domains Σq, as shown in Equation (9) [43]: where K f is bulk modulus of fluid; V d 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 P f and the length of the connecting line l i 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 P f1 and P f2 in the two domains before bond failure to express the fluid pressure P f in the new domain: 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]: where V f1 and V f2 are the fluid volume in the two domains before bond failure, respectively; V r 1 and V r 2 are the volume of the two domains under the fluid pressure of 0 MPa; ϕ is the porosity of the model. (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 Pf1 and Pf2 in the two domains before bond failure to express the fluid pressure P f ʹ in the new domain: 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]: where Vf1 and Vf2 are the fluid volume in the two domains before bond failure, respectively; Vr1 and Vr2 are the volume of the two domains under the fluid pressure of 0 MPa; φ is the porosity of the model.

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][57][58][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 Energies 2020, 13, 4718 7 of 26 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.

Laboratory Test
As one of the preferred candidate locations for geological disposal of high-level radioactive waste, Alxa is located in the west of Inner Mongolia Autonomous Region, with the geographic coordinates of 97 • 10 ~106 • 53 E and 37 • 24 ~42 • 47 N. The NRG-1 deep borehole in Alxa revealed that a large number of porphyritic granites are distributed from the surface to 600 m underground. The granitic cores from different depths of the borehole were processed into the standard specimens with a diameter of 50 mm and a height of 100 mm. Before the uniaxial compression tests, the average density of the specimens was 2650.0 kg/m 3 . During the testing, the MTS 815 servo-controlled hydraulic testing machine at the Institute of Crustal Dynamics, China Earthquake Administration was employed to compress the granites with an axial loading speed of 0.06 m/min.
The failure characteristics of Alxa porphyritic granites after the tests and the stress-strain curves are shown in Figures 3 and 4, respectively. It can be seen that all the specimens were cut by numerous splitting cracks, which conforms to the typical brittle fracturing pattern. The macro mechanical parameters of the specimens were extracted from the stress-strain curves, and the average values of the UCS, the elastic modulus (E) and the Poisson's ratio (ν) were 159.4 MPa, 48.0 GPa and 0.23, respectively. The tensile test is not conducted in this research; however, the average UTS of this kind of granites is about 9.5 MPa according to our previous studies [60,61].
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][57][58][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.

Laboratory Test
As one of the preferred candidate locations for geological disposal of high-level radioactive waste, Alxa is located in the west of Inner Mongolia Autonomous Region, with the geographic coordinates of 97°10′ ~ 106°53′ E and 37°24′ ~ 42°47′ N. The NRG-1 deep borehole in Alxa revealed that a large number of porphyritic granites are distributed from the surface to 600 m underground. The granitic cores from different depths of the borehole were processed into the standard specimens with a diameter of 50 mm and a height of 100 mm. Before the uniaxial compression tests, the average density of the specimens was 2650.0 kg/m 3 . During the testing, the MTS 815 servo-controlled hydraulic testing machine at the Institute of Crustal Dynamics, China Earthquake Administration was employed to compress the granites with an axial loading speed of 0.06 m/min.
The failure characteristics of Alxa porphyritic granites after the tests and the stress-strain curves are shown in Figures 3 and 4, respectively. It can be seen that all the specimens were cut by numerous splitting cracks, which conforms to the typical brittle fracturing pattern. The macro mechanical parameters of the specimens were extracted from the stress-strain curves, and the average values of the UCS, the elastic modulus (E) and the Poisson's ratio (ν) were 159.4 MPa, 48.0 GPa and 0.23, respectively. The tensile test is not conducted in this research; however, the average UTS of this kind of granites is about 9.5 MPa according to our previous studies [60,61].

Material Parameters Calibration
A rectangular numerical model with the same size as the test specimen was generated in PFC2D to emulate the uniaxial loading responses. The model was comprised of 10,412 circular particles with a radius distribution of 0.3-0.45 mm and was assigned the contact bonds of PBM. The effective

Material Parameters Calibration
A rectangular numerical model with the same size as the test specimen was generated in PFC2D to emulate the uniaxial loading responses. The model was comprised of 10,412 circular particles with a radius distribution of 0.3-0.45 mm and was assigned the contact bonds of PBM. The effective modulus and stiffness of particles can inherit from the corresponding parameters of PBM. After the initial parameters were given, the simulations of uniaxial compression and direct tension were carried out to calculate the UCS, E, ν and UTS. Then the relevant parameters were continuously adjusted by the comparison between the experimental results and the simulated results. The final calibrated parameters are listed in Table 1, by which the simulated macro mechanical properties are in good agreement with the actual values, as shown in Table 2. The stress-strain curves and the spatial distribution of cracks are presented in Figure 4, which also proves the accuracy of parametric calibration. carried out to calculate the UCS, E, ν and UTS. Then the relevant parameters were continuously adjusted by the comparison between the experimental results and the simulated results. The final calibrated parameters are listed in Table 1, by which the simulated macro mechanical properties are in good agreement with the actual values, as shown in Table 2. The stress-strain curves and the spatial distribution of cracks are presented in Figure 4, which also proves the accuracy of parametric calibration.  Most of the microcracks generated in the numerical granite under the compressive loading are the tension-type modes. Fewer microcracks formed in the middle of model under the tensile loading, and aggregated into a macro fracture.

Micro-Parameters of the Basic Particles Notations Values
Particle density (kg/m 3 ) 2650.0 Porosity (%) ϕ 10 Effective modulus of particle (GPa) E' 30.0 Contact normal to shear stiffness ratio k n /k s 1.5 Friction coefficient of particle f 0.8

Micro-Parameters of the Parallel Bonds Notations Values
Effective modulus of the parallel bonds (GPa) E c 30.0 Normal to shear stiffness ratio of the parallel bonds k n /k s 1.

Validation of Hydraulic Fracturing Model
The fluid parameters for hydraulic fracturing were decided by the actual working conditions. After constructing the fluid networks and adding the fluid parameters in the PBM, it was still vital to verify whether the fracturing simulation could describe the realistic fracturing behaviors, including the breakdown pressure and the propagation of hydraulic fractures. A classical equation for the breakdown pressure P wf when the fluid is injected into a circular elastic borehole was proposed by Haimson and Fairhurst [62] as following expression: where σ v is the confining pressure in vertical direction; σ h is the confining pressure in horizontal direction; P 0 is initial pore pressure; σ t is the UTS of the model. In order to facilitate the comparison of simulated and theoretical breakdown pressures and avoid the stress interference of multiple injection wells, the granite model with a single injection well was established, as shown in Figure 5a. The square fracturing model with a width and height of 1.0 m consisted of 11,778 circular particles with a radius distribution of 4-6 mm. Furthermore, the vertical and horizontal confining pressures were servo-controlled by moving the walls around the model. An injection well with a radius of 50 mm was excavated in the center of the model to inject viscous fluid. It is noted that the particles surrounding the injection well have been replaced by smaller particles to smooth the well surface and prevent stress concentration (Figure 5b). No fluid flow domains were covered on the border of the fracturing model to give the impermeable boundary conditions. The hydraulic apertures e 0 and e inf calculated by Equations (7) and (8) were 2.2 × 10 −6 and 2.2 × 10 −7 , respectively, which were input into the model as basic hydraulic property parameters ( Table 1). The fracturing fluid was the liquid water with a bulk modulus of 2.0 GPa and a viscosity of 1.0 × 10 −3 . The numerical model kept the vertical confining pressure σ v and the horizontal confining pressure σ h equal. Under the conditions of the initial pore pressure of 0 MPa and the fluid injection rate of 1.0 × 10 −5 m 3 /s, the hydraulic fracturing simulations with confining pressures from 5 MPa to 50 MPa (5-MPa intervals) were implemented.
Energies 2020, 13, x FOR PEER REVIEW 10 of 26  The simulated values and theoretical values of breakdown fractures are summarized in Figure  6, and the error between them is less than 20%. The main reason for the deviation is that the analytical model assumes that the entire reservoir is impermeable, while the modified hydraulic fracturing model in PFC considers the leakage of fracturing fluid to the surrounding domains.  The simulated values and theoretical values of breakdown fractures are summarized in Figure 6, and the error between them is less than 20%. The main reason for the deviation is that the analytical model assumes that the entire reservoir is impermeable, while the modified hydraulic fracturing model in PFC considers the leakage of fracturing fluid to the surrounding domains.
Energies 2020, 13, x FOR PEER REVIEW 10 of 26  The simulated values and theoretical values of breakdown fractures are summarized in Figure  6, and the error between them is less than 20%. The main reason for the deviation is that the analytical model assumes that the entire reservoir is impermeable, while the modified hydraulic fracturing model in PFC considers the leakage of fracturing fluid to the surrounding domains.    Figure 7 shows the representative borehole pressure histories, crack spatial distribution and corresponding fluid pressure field in hydraulic fracturing under different confining pressures. It is found that the borehole pressure increases rapidly at the beginning of injection, and as the borehole pressure approaches the breakdown pressure, the borehole pressure increases slowly and non-linearly, which reflects that the leakage of fracturing faster is greater when the borehole pressure is high. The leakage of fluid from the injection well and fractures to the surrounding fluid networks can also be seen from the figures of the fluid pressure distribution.
Energies 2020, 13, x FOR PEER REVIEW 11 of 26 Figure 7 shows the representative borehole pressure histories, crack spatial distribution and corresponding fluid pressure field in hydraulic fracturing under different confining pressures. It is found that the borehole pressure increases rapidly at the beginning of injection, and as the borehole pressure approaches the breakdown pressure, the borehole pressure increases slowly and nonlinearly, which reflects that the leakage of fracturing faster is greater when the borehole pressure is high. The leakage of fluid from the injection well and fractures to the surrounding fluid networks can also be seen from the figures of the fluid pressure distribution.  Before the borehole pressure reaches the breakdown pressure, a few cracks initiated around the wellbore. After increasing to the breakdown pressure, the borehole pressure dropped sharply to the fracture propagation pressure and remained unchanged. As a result, the number of microcracks grew dramatically along certain directions. Hubbert and Willis [63] proposed that the initiation and propagation of deep hydraulic fracturing fractures follow the direction of maximum principal stress. From Figure 7, we can see that the macro hydraulic fractures propagated along the different Before the borehole pressure reaches the breakdown pressure, a few cracks initiated around the wellbore. After increasing to the breakdown pressure, the borehole pressure dropped sharply to the fracture propagation pressure and remained unchanged. As a result, the number of microcracks grew dramatically along certain directions. Hubbert and Willis [63] proposed that the initiation and propagation of deep hydraulic fracturing fractures follow the direction of maximum principal stress.
From Figure 7, we can see that the macro hydraulic fractures propagated along the different directions under various confining pressures. The difficulty in predicting the propagation direction stems from the fact that the micro defects on the surface of the injection well have more significant impact on the initiation and propagation of hydraulic fractures under hydrostatic pressure conditions.
Another noteworthy issue is that almost all the hydraulic fractures acquired from the fluid-mechanical coupling algorithm are made up of tensile cracks, which seems to be distinct from lots of shear-type seismic events observed in hydraulic fracturing tests for granite [64,65]. Al-Busaidi et al. [35] have demonstrated that the shear failures recorded in the laboratory tests are caused by the slippage of the preexisting cracks in the specimens, and the mechanism of injection-induced fracture is predominantly tensile failure. In brief, the simulation results of granite fracturing using the modified fluid-mechanical coupling algorithm are accurate.

Modeling Scenarios
The multi-well hydraulic fracturing behaviors of reservoir rocks not only affected by the influencing factors of single-well fracturing (e.g., the mechanical properties of rocks, fluid injection rate, fluid viscosity and in-situ stress), but also by the influencing factors unique to multi-well fracturing (e.g., the number of injection wells, injection sequence and well spacing). Taking the two-well hydraulic fracturing process as the research object, a larger numerical model with a width of 1.5 m and a height of 1.0 m, which was an aggregation of 17,700 circular particles, was constructed based on the fluid-solid coupling method introduced in Section 2 to get enough space to propagate for the fractures (Figure 8a). Two injection wells (Well No.1 and Well No.2) with the radius of 50 mm and the well spacing of L are symmetrical about the central axis of the model. After the calibration parameters in Table 1 were given in the two-well fracturing model, the influences of injection sequence and well spacing under different in-situ stress states on the breakdown pressure, fluid pressure distribution and hydraulic fracture propagation were discussed. directions under various confining pressures. The difficulty in predicting the propagation direction stems from the fact that the micro defects on the surface of the injection well have more significant impact on the initiation and propagation of hydraulic fractures under hydrostatic pressure conditions. Another noteworthy issue is that almost all the hydraulic fractures acquired from the fluidmechanical coupling algorithm are made up of tensile cracks, which seems to be distinct from lots of shear-type seismic events observed in hydraulic fracturing tests for granite [64,65]. Al-Busaidi et al. [35] have demonstrated that the shear failures recorded in the laboratory tests are caused by the slippage of the preexisting cracks in the specimens, and the mechanism of injection-induced fracture is predominantly tensile failure. In brief, the simulation results of granite fracturing using the modified fluid-mechanical coupling algorithm are accurate.

Modeling Scenarios
The multi-well hydraulic fracturing behaviors of reservoir rocks not only affected by the influencing factors of single-well fracturing (e.g., the mechanical properties of rocks, fluid injection rate, fluid viscosity and in-situ stress), but also by the influencing factors unique to multi-well fracturing (e.g., the number of injection wells, injection sequence and well spacing). Taking the twowell hydraulic fracturing process as the research object, a larger numerical model with a width of 1.5 m and a height of 1.0 m, which was an aggregation of 17,700 circular particles, was constructed based on the fluid-solid coupling method introduced in Section 2 to get enough space to propagate for the fractures (Figure 8a). Two injection wells (Well No.1 and Well No.2) with the radius of 50 mm and the well spacing of L are symmetrical about the central axis of the model. After the calibration parameters in Table 1 were given in the two-well fracturing model, the influences of injection sequence and well spacing under different in-situ stress states on the breakdown pressure, fluid pressure distribution and hydraulic fracture propagation were discussed.

The Influence of Injection Sequence under Different In-Situ Stress Conditions
The well spacing L of the two-well fracturing model was set to 0.5 m. Three injection sequence schemes were designed for two water injection wells, which are: Scheme 1-simultaneous injection in Well No.

The Influence of Injection Sequence under Different In-Situ Stress Conditions
The well spacing L of the two-well fracturing model was set to 0.5 m. Three injection sequence schemes were designed for two water injection wells, which are: Scheme 1-simultaneous injection in and Scheme 3-first injection in Well No.2 and second injection in Well No.1. The fluid injection rate for each well was always 1.0 × 10 −5 m 3 /s. If any hydraulic fracture propagates to the boundary of the model, the fluid injection and numerical simulation stop immediately, regardless of the shut-in stage in the field fracturing process. In the asynchronous fracturing of two injection wells, the second injection was conducted after the pore fluid pressure produced by the first injection was removed, so as to eliminate the effect of residual pore pressure.
The influence of different in-situ stresses was also considered in these injection schemes. When the model was under hydrostatic pressure, the vertical and horizontal confining pressures were fixed at 30 MPa; when the maximum in-situ stress was in the horizontal direction, the vertical and horizontal confining pressures were fixed at 15 MPa and 30 MPa, respectively. The three injection schemes have been repeatedly executed under each in-situ stress condition.  (Figure 9). Similar to the fracturing of the single well, microcracks initially originated at the micro defects around the boreholes, and increased rapidly as the borehole pressures reduced to the fracture propagation pressure. However, the hydraulic fractures between the two wells were deflected at a large angle and close to each other during the propagations along the initial microcrack directions. Especially from the corresponding distribution of pore pressure field, it can be seen that the pore pressure field in the surrounding rock resulted from the fluid seepage in the two fractures will be superimposed on the tips. The other two hydraulic fractures were also considerably distorted until they propagated to the calculation boundary. In the end, four fracture branches formed in the model. for each well was always 1.0 × 10 −5 m 3 /s. If any hydraulic fracture propagates to the boundary of the model, the fluid injection and numerical simulation stop immediately, regardless of the shut-in stage in the field fracturing process. In the asynchronous fracturing of two injection wells, the second injection was conducted after the pore fluid pressure produced by the first injection was removed, so as to eliminate the effect of residual pore pressure. The influence of different in-situ stresses was also considered in these injection schemes. When the model was under hydrostatic pressure, the vertical and horizontal confining pressures were fixed at 30 MPa; when the maximum in-situ stress was in the horizontal direction, the vertical and horizontal confining pressures were fixed at 15 MPa and 30 MPa, respectively. The three injection schemes have been repeatedly executed under each in-situ stress condition.  (Figure 9). Similar to the fracturing of the single well, microcracks initially originated at the micro defects around the boreholes, and increased rapidly as the borehole pressures reduced to the fracture propagation pressure. However, the hydraulic fractures between the two wells were deflected at a large angle and close to each other during the propagations along the initial microcrack directions. Especially from the corresponding distribution of pore pressure field, it can be seen that the pore pressure field in the surrounding rock resulted from the fluid seepage in the two fractures will be superimposed on the tips. The other two hydraulic fractures were also considerably distorted until they propagated to the calculation boundary. In the end, four fracture branches formed in the model.   For the Scheme 2 of injection sequence, the breakdown pressures of Well No.1 and Well No.2 were 62.7 MPa and 58.8 MPa (Figure 10). After the first injection and the second injection, both hydraulic fractures in the left and right part of the model extended straight to the boundary at a high angle without distortion. Although the microcracks between the two injection wells were also initiated in the model, they did not propagate further.
For the Scheme 3 of injection sequence, the breakdown pressures of the two wells were 64.3 MPa and 58.3 MPa (Figure 11). The final fracture propagation pattern was the same as that in the Scheme 2. The above simulation results may indicate that under the hydrostatic pressure condition, obvious differences between simultaneous injection and asynchronous injection exist in the fracturing behaviors, but the injection sequences in asynchronous fracturing has little effect on the generation of cracks. For the Scheme 2 of injection sequence, the breakdown pressures of Well No.1 and Well No.2 were 62.7 MPa and 58.8 MPa (Figure 10). After the first injection and the second injection, both hydraulic fractures in the left and right part of the model extended straight to the boundary at a high angle without distortion. Although the microcracks between the two injection wells were also initiated in the model, they did not propagate further.
For the Scheme 3 of injection sequence, the breakdown pressures of the two wells were 64.3 MPa and 58.3 MPa (Figure 11). The final fracture propagation pattern was the same as that in the Scheme 2. The above simulation results may indicate that under the hydrostatic pressure condition, obvious differences between simultaneous injection and asynchronous injection exist in the fracturing behaviors, but the injection sequences in asynchronous fracturing has little effect on the generation of cracks.        In-situ differential stress gave rise to the propagations for the four hydraulic fractures along the horizontal direction in the simultaneous fracturing ( Figure 12). The fractures between the two wells were not at the identical height, yet had a tendency to deflect and merge.
After the fracturing of Well No.1 in the Scheme 2, two hydraulic fractures extended along the horizontal direction to both sides of the borehole (Figure 13). The propagation pattern of the left hydraulic fracture was entirely consistent with that in simultaneous fracturing. Interestingly, the right hydraulic fracture was not bent and extended straight through the bottom of Well No.2. From the borehole pressure history and pore pressure distribution, it seems that Well No.2 was isolated from the hydraulic fracture as before.  In-situ differential stress gave rise to the propagations for the four hydraulic fractures along the horizontal direction in the simultaneous fracturing ( Figure 12). The fractures between the two wells were not at the identical height, yet had a tendency to deflect and merge.
After the fracturing of Well No.1 in the Scheme 2, two hydraulic fractures extended along the horizontal direction to both sides of the borehole (Figure 13). The propagation pattern of the left hydraulic fracture was entirely consistent with that in simultaneous fracturing. Interestingly, the right hydraulic fracture was not bent and extended straight through the bottom of Well No.2. From the borehole pressure history and pore pressure distribution, it seems that Well No.2 was isolated from the hydraulic fracture as before.  The breakdown pressure, fracturing propagation pressure, crack number and hydraulic fracture branches after the fracturing under different in-situ stress states are listed in Table 3. In any injection sequence scheme, before the borehole pressure increased to the breakdown pressure, the fluid injected into the wells had no time to migrate in large area. Consequently, no evidence was found for the clear correlation between the breakdown pressure and injection sequence. As more fluid seeped The breakdown pressure, fracturing propagation pressure, crack number and hydraulic fracture branches after the fracturing under different in-situ stress states are listed in Table 3. In any injection sequence scheme, before the borehole pressure increased to the breakdown pressure, the fluid injected into the wells had no time to migrate in large area. Consequently, no evidence was found for the clear correlation between the breakdown pressure and injection sequence. As more fluid seeped into the rock matrix in the simultaneous injection, the effective stress between particles was markedly weakened, making the fluid pressure required for hydraulic fracture propagation lower than that in asynchronous fracturing. In the process of asynchronous fracturing, if the hydraulic fractures induced in the first injection have penetrated into the well for the second injection, the fluid pressure in the second injection well was hold at the fracture propagation pressure and the fractures extended tardily. Overall, the propagations of the hydraulic fractures between the two wells can be divided into four modes: (1) extending linearly and directly penetrating the two wells; (2) extending linearly at different heights without connection; (3) extending at certain deflection angles and merging at the tips; (4) suspending the extension after initiation. Mode (2) and (3) are suggested to improve the complexity of fracture networks. By comparing the hydraulic fracture propagation modes for the three schemes, both simultaneous fracturing and asynchronous fracturing under different in-situ stresses were observed to have the ability to form multiple fracture branches, which depended on the superposition of fluid pressures. Under the hydrostatic pressure condition, the initial expansion direction of hydraulic fractures in the asynchronous fracturing is affected by the micro-defects around the borehole, largely leading to no fractures between the two wells. In contrast, simultaneous fracturing reduces the effective stress of the rock matrix and forces the hydraulic fractures to extend to this region. When the maximum in-situ stress is in the horizontal direction, the initial expansion direction of hydraulic fractures is also along the horizontal direction. In this case, simultaneous fracturing intensifies the coalescence of the fractures between the two wells and simplifies the fracture networks. Therefore, selecting the reasonable injection sequence according to the in-situ stress condition is necessary for stimulating reservoir to produce more hydraulic fractures.

The Influence of Well Spacing under Different In-Situ Stress Conditions
The two-well fracturing models with the well spacing L of 0.3 m, 0.4 m and 0.6 m were established respectively, and the micro mechanical parameters, permeability, injection rate and in-situ stress conditions consistent with those in Section 4.2 were allocated to the models. Simultaneous fracturing was conducted on each model without considering the influence of injection sequence. Figures 15  and 16 show the fracturing simulation results of the granite models with diverse well spacing under different in-situ stress conditions. Refer to Figures 9 and 12 for the simultaneous fracturing results when the well spacing is 0.5 m.  Counting the hydraulic fracture extending to the right edge, three fracture branches were developed in the model. When L = 0.5 m, the hydraulic fractures between the two wells were also distorted, but did not meet before the simulation stopped as a consequence of the larger well spacing. When the well spacing grew to 0.6 m, only two hydraulic fractures extended to the left and right edges of the model, and the microcracks between the two wells stagnated.
Under the in-situ differential stress condition, each model after fracturing had the similar propagation pattern for fractures. Initially, four fractures extended horizontally from both sides of the two injection wells. Then the two hydraulic fractures between the wells merged into one fracture. Finally, three major fracture branches crossed the whole model. With the increase of well spacing, the micro defects around the wells control the crack initiation, thereby creating the hydraulic fractures at the different heights to merge after deflection.   Counting the hydraulic fracture extending to the right edge, three fracture branches were developed in the model. When L = 0.5 m, the hydraulic fractures between the two wells were also distorted, but did not meet before the simulation stopped as a consequence of the larger well spacing. When the well spacing grew to 0.6 m, only two hydraulic fractures extended to the left and right edges of the model, and the microcracks between the two wells stagnated.
Under the in-situ differential stress condition, each model after fracturing had the similar propagation pattern for fractures. Initially, four fractures extended horizontally from both sides of the two injection wells. Then the two hydraulic fractures between the wells merged into one fracture. Finally, three major fracture branches crossed the whole model. With the increase of well spacing, the micro defects around the wells control the crack initiation, thereby creating the hydraulic fractures at the different heights to merge after deflection.
The representative fluid pressures and crack information for the models with different well spacing after the fracturing are listed in Table 4. Before the continuous propagation of hydraulic fractures, there were too few microcracks to build the hydraulic connection between the two wells. Not surprisingly, the breakdown pressures under the same in-situ stress condition for these models remained steady as the well spacing changes. By contrast, a positive correlation was found between well spacing and fracture propagation pressure under the hydrostatic pressure condition. In the fracturing of the models with small well spacing, due to the superposition of the leaking fluid in the middle of the model, the effective stress of the rock matrix in this region and the fluid pressure required for hydraulic fracture propagation were correspondingly reduced. Obviously, this superimposed effect of fluid was suppressed when the well spacing was set larger. In addition, this phenomenon of fracture propagation pressure had not been monitored when the hydraulic fractures were controlled by the maximum horizontal in-situ stress. Taken together, these results suggest that there is an association between well spacing and fracture branches. Too small or too large well spacing makes it difficult to construct complex fracture networks in the reservoir. Relatively speaking, to ensure the hydraulic fractures to extend in the area between the two wells, the well spacing under the hydrostatic pressure may be smaller than that under the in-situ differential stress.

Analysis of Stress Shadow Effect
The fluid net pressure during the propagation of hydraulic fractures has contributed to the increase in the stress of the rock matrix around the fractures along the height direction. This important behavior of fractures is defined as "stress shadow effect", which has non-negligible impacts on the extension directions, opening degrees and shapes [66][67][68]. Especially in the situation of multiple injection wells in the reservoir, the interaction of the stress shadows for the adjacent hydraulic fractures becomes more intense and complicated.
However, discrete particle elements in the DEM suppress the direct expression of continuum physical quantities including stress and strain rate. The measurement circle is a monitoring mechanism built in PFC software to describe these quantities in a specified circular area by tracking the forces and displacements of particles and their related contacts. Therefore, the measurement circles densely distributed in the model realize the transition from the discrete physical quantities in the micro scale Energies 2020, 13, 4718 21 of 26 to the continuous ones in the macro scale. In order to obtain the rich data of measurement points, a total of 14,751 measurement circles with a radius of 0.01 m were regularly arranged in 99 rows and 149 columns to cover each two-well fracturing model (Figure 8b). The stress distributions of the models after fracturing were compared with their initial states without fluid injection to obtain the final stress increments.
Typical vertical stress increment distributions (∆σ y ) under different injection sequence and well spacing conditions are presented in Figures 17 and 18, respectively. All the models demonstrated a noticeable rise in the vertical stress on the upper and lower outer sides of the hydraulic fractures, which was in line with the basic law of stress shadow effect and proved the reliability of the modified fluid-mechanical coupling algorithm again. These vertical stress increments gradually decreased with the increasing distances to the fractures. In the area near the model boundaries, the vertical stress was roughly unchanged. Nevertheless, in the small area closest to the hydraulic fractures, the vertical stress was reduced rather than increased because of the leakage of fluid into the rock matrix. The injection of fluid led to an evident increase of vertical stress in the area around the wells without hydraulic fractures as well, which was higher than that caused by the hydraulic fractures.
Energies 2020, 13, x FOR PEER REVIEW 21 of 26 distributions of the models after fracturing were compared with their initial states without fluid injection to obtain the final stress increments. Typical vertical stress increment distributions (Δσy) under different injection sequence and well spacing conditions are presented in Figures 17 and 18, respectively. All the models demonstrated a noticeable rise in the vertical stress on the upper and lower outer sides of the hydraulic fractures, which was in line with the basic law of stress shadow effect and proved the reliability of the modified fluid-mechanical coupling algorithm again. These vertical stress increments gradually decreased with the increasing distances to the fractures. In the area near the model boundaries, the vertical stress was roughly unchanged. Nevertheless, in the small area closest to the hydraulic fractures, the vertical stress was reduced rather than increased because of the leakage of fluid into the rock matrix. The injection of fluid led to an evident increase of vertical stress in the area around the wells without hydraulic fractures as well, which was higher than that caused by the hydraulic fractures. More concretely, the injection sequence and well spacing are the dominant factors in the stress shadow effect. For the simultaneous fracturing when σv = 15 MPa and σh = 30 MPa (Figure 17a), no expected superposition of the stress shadow effects existed between the two adjacent hydraulic fractures. Under such superposition effect, the two parallel hydraulic fractures would repel each other and progressively separate [69]. A possible explanation for this is that the tensile stress fields at the fracture tips and the fluid leakage effect both acted on the area between the two fractures, which diminished the vertical stress and promoted the mutual attraction of fractures. Excluding this area, the superposition effect of stress shadows was widely scattered in the middle of the model. For the first injection in Well No.1 (Figure 17b), a large range of vertical compressive stress declined at the tip of fracture on the right side of Well No.2, resulting from the concentration of tensile stress and the climb of pore fluid pressure. The initial compressive stress state in the model have been altered into the tensile stress state near the fracture tip, whereas the decrease of vertical stress in the region far from the tip may not satisfy the alteration of stress state. For the second injection in Well No.2 ( Figure  17c), the stress shadow effect of the hydraulic fractures around Well No.1 was weakened stemming from the removal of preceding fluid and the primary vertical stress increment encompassed Well distributions of the models after fracturing were compared with their initial states without fluid injection to obtain the final stress increments. Typical vertical stress increment distributions (Δσy) under different injection sequence and well spacing conditions are presented in Figures 17 and 18, respectively. All the models demonstrated a noticeable rise in the vertical stress on the upper and lower outer sides of the hydraulic fractures, which was in line with the basic law of stress shadow effect and proved the reliability of the modified fluid-mechanical coupling algorithm again. These vertical stress increments gradually decreased with the increasing distances to the fractures. In the area near the model boundaries, the vertical stress was roughly unchanged. Nevertheless, in the small area closest to the hydraulic fractures, the vertical stress was reduced rather than increased because of the leakage of fluid into the rock matrix. The injection of fluid led to an evident increase of vertical stress in the area around the wells without hydraulic fractures as well, which was higher than that caused by the hydraulic fractures. More concretely, the injection sequence and well spacing are the dominant factors in the stress shadow effect. For the simultaneous fracturing when σv = 15 MPa and σh = 30 MPa (Figure 17a), no expected superposition of the stress shadow effects existed between the two adjacent hydraulic fractures. Under such superposition effect, the two parallel hydraulic fractures would repel each other and progressively separate [69]. A possible explanation for this is that the tensile stress fields at the fracture tips and the fluid leakage effect both acted on the area between the two fractures, which diminished the vertical stress and promoted the mutual attraction of fractures. Excluding this area, the superposition effect of stress shadows was widely scattered in the middle of the model. For the first injection in Well No.1 (Figure 17b), a large range of vertical compressive stress declined at the tip of fracture on the right side of Well No.2, resulting from the concentration of tensile stress and the climb of pore fluid pressure. The initial compressive stress state in the model have been altered into the tensile stress state near the fracture tip, whereas the decrease of vertical stress in the region far from the tip may not satisfy the alteration of stress state. For the second injection in Well No.2 ( Figure  17c), the stress shadow effect of the hydraulic fractures around Well No.1 was weakened stemming from the removal of preceding fluid and the primary vertical stress increment encompassed Well  (Figure 17a), no expected superposition of the stress shadow effects existed between the two adjacent hydraulic fractures. Under such superposition effect, the two parallel hydraulic fractures would repel each other and progressively separate [69]. A possible explanation for this is that the tensile stress fields at the fracture tips and the fluid leakage effect both acted on the area between the two fractures, which diminished the vertical stress and promoted the mutual attraction of fractures. Excluding this area, the superposition effect of stress shadows was widely scattered in the middle of the model. For the first injection in Well No.1 (Figure 17b), a large range of vertical compressive stress declined at the tip of fracture on the right side of Well No.2, resulting from the concentration of tensile stress and Energies 2020, 13, 4718 22 of 26 the climb of pore fluid pressure. The initial compressive stress state in the model have been altered into the tensile stress state near the fracture tip, whereas the decrease of vertical stress in the region far from the tip may not satisfy the alteration of stress state. For the second injection in Well No.2 (Figure 17c), the stress shadow effect of the hydraulic fractures around Well No.1 was weakened stemming from the removal of preceding fluid and the primary vertical stress increment encompassed Well No.2. As a consequence, the new hydraulic fracture did not deflect to the old one, but moved away. The comparison of the fracturing under different injection sequences reveals that the influence range of stress shadow effect in simultaneous fracturing exceeds that in asynchronous fracturing.
When σ v = 30 MPa and σ h = 30 MPa, under the condition of small well spacing (Figure 18a), the superposition of stress shadow effect between the two injection wells enhanced greatly and restricted the initiation of hydraulic fractures in other directions. As the well spacing was slightly increased (Figure 18b), the superposition range of stress shadow effect reduced and the high vertical stress on the upper and lower sides of the wells prohibited the fractures in the middle of the model from twisting outwards. When L = 0.6 m (Figure 18c), no interaction between the stress shadow effects was produced by different fractures. Therefore, the superposition degree of stress shadow effects between the two wells is in inverse proportion to their well spacing.

Discussion and Conclusions
The two-dimensional fluid-mechanical coupling algorithm in PFC has been extensively used to study the flow of fracturing fluid in rock, hydraulic fracture propagation and induced seismicity [33,[35][36][37][38][39][40][41][42][43][44][45][46][47][53][54][55]. In theory, it is practicable to conduct the two-dimensional simplified treatment for the relatively homogeneous and isotropic rocks in the macro scale, for the reason that the two-dimensional models can be regarded as a section of the real three-dimensional domains. Furthermore, two-dimensional simulation has the advantages of simple pre-processing, efficient calculation and intuitive results presentation in solving certain specific problems compared with three-dimensional simulation.
In this paper, the two-dimensional numerical models of Alxa porphyritic granites were strictly generated in PFC based on the PBM and the modified fluid-mechanical coupling algorithm. The micro input parameters in these models were effectively calibrated by the laboratory test results. Moreover, the accuracy of injection-induced fracturing was validated against the analytical solutions. A series of simulations for hydraulic fracturing in two-well granite models were performed to investigate the influences of injection sequence and well spacing on the breakdown pressure, fracture propagation pressure, fracture propagation pattern and stress shadow. The main results of this study are summarized as follows: (1) The injection sequence and well spacing have significant effects on the fracture propagation pressure instead of the breakdown pressure. Under the condition of hydrostatic pressure, lower fracture propagation pressure is required by the simultaneous fracturing models or small well spacing models. But these effects generally decay as the vertical in-situ stress declines. (2) Adjusting the injection sequence is a feasible way to control the propagation direction of hydraulic fractures. In the state of hydrostatic pressure, simultaneous fracturing allows more fracturing fluid to penetrate into the surrounding rock and compels the fractures to deflect to the horizontal direction, which is beneficial for the formation of complex fracture networks between the injection wells. The initial expansion of the fractures in asynchronous fracturing is easily disturbed by the micro defects around the boreholes, largely engendering no fractures between these wells. When the maximum in-situ stress is in the horizontal direction, in contrast to the asynchronous fracturing, simultaneous fracturing intensifies the coalescence of horizontal fractures. (3) The reasonable well spacing in the reservoir affects the growth of fracture branches. Excessively small or large well spacing limits the generation of multiple fracture branches in fracture networks. Comparatively speaking, for getting more fractures between the injection wells, the expected well spacing under the hydrostatic pressure may be smaller than that under the in-situ differential stress. (4) The stress shadow effects of hydraulic fractures are also closely related to the injection sequence and well spacing. In the case of simultaneous fracturing or small well spacing, the superposition of the stress shadow effects in the middle of the model strengthens their impact. When the fluid injection decreases or the well spacing increases, this superposition will be suppressed.
In summary, the findings of this study not only reveal the propagation and interaction mechanisms of hydraulic fractures in multi-well fracturing, but also provide valuable reference for the optimization of fracturing technology in field production. However, the two-dimensional models still have an inevitable defect, that is, the hydraulic fractures are forced to expand in the selected section, which fails to represent the three-dimensional expansion patterns. Further works associated with the two-dimensional and three-dimensional multi-well fracturing simulations, such as the interaction between hydraulic and natural fractures, the temperature effect of fracturing fluid, and the proppant migration, should be involved in future research.