*4.1. Modeling Scenarios*

The multi-well hydraulic fracturing behaviors of reservoir rocks not only a ffected 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 ge<sup>t</sup> 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 di fferent in-situ stress states on the breakdown pressure, fluid pressure distribution and hydraulic fracture propagation were discussed.

**Figure 8.** Granite model with two injection wells for hydraulic fracturing: (**a**) the model size and applied boundary conditions; and (**b**) arrangemen<sup>t</sup> of measurement circles.

#### *4.2. The Influence of Injection Sequence under Di*ff*erent 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.1 and Well No.2; Scheme 2—first injection in Well No.1 and second injection in Well No.2; 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−<sup>5</sup> m<sup>3</sup>/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 e ffect of residual pore pressure.

The influence of di fferent 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.

Figures 9–11 summarize the borehole pressure histories, crack spatial distribution and corresponding fluid pressure field for the fracturing of two-well granite model with the well spacing of 0.5 m under the hydrostatic pressure condition (<sup>σ</sup>v = σh = 30 MPa).

For the Scheme 1 of injection sequence, the borehole pressure curves recorded in Well No.1 and Well No.2 were nearly coincident in the process of simultaneous injection, and the breakdown pressures of the two wells were 62.0 MPa and 57.7 MPa (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.

**Figure 9.** The borehole pressure histories, crack spatial distribution and corresponding fluid pressure field in the simultaneous fracturing under the hydrostatic pressure condition (<sup>σ</sup>v = σh = 30 MPa).

**Figure 10.** The borehole pressure histories, crack spatial distribution and corresponding fluid pressure field in the asynchronous fracturing (first injection: Well No.1 and second injection: Well No.2) under the hydrostatic pressure condition (<sup>σ</sup>v = σh = 30 MPa).

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.

**Figure 11.** The borehole pressure histories, crack spatial distribution and corresponding fluid pressure field in the asynchronous fracturing (first injection: Well No.2 and second injection: Well No.1) under the hydrostatic pressure condition (<sup>σ</sup>v = σh = 30 MPa).

Figures 12–14 present the borehole pressure histories, crack spatial distribution and corresponding fluid pressure field under the conditions of the three injection schemes when σv = 15 MPa and σh = 30 MPa. The breakdown pressures of the two wells were slightly reduced compared with those under hydrostatic pressure.

**Figure 12.** The borehole pressure histories, crack spatial distribution and corresponding fluid pressure field in the simultaneous fracturing under the in-situ differential stress condition (<sup>σ</sup>v = 15 MPa and σh = 30 MPa).

**Figure 13.** The borehole pressure histories, crack spatial distribution and corresponding fluid pressure field in the asynchronous fracturing (first injection: Well No.1 and second injection: Well No.2) under the in-situ differential stress condition (<sup>σ</sup>v = 15 MPa and σh = 30 MPa).

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, ye<sup>t</sup> 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. Stopping the injection in Well No.1 and continuing to fracture the Well No.2, the other two hydraulic fractures were produced along the horizontal direction. Due to the interference of the existing hydraulic fractures, the initiation position of the new fracture in the right part of the model has been changed in contrast to the fracture simulated by the Scheme 1. At about 13 s of the simulation time, a fluid pressure of 1.8 MPa was observed in Well No.1, which implies that the new hydraulic fracture between the injection wells extended from Well No.2 has been connected with Well No.1. A large quantity of fracturing fluid flowed into the Well No.1 through this fracture, and further migrated to fill the hydraulic fractures generated by the first fracturing. The four fracture branches in the model aggravated the leakage of fluid to the surrounding formation and prevented the borehole pressure of Well No.1 from increasing.

**Figure 14.** The borehole pressure histories, crack spatial distribution and corresponding fluid pressure field in the asynchronous fracturing (first injection: Well No.2 and second injection: Well No.1) under the in-situ differential stress condition (<sup>σ</sup>v = 15 MPa and σh = 30 MPa).

The breakdown pressure of the Well No.2 was 43.3 MPa in the Scheme 3 (Figure 14). Unlike the fracturing in the Scheme 1 and Scheme 2, the hydraulic fracture between the two wells directly penetrated into the Well No.1 during the first injection, so that the borehole pressure of Well No.1 gradually rose to approximately 17 MPa. It is noteworthy that a new hydraulic fracture emerged from the outside of Well No.1 and propagated to the left edge of the model without reaching the previous breakdown pressures. This may be attributed to two aspects: on the one hand, the seepage of fracturing fluid has increased the pore fluid pressure and decreased the effective stress of the granite model; on the other hand, the rate of the fluid injection from Well No.2 into Well No.1 through the fracture was relatively low, which provided enough time for the stress of the solid framework around the Well No.1 to adjust. Additionally, the fluid pressure was found to decay steeply when flowing in the fractures on the left side of Well No.2. In the absence of sufficient fluid pressure, the propagation speed of the corresponding fractures was slow. During the second injection, the existing hydraulic fractures induced the fracturing fluid leak into the surrounding rock and limited the growth of borehole pressures for two wells. For example, the borehole pressure of Well No.1 was maintained near the fracture propagation pressure, which promoted the hydraulic fracture in the left part of the model to extend to the edge, but inhibited the generation of new hydraulic fracture.

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.


**Table 3.** Summary of the Fracturing Characteristics for the Two-Well Granite Model with the Well Spacing of 0.5 m under Different Injection Sequences.

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.

#### *4.3. The Influence of Well Spacing under Di*ff*erent 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.

**Figure 15.** The crack spatial distributions and corresponding fluid pressure fields for the two-well granite models with the well spacing of 0.3 m, 0.4 m and 0.6 m in the simultaneous fracturing under the hydrostatic pressure condition (<sup>σ</sup>v = σh = 30 MPa).

**Figure 16.** The crack spatial distributions and corresponding fluid pressure fields for the two-well granite models with the well spacing of 0.3 m, 0.4 m and 0.6 m in the simultaneous fracturing under the in-situ differential stress condition (<sup>σ</sup>v = 15 MPa and σh = 30 MPa).

Under the hydrostatic pressure condition, the hydraulic fractures in these models exhibit the interesting propagation patterns. When L = 0.3 m, plenty of microcracks initiated around the Well No.1 and Well No.2. Two hydraulic fractures extended along the horizontal direction and merged into one fracture in the middle of the model. The fracturing fluid injected into Well No.1 flowed towards Well No.2 through the hydraulic fracture and assisted the hydraulic fracture on the right side of Well No.2 to extend to the model boundary, while other microcracks near Well No.1 were not allowed to further extend owing to insufficient fluid pressure. When L = 0.4 m, the hydraulic fractures initiated from the micro defects on the surface of the two wells extended along the initial directions of different heights and deflected to each other. Until the fracture from Well No.1 reached the bottom of Well No.2 and the fracture from Well No.2 was arrested by the fracture from Well No.1, the two wells were more closely connected. 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 di fferential 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 di fferent heights to merge after deflection.

The representative fluid pressures and crack information for the models with di fferent 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 e ffective stress of the rock matrix in this region and the fluid pressure required for hydraulic fracture propagation were correspondingly reduced. Obviously, this superimposed e ffect 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 sugges<sup>t</sup> that there is an association between well spacing and fracture branches. Too small or too large well spacing makes it di fficult 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 di fferential stress.


**Table 4.** Summary of the Fracturing Characteristics for the Two-Well Granite Model with Di fferent Well Spacing under the Simultaneous Fracturing.

#### **5. Analysis of Stress Shadow E** ff**ect**

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 e ffect", which has non-negligible impacts on the extension directions, opening degrees and shapes [66–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

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 ( <sup>Δ</sup>σy) under di fferent 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 e ffect 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.

**Figure 17.** Vertical stress increment distributions for the two-well granite model with the well spacing of 0.5 m under di fferent injection sequences when σv = 15 MPa and σh = 30 MPa: (**a**) simultaneous fracturing; (**b**) first injection in Well No.1; and (**c**) second injection in Well No.2.

**Figure 18.** Vertical stress increment distributions for the two-well granite models with the different well spacing in the simultaneous fracturing when σv = 30 MPa and σh = 30 MPa: (**a**) L = 0.3 m; (**b**) L = 0.5 m; and (**c**) L = 0.6 m.

More concretely, the injection sequence and well spacing are the dominant factors in the stress shadow e ffect. For the simultaneous fracturing when σv = 15 MPa and σh = 30 MPa (Figure 17a), no expected superposition of the stress shadow e ffects existed between the two adjacent hydraulic fractures. Under such superposition e ffect, 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 e ffect 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 e ffect 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 e ffect 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 di fferent injection sequences reveals that the influence range of stress shadow e ffect 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 e ffect 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 e ffect 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 e ffects was produced by di fferent fractures. Therefore, the superposition degree of stress shadow e ffects between the two wells is in inverse proportion to their well spacing.

#### **6. 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–47,53–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, e fficient 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 e ffectively 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:


expected well spacing under the hydrostatic pressure may be smaller than that under the in-situ di fferential stress.

(4) The stress shadow e ffects 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 e ffects 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 e ffect of fracturing fluid, and the proppant migration, should be involved in future research.

**Author Contributions:** Conceptualization, L.Z. and J.Z.; methodology, S.W.; software, J.Z.; validation, S.W. and Z.H.; investigation, S.W.; resources, L.Z.; data curation, Z.H.; writing—original draft preparation, S.W.; writing—review and editing, J.Z.; visualization, S.W.; supervision, L.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China (Grant Nos. 41672321, 41972287, 41572312), China Postdoctoral Science Foundation (Grant Nos. 2018M630204, 2019T120133) and National Key Research and Development Program (Grant No.2018YFB1501801).

**Conflicts of Interest:** The authors declare no conflict of interest.
