Research on the Interaction Mechanism of Multi-Fracture Propagation in Hydraulic Fracturing

Abstract

During the hydraulic-fracturing process, stress interference occurs among multiple wells and fractures, potentially affecting the trajectory of hydraulic fracture propagation. Previous studies have largely overlooked the influence of proppant support stresses on the trajectories of fracture propagation. This paper establishes a mathematical model, grounded in the boundary element method, designed to compute the propagation of multiple fractures, considering both proppant support on the fracture surface and dynamic perturbations within the local stress field. The findings of this research reveal that the stress field induced by hydraulic fracturing exhibits dynamic evolution characteristics, necessitating a comprehensive study of the fracture initiation and extension across the entire fracturing time domain. The effect of the residual fracture width under proppant action on the in situ stress field cannot be ignored. During simultaneous fracturing, hydraulic fractures are inclined to propagate in the direction of the maximum horizontal principal stress, particularly as the in situ differential stress escalates. Staggered fracturing between wells has been proven to be more effective than head-to-head fracturing. Simply increasing the well spacing cannot solve the problem of inter-well fracture interference. In zipper fracturing, adjusting the fracturing sequence can inhibit the fracture intersections between wells, thereby controlling the trajectory of fracture propagation. The aforementioned research has considerable significance in guiding the control of fracture morphology during hydraulic-fracturing processes.

1. Introduction

China is not only a major producer but also a significant consumer of petroleum. Data indicate that the country’s external dependence on crude oil has steadily increased, rising from 60% in 2015 to 73.6% in 2020 [1]. Enhancing the efficiency of oil and gas development is not only a crucial strategic requirement for ensuring national energy security but also a fundamental necessity for sustaining the rapid development of the country’s economy and infrastructure. China possesses abundant unconventional oil and gas resources. However, these resources often exhibit low porosity and permeability, necessitating hydraulic fracturing to establish channels for the flow of oil and gas within the reservoir. In addition, acidification is also an important method to improve oil recovery performance [2], which often improves reservoir permeability through dissolution. During hydraulic fracturing, the interaction between fractures can influence their respective initiation and propagation directions. Currently, various interference effects of hydraulic fractures (HFs) exist in techniques such as simultaneous fracturing [3], segmented multi-cluster fracturing [4], zipper fracturing [5], and infill well fracturing [6], significantly impacting the productivity of oil and gas wells [7,8].

During segmented multi-cluster fracturing of a single well, hydraulic-fracturing fluid is simultaneously injected into each cluster fracture within each fracturing segment. Due to the mutual interference between cluster fractures, the direction of outward propagation of the peripheral fractures may change [9]. Additionally, due to the uneven fluid injection between clusters, there is competition for initiation between cluster fractures. Due to the heterogeneity of the formation conditions, the amount of fluid entering the fractures during the fracturing process is not exactly the same. The heterogeneity of the fluid inflow results in the different opening degrees and expansion capabilities of fractures between clusters. As fracturing construction progresses, cracks with strong expansion capabilities will quickly form macro-scale cracks, inhibiting the cracking of other cracks on both sides. The above is the competitive expansion mechanism between cracks. Fractures that do not have initiation advantage may stop propagating under induced stress or intersect and merge with other HFs [10]. Olson [11] pointed out that to generate vertical longitudinal fractures perpendicular to the wellbore, there must be sufficient spacing between the fractures during fracturing to avoid mutual interference during expansion. Conversely, if complex fracture morphologies are desired, the entry points for fractures should be arranged within the range of induced stress from other fractures. Simultaneous fracturing refers to the simultaneous perforation and fracturing of two parallel horizontal wells. Li et al. [12], through simulation, demonstrated that there is a tendency for intersecting fractures between wells during their mutual approach in simultaneous fracturing. Yang et al. [13] studied the fracturing efficiency during simultaneous and sequential fracturing of vertical wells, finding that the three-dimensional fracture area produced by simultaneous fracturing is much smaller than that produced by sequential fracturing. In zipper fracturing, fractures on two parallel horizontal wells are sequentially opened. Cheng et al. [14], based on the stiffness of the HF support, established a mathematical model for fracture expansion during sequential fracturing. Guo et al. [5], through numerical simulation, proposed that the higher the brittleness of the reservoir rocks, the longer the length of the fracture expansion during simultaneous and zipper fracturing. Under the same reservoir conditions, the length of the fracture expansion during zipper fracturing is greater than that during simultaneous fracturing. Saberhosseini et al. [15], using ABAQUS extended finite element models, studied the influence of the stitching positions on the fracture expansion in zipper fracturing, suggesting that intersections are prone to occur when fractures between wells approach each other. It is necessary to control the injection volume during fracturing to prevent deviations from the original direction or intersections when multiple HFs expand. During the expansion of multiple fractures, issues such as fracture initiation and expansion are caused by changes in the local stress field, interference and the connectivity of fractures between wells. These conditions can usually be observed from microseismic maps. When the microseismic points overlap, it is considered that there is inter-well fracture interference [6,16].

Usually, in the field-fracturing operation, there are other horizontal wells that have completed fracturing operation around the fractured well. At the same time, the fracture is not completely closed after hydraulic fracturing, and there is still a certain residual fracture width under the action of the proppant, which will continue to affect the stress field around the HF. Therefore, it is necessary to consider the influence of the residual crack width of the HF formed in the previous fracturing operation on the induced stress field and fracture propagation trajectory. However, in past numerical simulations, there has been a lack of research on the mechanism of fracture interaction in the coexistence of multi-process fracturing operations. Meanwhile, during the hydraulic-fracturing process involving multiple wells and fractures, the in situ stress field undergoes continuous alterations. Previous studies have overlooked the interference of the spatiotemporal evolution characteristics of induced stress on fracture trajectories. In this paper, a computational model for the propagation of multiple fractures during hydraulic fracturing is developed, and based on this model, the interaction mechanism between multi-wells and multi-fractures is studied. Compared with the previous calculation models based on the boundary element method, this model considers the proppant supporting force and is designed to simulate the propagation behavior of HFs after the local stress field changes due to large-scale hydraulic fracturing.

In conclusion, this paper aims to establish a model that considers the dynamic changes in the horizontal stress field for the interaction of multiple HFs and to investigate the interaction mechanisms between HFs. The rest of the paper is structured as follows. Section 2 is dedicated to detailing the mathematical calculation model. Section 3 conducts the numerical simulation analysis and provides the discussion, exploring the spatiotemporal evolution patterns of induced stress fields and the characteristics of fracture expansion trajectories. Section 4 presents the conclusions drawn from the study.

2. Numerical Model

2.1. Model Assumption

In practice, the propagation of HFs is influenced by the interplay of various physical fields, making it challenging to create a comprehensive simulation model that incorporates all these factors. In this paper, a calculation model of plane fracture propagation under reservoir conditions is established. To streamline the computational model and concentrate on the morphology of fractures and stress interference in horizontal wells, we introduce the following assumptions to construct our numerical model. (1) Isothermal conditions are assumed, the effect of temperature on rock deformation is not considered. (2) The model is considered homogeneous and isotropic, thus excluding the variability in rock mechanical properties and the abnormal pressure resulting from geological tectonic movements. (3) The target formation is treated as a linearly elastic medium, where the relationship between the load and the deformation adheres to Hooke’s law [17]. (4) The fracture is modeled according to the KGD (the fracture length is much greater than the fracture height) [18] model, which posits that the HF is a vertical fracture with a constant height, and the horizontal plane meets the plane strain condition for fractures that are short relative to their height [19]. (5) It is assumed that the study area is a homogeneous infinite reservoir. (6) The fracture propagation process is governed by the principles of linear elastic fracture mechanics [20], employing the linear theory of elasticity to analyze cracks, with characteristic parameters such as the stress intensity factor used to assess crack growth. (7) The fluid within the fracture is incompressible and belongs to the Poiseuille flow between two parallel plates [21].

2.2. Calculation Model of the Supporting Force from the Proppant on the Fracture Surface

In the hydraulic-fracturing process, when the fracturing is stopped, the fracture width will be relatively reduced under the original ground stress, but the fracture surface will not be completely closed due to the proppant support, and there will be a certain residual fracture width. It is assumed that the total stress of the fracture surface is compressive stress. The stress and deformation model of the fracture surface, considering the supporting stiffness of the proppant in the fracture, is shown in Figure 1.

The normal and tangential discontinuous displacements of the fracture element satisfy the following equation [22]:

$$\left\{\begin{array}{cc}{\sigma }_n=-K_n\Delta D_n;& \Delta D_n<{\left[D_n\right]}_0\\ {\sigma }_s=-K_s\Delta D_s;& \Delta D_s<{\left[D_s\right]}_0\end{array}\right.$$

(1)

In the equation, s − n represents the local coordinate system of the fracture plane; ${\sigma }_n,{\sigma }_s$ are the respective support stresses acting on the fracture plane, Pa; $K_n,K_s$ are the normal and tangential stiffness of the proppant-supported fracture (supported fracture), Pa/m; $\Delta D_n,\Delta D_s$ represent the change in the normal and tangential discontinuous deformation of the fracture plane, m; and ${\left[D_n\right]}_0,{\left[D_s\right]}_0$ are the historical maximum normal and tangential discontinuous displacement of the fracture plane, m.

In order to establish a nonlinear calculation model for proppant stiffness, the proppant stiffness $K_n$ [23] can be expressed in the following form:

$$K_n=K_{n0}/{\left(1-\frac{∆D_n}{{\left[D_n\right]}_0}\right)}^2$$

(2)

In the equation, $K_{n0}$ represents the support stiffness during the original linear deformation, Pa/m; and $K_n$ denotes the nonlinear support stiffness, Pa/m.

2.3. Deformation of Hydraulic Fracture

It is assumed that there are m HFs in the reservoir. During the hydraulic-fracturing process, after a finite number of extensions, any HF contains N_m discrete fracture elements. In this case, N is used to represent the total number of fracture elements in m fractures, $N=N_1+N_2+\dots +N_m$. At the initial time, the fracture surfaces are only affected by the in situ stress and the pressure inside the fracture. The deformation equation of the fracture faces is then expressed as [10]:

$$\begin{array}{c}\left\{\begin{array}{c}{\displaystyle \sum _{j=1}^N}\left(A_{ss}^{i,j}{\left[D_s^j\right]}_0+A_{sn}^{i,j}{\left[D_n^j\right]}_0\right)={\sigma }_s^i=-\frac{1}{2}\left({\sigma }_h-{\sigma }_H\right)\mathit{sin}2{\theta }_i;\\ {\displaystyle \sum _{j=1}^N}\left(A_{ns}^{i,j}{\left[D_s^j\right]}_0+A_{nn}^{i,j}{\left[D_n^j\right]}_0\right)={\sigma }_n^i=p_f-{\sigma }_h{\mathit{sin}}^2{\theta }_i-{\sigma }_H{\mathit{cos}}^2{\theta }_i;\end{array}\right.\\ \left(i=\mathrm{1,2}\dots ,N\right)\end{array}$$

(3)

In the equation, ${\left[D_s^j\right]}_0$ represents the historical maximum tangential discontinuous displacement of the j-th node element during hydraulic fracturing, m; and ${\left[D_n^j\right]}_0$ represents the historical maximum normal discontinuous displacement (fracture width) of the j-th node element during hydraulic fracturing, m.

Solving the above system of equations yields the discontinuous displacement quantities ${\left[D_s^j\right]}_0$ and ${\left[D_n^j\right]}_0$ for the fracture elements. Based on the discontinuous displacement of the fracture elements, the total stresses ${\sigma }_{xx}$, ${\sigma }_{xy}$ and ${\sigma }_{yy}$ acting on the fracture node elements are calculated considering the effects of the in situ stress, fluid inside the fracture, and the induced stress between multiple fractures. The forces on the fracture surface in the global coordinates are then transformed into a local coordinate system along the fracture surface and perpendicular to the fracture surface (s − n coordinate system). The tangential and normal forces on the fracture plane can be expressed in the local coordinate system as follows:

$$\left\{\begin{array}{c}{\sigma }_s^i=\left({\left[{\sigma }_{yy}\right]}^i-{\left[{\sigma }_{xx}\right]}^i\right)sin{\beta }_i·cos{\beta }_i+{\left[{\sigma }_{xy}\right]}^i\left({cos}^2{\beta }_i-{sin}^2{\beta }_i\right);\\ {\sigma }_n^i={\left[{\sigma }_{xx}\right]}^i{sin}^2{\beta }_i-2{\left[{\sigma }_{xy}\right]}^{i}sin{\beta }_i·cos{\beta }_i+{\left[{\sigma }_{yy}\right]}^i{cos}^2{\beta }_i;\end{array}\right.$$

(4)

Judge the forces acting on the fracture element nodes at this point. In this study, it is assumed that after hydraulic fracturing, the HFs are completely filled with proppant along the length of the fracture. Fracture elements from the k-th to the (k + m)-th experience a total stress in the form of compressive stress. Under the combined effects of the in situ stress and proppant, the fracture faces undergo a closure of $\Delta D_n$. Meanwhile, the total stress on the fracture faces of other fracture elements results in an outward tensile force. The overall deformation of all the discrete fracture elements can be comprehensively represented in the following form:

$$\begin{array}{c}\left\{\begin{array}{c}{\displaystyle \sum _{j=k}^{k+m}}\left[A_{ss}^{i,j}\left({\left[D_s^j\right]}_0+\Delta D_s^j\right)+A_{sn}^{i,j}\left({\left[D_n^j\right]}_0+\Delta D_n^j\right)\right]+{\displaystyle \sum _{j=1}^{k-1}}\left(A_{ss}^{i,j}{D}_s^j+A_{sn}^{i,j}{D}_n^j\right)+{\displaystyle \sum _{j=k+m+1}^N}\left(A_{ss}^{i,j}{D}_s^j+A_{sn}^{i,j}{D}_n^j\right)={\sigma }_s^i;\\ {\displaystyle \sum _{j=k}^{k+m}}\left[A_{ns}^{i,j}\left({\left[D_s^j\right]}_0+\Delta D_s^j\right)+A_{nn}^{i,j}\left({\left[D_n^j\right]}_0+\Delta D_n^j\right)\right]+{\displaystyle \sum _{j=1}^{k-1}}\left(A_{ns}^{i,j}{D}_s^j+A_{nn}^{i,j}{D}_n^j\right)+{\displaystyle \sum _{j=k+m+1}^N}\left(A_{ns}^{i,j}{D}_s^j+A_{nn}^{i,j}{D}_n^j\right)={\sigma }_n^i;\end{array}\right.\\ \end{array}$$

(5)

In the equation, $\Delta D_s^j$ represents the change in the tangential discontinuous displacement on the j-th node element, m; and $\Delta D_n^j$ represents the change in the normal discontinuous displacement on the j-th node element, m.

For the fracture element i undergoing compression during hydraulic fracturing, the total stress it encounters is a result of the interplay among the horizontal principal stress, the pressure within the fracture, and the counteracting force exerted by the proppant:

$$\left\{\begin{array}{c}{\sigma }_s^i=-\frac{1}{2}\left({\sigma }_h-{\sigma }_H\right)\mathit{sin}2{\theta }_i-{\left({\sigma }_s^i\right)}_{stiff};\\ {\sigma }_n^i=p_i-{\sigma }_h{\mathit{sin}}^2{\theta }_i-{\sigma }_H{\mathit{cos}}^2{\theta }_i-{\left({\sigma }_n^i\right)}_{stiff};\end{array}\right. \left(i=k, k+1, \dots , k+m\right)$$

(6)

In the equation, ${\left({\sigma }_s^i\right)}_{stiff}$ represents the force exerted by the proppant normal to the fracture surface, Pa; and ${\left({\sigma }_n^i\right)}_{stiff}$ represents the force exerted by the proppant tangential to the fracture surface, Pa.

For fractures that do not shrink during hydraulic fracturing, they experience the influence of the horizontal principal stress and the pressure from the fluid inside the fracture:

$$\left\{\begin{array}{c}{\sigma }_s^i=-\frac{1}{2}\left({\sigma }_h-{\sigma }_H\right)\mathit{sin}2{\theta }_i;\\ {\sigma }_n^i=p_i-{\sigma }_h{\mathit{sin}}^2{\theta }_i-{\sigma }_H{\mathit{cos}}^2{\theta }_i;\end{array}\right. \left(i=1, 2,\dots ,k-1, k+m+1,\dots ,N\right)$$

(7)

Introduce the total discontinuous displacements of the fracture elements, $D_s^j$, $D_n^j$. The total fracture discontinuous displacement is equal to the sum of the maximum fracture width and fracture surface shrinkage during fracture propagation:

$$\left\{\begin{array}{cc}{D}_s^j={\left[D_s^j\right]}_0+\Delta D_s^j, D_n^j={\left[D_n^j\right]}_0+\Delta D_n^j;\hfill & \hfill j=k, k+1, \dots , k+m\\ D_s^j=D_s^j, D_n^j=D_n^j;\hfill & \hfill j=1, 2,\dots ,k-1, k+m+1,\dots ,N\end{array}\right.$$

(8)

Taking into account the discontinuous displacement of the fractures, the stress on the fracture faces is iteratively computed until the discontinuous displacement reaches equilibrium under the influence of the proppant, i.e., $\Delta D^{j+1}-\Delta D^j<\epsilon$ (allowable error). It is assumed that the compressed fracture element nodes are k′ to k′ + m′, and updating yields the final discontinuous displacement of the fracture elements ${\left[D_s^j\right]}^{\prime }$, ${\left[D_s^j\right]}^{\prime }$. According to the amount of discontinuous displacement, the fracture criterion in Section 2.4 is used to determine whether the crack tip is ruptured and extended, where:

$$\left\{\begin{array}{cc}{\left[D_s^j\right]}^{\prime }={\left[D_s^j\right]}_0+\Delta D_s^j, {\left[D_n^j\right]}^{\prime }={\left[D_n^j\right]}_0+\Delta D_n^j ;\hfill & \hfill j=k\prime , k\prime +1, \dots , k\prime +m\prime \\ {\left[D_s^j\right]}^{\prime }=D_s^j, {\left[D_n^j\right]}^{\prime }=D_n^j;\hfill & \hfill j=1, 2,\dots ,k^{\prime }-1, k\prime +m^{\prime }+1,\dots ,N\end{array}\right.$$

(9)

For the problem of multiple fracture expansion, when there are K fracture tips that are ruptured and extended, corresponding to each of them, K fracture units are added:

$$\left\{\begin{array}{c}{\displaystyle \sum _{j=k^{\prime }}^{k^{\prime }+m^{\prime }}}\left[A_{ss}^{i,j}\left({\left[D_s^j\right]}_0+\Delta D_s^j\right)+A_{sn}^{i,j}\left({\left[D_n^j\right]}_0+\Delta D_n^j\right)\right]+{\displaystyle \sum _{j=1}^{k^{\prime }-1}}\left(A_{ss}^{i,j}{D}_s^j+A_{sn}^{i,j}{D}_n^j\right)+{\displaystyle \sum _{j=k^{\prime }+m^{\prime }+1}^{N+K}}\left(A_{ss}^{i,j}{D}_s^j+A_{sn}^{i,j}{D}_n^j\right)={\sigma }_s^i;\\ {\displaystyle \sum _{j=k^{\prime }}^{k^{\prime }+m^{\prime }}}\left[A_{ns}^{i,j}\left({\left[D_s^j\right]}_0+\Delta D_s^j\right)+A_{nn}^{i,j}\left({\left[D_n^j\right]}_0+\Delta D_n^j\right)\right]+{\displaystyle \sum _{j=1}^{k^{\prime }-1}}\left(A_{ns}^{i,j}{D}_s^j+A_{nn}^{i,j}{D}_n^j\right)+{\displaystyle \sum _{j=k^{\prime }+m^{\prime }+1}^{N+K}}\left(A_{ns}^{i,j}{D}_s^j+A_{nn}^{i,j}{D}_n^j\right)={\sigma }_n^i ;\end{array}\right.$$

(10)

At this time, for the pressure fracture unit:

$$\left\{\begin{array}{c}{\sigma }_s^i=-\frac{1}{2}\left({\sigma }_h-{\sigma }_H\right)\mathit{sin}2{\theta }_i-K_s\Delta D_s^i;\\ {\sigma }_n^i=p_i-{\sigma }_h{\mathit{sin}}^2{\theta }_i-{\sigma }_H{\mathit{cos}}^2{\theta }_i-K_n\Delta D_n^i;\end{array}\right. \left(i=k\prime , k\prime +1, \dots , k\prime +m\prime \right)$$

(11)

For unpressurized fractures:

$$\left\{\begin{array}{c}{\sigma }_s^i=-\frac{1}{2}\left({\sigma }_h-{\sigma }_H\right)\mathit{sin}2{\theta }_i;\\ {\sigma }_n^i=p_i-{\sigma }_h{\mathit{sin}}^2{\theta }_i-{\sigma }_H{\mathit{cos}}^2{\theta }_i;\end{array}\right. \left(i=\begin{array}{c}1, 2,\dots ,k^{\prime }-1,\\ k^{\prime }+m^{\prime }+1,\dots ,N+K\end{array}\right)$$

(12)

Iterate the calculation until there is no fracture propagation at the HF tip or until the predetermined time step is reached.

2.4. The Initiation and Propagation of Fracture

The stress intensity factor (SIF) serves as a crucial parameter for characterizing the stress field at the fracture tip and is determined by the value of the discontinuous displacement on the fracture surface. The SIF at point a on the fracture tip can be illustrated as follows [24]:

$$\left\{\begin{array}{c}{K}_{\mathsf{{\rm I}}}=\frac{0.806E\sqrt{\pi }}{4\left(1-{\upsilon }^2\right)\sqrt{2a}}{D}_n\\ K_{\mathsf{{\rm I}}\mathsf{{\rm I}}}=\frac{0.806E\sqrt{\pi }}{4\left(1-{\upsilon }^2\right)\sqrt{2a}}{D}_s\end{array}\right.$$

(13)

In the equation, $K_{\mathsf{{\rm I}}}$ and $K_{\mathsf{{\rm I}}\mathsf{{\rm I}}}$ are the stress intensity factor of the fracture tip, MPa∙m^0.5; E is Young’s modulus, MPa; ν is Poisson’s ratio; and a is the perforating depth, m.

Simultaneously, the maximum circumferential stress (MCS) criterion is employed to determine the initiation and propagation direction of the fracture [25]. The fracture deflection angle (θ₀) can be calculated by the following formula [10]:

$${\theta }_0=arcsin\left[\frac{K_{\mathsf{{\rm I}}\mathsf{{\rm I}}}\left(K_{\mathsf{{\rm I}}}+3\sqrt{{K_{\mathsf{{\rm I}}}}^2+8{K_{\mathsf{{\rm I}}\mathsf{{\rm I}}}}^2}\right)}{{K_{\mathsf{{\rm I}}}}^2+9{K_{\mathsf{{\rm I}}\mathsf{{\rm I}}}}^2}\right]$$

(14)

The instability of the HF can be judged by the mixed stress intensity factor:

$$K_e=\cos \frac{{\theta }_0}{2}\left[K_{\mathsf{{\rm I}}}{\mathrm{c}\mathrm{o}\mathrm{s}}^2\frac{\theta }{2}-\frac{3}{2}{K}_{\mathsf{{\rm I}}\mathsf{{\rm I}}}\sin {\theta }_0\right]$$

(15)

where $K_e$ is the fracture toughness, MPa∙m^1/2.

When ${K_{IC}>K}_e$, the fracture tip is unstable.$K_{IC}$ is the fracture toughness, MPa∙m^1/2.

2.5. Numerical Calculation Flow

Figure 2 presents the flow chart detailing the solution procedure. The steps involved are as follows. ① Input the initial parameters, which include the maximum and minimum horizontal stresses, Young’s modulus, Poisson’s ratio, fracture toughness, the number of fractures, fracture spacing, injection pressure, perforation depth, and the angle of perforation. ② Hydraulic fracture discretization and establish a calculation model of the supporting force from the proppant (Equations (1) and (2)). ③ Establish an initial deformation control equation for the discrete crack element (Equation (3)). ④ Calculate the $D_n$, $D_s$ of the HF and induced stress ${\sigma }_{xx}$, ${\sigma }_{xy}$, ${\sigma }_{yy}$. ⑤ Establish a deformation control equation for multiple cracks, considering the supporting force from the proppant (Equations (5)–(8)). ⑥ Determine whether the fracture tip is unstable and ruptured (Equation (5)); if the answer is yes, calculate the crack deflection angle and increase the number of fracture units, and then go to step ⑤; otherwise, end the calculation.

2.6. Model Validation

At the outset, it is crucial to confirm the stress caused by fracture deformation. Our model assumes plane strain conditions and maintains a constant fracture height [20]. The analytical solution for the stress field surrounding the HF can be expressed as follows [26]:

$$\left\{\begin{array}{c}{\sigma }_x=-P_f\frac{r}{a}{\left(\frac{a^2}{r_1{r}_2}\right)}^{\frac{3}{2}}sin\theta sin\left[\frac{3}{2}\left({\theta }_1+{\theta }_2\right)\right]-P_f\left[\frac{r}{\sqrt{r_1{r}_2}}cos\left(\theta -\frac{{\theta }_1+{\theta }_2}{2}\right)-1\right]\\ {\sigma }_y=P_f\frac{r}{a}{\left(\frac{a^2}{r_1{r}_2}\right)}^{\frac{3}{2}}sin\theta sin\left[\frac{3}{2}\left({\theta }_1+{\theta }_2\right)\right]-P_f\left[\frac{r}{\sqrt{r_1{r}_2}}cos\left(\theta -\frac{{\theta }_1+{\theta }_2}{2}\right)-1\right]\\ {\sigma }_{xy}=-P_f\frac{r}{a}{\left(\frac{a^2}{r_1{r}_2}\right)}^{\frac{3}{2}}sin\theta cos\left[\frac{3}{2}\left({\theta }_1+{\theta }_2\right)\right]\end{array}\right.$$

(16)

where P_f denotes the pressure inside the HF, Pa; r, r₁, r₂, respectively, represent the distance from the fracture to point Q, m; and θ, θ₁, θ₂ represent the angle between the distance line and X-axis. The mechanical model is depicted in Figure 3a.

Initial simulation parameters: the perforating depth (a) is 1 m, the pressure inside the HF (P_f) is −3 MPa, the coordinates of point Q: x = 0.5 m, and y increases from 0 m to 10 m. Calculate the induced stress at point Q for varying vertical distances (coordinate y) from the fracture plane. The relevant description can be found in our previous research results [27]. In Figure 3b, the solid line represents analytical results and the dotted line represents the simulation results. The error between the numerical and analytical solutions is less than 5%.

Secondly, verify the calculation model for the supporting force from the proppant. Take $D_{n0}=1.3\times {10}^{-4}\mathrm{m}$, $K_{n0}=50 \mathrm{G}\mathrm{P}\mathrm{a}·{\mathrm{m}}^{-1}$. The stress–strain behavior of the fracture surface under compression is illustrated in Figure 4a. During the initial compression phase of the fracture surface, due to the loose contact between the proppant particles, the fracture surface tends to contract easily under relatively low pressure. However, when the contraction of the fracture surface reaches 1.2 × 10⁻⁴ m, the proppant is compacted and deformed. At this point, the stiffness of the fracture surface support tends to infinity with the increase in the contraction of the fracture width. The relationship between the variables is compared with the experimental results obtained on granite with diamond particles as proppants (Figure 4b), and the curve patterns show a good match.

3. Results and Discussion

The interaction between fractures during hydraulic fracturing is reflected in fracturing technologies such as single-well staged fracturing, multi-well simultaneous fracturing, and zipper fracturing. This section will carry out numerical simulation research on the above three fracturing technologies, analyzing the trajectories of fracture expansion and the variation law of the induced stress field. The initial computational models are illustrated in Figure 5. Figure 5a represents the initial fracture geometry model during segmented fracturing. Figure 5b represents the initial fracture geometry model during synchronous fracturing and zipper fracturing.

3.1. Space-Time Evolution Mechanism of Induced Stress Field between Fractures

The induced stress field during the hydraulic-fracturing process is not static; it evolves continuously with the expansion of the fractures. Mastering the evolutionary mechanism of the stress field is beneficial to fundamentally revealing the evolution mechanism of fractures during hydraulic fracturing. In this section, a numerical model is established by taking the fracture propagation of segmented fracturing in a single-level well as an example to study the spatiotemporal evolution mechanism of the induced stress field during hydraulic fracturing, as illustrated in Figure 5a. The hydraulic-fracturing operations are conducted in a sequential manner: first fracturing HF1 until it reaches the target length, then isolating it with a packer, before proceeding to fracturing HF2. This sequence is repeated to complete the pressure operations for three HFs. The initial simulation parameters utilized in this section are detailed in

Table 1. Initial simulation parameters.

Parameters	Value	Parameters	Value
Young’s modulus (MPa)	30,000	Minimum horizontal stress (MPa)	40
Poisson’s ratio	0.25	Fracture toughness (MPa·m1/2)	2.5
Maximum horizontal stress (MPa)	40	Injection pressure (MPa)	44
Perforating depth (m)	1	Fracture spacing (m)	10

. The simulation results are presented in Figure 6, which respect the induced stress field formed by the fracture surface expansion and shear deformation during fracturing. Figure 6a shows the simulation results at the completion of HF1, and Figure 6b shows the results when all three fractures have completed hydraulic fracturing. The letters A, B, C, D, and E in Figure 6 represent five observation points within the domain of the HF expansion, enabling the analysis of the variation patterns of the induced stress fields at the observation points during the fracturing process.

The induced stress values at each observation point during the dynamic expansion of the HFs were extracted, and the curves are plotted in Figure 7, where m_i represents the completion of the previous i stages of fracturing (i = 1, 2, 3); n_i represents the start of stage i + 1 fracturing, when the fracturing and plugging operation in the previous stage of the horizontal well has been completed; A, B, C, D, E correspond to five observation points within the influence range of induced stress by HF; and Δσ_xx represents the induced stress drop caused by the well shut-in process. t is the dimensionless time, where it represents the iterative step of the computation.

In Figure 7a, point A is situated between HF1 and HF2. Following the propagation of HF1, it generate induced tensile stress at point A. When HF1 propagates to a certain length, the induced stress at point A transitions between tensile and compressive stress. The compressive stress mainly comes from HF1, causing the induced stress at point A to continuously rise to 1.9 MPa until the completion of the expansion of fracture 1 (dimensionless time step t = 25). Subsequently, after shutting in the well and initiating the second stage of fracturing (t = 26), following the conclusion of the first stage of fracturing, and given the combined effect of the in situ stress and proppant support within the fracture, the fracture surface of HF1 shrinks, leading to a corresponding decrease in the induced stress at point A by Δσ_xx = 1 MPa. As HF2 propagates, the induced compressive stress at point A begins to increase continuously, ultimately reaching 3.1 MPa (t = 51). After the second stage of fracturing stops, the induced compressive stress at point A decreases (t = 52). During the propagation of HF3, the induced stress at point A is relatively small due to the greater distance. Point B is located between HF2 and HF3, and because points A and B are symmetrically distributed with respect to HF1, 2 and 3, the induced stresses at points A and B are roughly equal at 1.8 MPa after the completion of all three fracturing stages. However, point A experiences larger fluctuations in the induced compressive stress relative to point B throughout the fracturing stages. When there is significant fluctuation in the in situ stress, it is more likely to disrupt the reservoir stress balance, potentially activating natural fractures. Therefore, when there is a natural fracture group in the reservoir, it is easier to awaken the natural fracture group by fracturing the area near the natural fracture group first than by fracturing the peripheral area of natural fracture group first and then fracturing the center area of natural fracture group. During the fracturing process, point C experiences a maximum value of induced tensile stress (point Q). This is because point C is located on the extension line of HF2, and as the fracture tip approaches point C, the induced tensile stress rapidly increases. After the fracture tip passes through point C, the induced stress quickly transitions from tensile to compressive.

In Figure 7b, the coordinates of observation points A, D, and E vary along the vertical axis. Due to the consistent location of points D and E in the forward direction of the inclined fracture, they consistently experience induced stress dominated by tensile stress. During the well shut-in process, the induced stress will also experience a decrease. It should be noted that the induced stress generated at point D increases by Δσ′_xx = 0.75 MPa when HF2 is shut in after fracturing and begins the next stage of fracturing. This is because during the propagation of HF2, the residual width of fracture 1 is compressed, reducing the impact of the residual width on the induced stress at point D. When HF2 stops fracturing and the well is shut in, the compressive effect on HF1 is reduced, resulting in a sudden increase in the induced tensile stress at point D. It is evident that due to the fracture interaction, when the spacing between adjacent fractures is small, the deformation of the fractures reduces the induced tensile stress in the regions near their respective fracture tips.

In this section, the induced stress generated at fixed points in the study area during the staged fracturing process is analyzed. The findings indicate that, as the fractures propagate, the dynamic fluctuations of the induced stress field should not be overlooked.

When weak structural planes, including joints and micro-fractures, are present within the zone affected by the induced stress, the weak structural plane will cause damage and micro-fracturing with the propagate of HFs. Furthermore, with the dynamic fluctuations in the local induced stress field, the fractures exhibit different initiation capabilities at various fracturing moments.

3.2. Mechanism of Fracture Interaction in Synchronous Fracturing

In the well factory mode, it is not uncommon for multiple wells to experience simultaneous fracturing, necessitating the analysis of the fracture expansion trajectories during synchronous-fracturing processes. The physical model is depicted in Figure 5b. It is assumed that in the multi-cluster-fracturing process, each cluster forms an HF. The minimum and maximum horizontal principal stress are 40 MPa, respectively, the pressure inside the artificial fracture is 42 MPa, the fracture cluster spacing is 40 m, and the well spacing is 100 m. Other parameters are provided in Table 1.

Firstly, the interaction mechanism under the influence of the relative positions of the fractures between wells is studied. The lateral spacing between adjacent well fractures is varied at 0 m, 10 m, and 20 m. Figure 8 shows the simulation results, where the stress field in the figure corresponds to the induced stress generated by HFs, with the units in Pa. As can be observed in the figure, the fracture tips tend to intersect at the top when the fractures are in close proximity to one another. It can be seen from the figure that the cracks are easy to intersect in the arrangement mode shown in Figure 8a. In actual fracture perforation, the fractures between wells may not be perfectly aligned. The relative positions between the fractures are closer to those depicted in Figure 8b,c. In the initial stage, the propagation of fracture is mainly influenced by inter-cluster interference. After the fracture extension length reaches 25 m, the trajectory is affected by neighboring well fractures, ultimately causing the tips of the fractures to approach each other.

From this, it can be seen that whether the fractures between wells are arranged in a head-to-head manner or stagger, the fracture tips tend to intersect when they approach each other. When the fracture spacing is not large enough, changing the pattern of the arrangement cannot fundamentally solve the issue of intersecting fractures between wells.

Secondly, an analysis is conducted to examine the impact of the horizontal principal stress difference on the propagation trajectory of the fractures. The pressure in the HF is 42 MPa and the maximum horizontal principal stress is 40 MPa. The effect of the horizontal principal stress difference on the fracture propagation trajectory is simulated by changing the minimum horizontal principal stress. The simulation outcomes are depicted in Figure 9, where the stress field distribution (σ_xy) represents the induced stress field resulting from the HFs, with the units expressed in Pa.

Figure 9 demonstrates that variations in the horizontal principal stress significantly influence the fracture deviation during concurrent fracturing. The greater the stress differential, the lower the probability of fracture intersections occurring during synchronous fracturing. In the actual synchronous-fracturing process, in order to increase the volume of the transformation as much as possible, it is often necessary to stagger the adjacent well fracture so that they do not intersect easily when the fracture tips are approaching. However, for reservoirs with low horizontal principal stress differences, because the fractures between wells are easy to intersect when they are close to each other, it is necessary to relatively increase the fracture spacing.

3.3. Study on Fracture Propagation Trajectory of Zipper Fracturing

Compared to simultaneous fracturing, when it comes to zipper fracturing, the expansion of HFs does not occur simultaneously but rather one by one in a sequential manner. The characteristics of the induced stress variation during the sequential-fracturing process have been analyzed in Section 3.1. This section primarily focuses on simulating and studying the morphology of fracture expansion in the zipper-fracturing process. Figure 5b is the physical model, where two horizontal wells are oriented perpendicular to the direction of the maximum horizontal principal stress. The numbers 1 to 5 of HF1 to HF5 in Figure 5b indicate the fracturing sequence, with HF1 first being fractured, followed by HF2, and so on. In this example, HF1 is first opened, followed by the fracturing of HF2, and so on, until a total of seven fractures are opened. The rock mechanics parameters are provided in Table 1, the maximum horizontal principal stress is 40 MPa, the minimum horizontal principal stress is 36 MPa, and the internal pressure is 40 MPa. The well spacing (d4) is 80 m, and the fracture spacings (d5 and d6) are 20 m and 40 m, respectively.

Figure 10 shows the simulation results, where the stress field distribution (σ_yy) represents the induced stress field produced by the HFs, with the units in Pa.

Figure 10a illustrates the induced stress distribution when HF1 is about to stop fracturing. It is evident that during the process of zipper fracturing, the first HF propagates linearly, following the orientation of the maximum principal stress, which is dictated by the horizontal principal stress. After the completion of fracturing, the local stress field distribution changes. Subsequently, the propagation of HF2 is simulated on this basis, as shown in Figure 10b, where HF2 deviates toward HF1 under the influence of the local stress field. The simulation results for fractures 3, 4, 5, 6, and 7 are presented in Figure 10c,d, showing that subsequent fractures tend to close to each other. This phenomenon of proximity to each other during expansion can easily cause the intersection of HFs, resulting in the formation of inter-well interference during the fracturing process.

3.4. Sensitivity Analysis of Fracture Propagation Trajectory

To examine the impact of the well spacing and fracture spacing on the trajectory of fracture propagation during zipper fracturing, a fracturing simulation model was established with varying well spacing. To assess the effect of the well spacing on the trajectory of fracture propagation during the execution of multi-well zipper fracturing, based on the model established in Section 3.3 and the simulation parameters shown in Figure 10, numerical simulation was conducted by changing the well spacing. The well spacing is 60 m, 80 m and 120 m, respectively.

The simulation outcomes are depicted in Figure 11. The stress field is the induced stress field σ_yy (unit Pa) generated by the HFs. From the graph, it is evident that the fractures still exhibit a tendency to close to each other under different well spacings, and even doubling the well spacing does not effectively inhibit fracture deflection.

To address the issue of fracture intersection arising from the interaction between fractures during the process of zipper fracturing, the fracturing sequence is adjusted to simulate the fracture propagation trajectory. Figure 12a illustrates the fracture morphology observed during conventional zipper fracturing.

The hydraulic-fracturing sequence was adjusted to fracture 1-3-2-5-4-7-6 successively. The simulation results are presented in Figure 12b, where the stress distribution represents the induced stress σ_yy generated by the HFs. It can be observed from the graph that adjusting the fracturing sequence can make the fracture less prone to deflection when approaching adjacent well fractures. However, subsequent HF4, 5, 6, and 7 still did not maintain straight-line expansion, mainly due to the deflection of HF5. Therefore, the number of HFs should not exceed 5 in the control of the fracture trajectory of zipper fracturing.

3.5. Discussion

Hydraulic fracturing is an important method for unconventional oil and gas extraction. The macroscopic permeability of the reservoir can be improved by forming permeability channels with high conductivity within the reservoir. The oil and gas flow channels formed by fracturing are usually called hydraulic fractures. The trajectory and morphology of fractures are key factors affecting the reservoir reconstruction effect. Previous studies attributed the initiation and propagation of hydraulic fractures to local stress fields [29]. However, during the fracturing process, the in situ stress field is not constant [30]. As hydraulic fractures propagate, an induced stress field is formed inside the reservoir, changing the distribution of the original stress field, resulting in a more complex fracture expansion trajectory. In response to these problems, this paper established a multi-well fracture propagation model that considers the spatiotemporal evolution of the stress field. First, the spatiotemporal evolution law of the induced stress field during the fracturing process was analyzed through numerical simulation. Based on this evolution pattern, the interaction between fractures during staged fracturing, simultaneous fracturing and zipper fracturing was studied. The effects of stress fields and suturing methods on the crack propagation trajectories were studied. Compared with previous studies, this study is innovative in two aspects:

• A stress field calculation model, considering the supporting effect of the proppant, was established. After fracturing, the hydraulic fracture is not completely closed with the support of the proppant but will maintain a certain residual width, affecting the in situ stress field distribution. Previous multi-well fracturing numerical models did not consider the impact of these residual widths on the in situ stress field, resulting in deviations between the simulation results and the actual situation [31,32]. Especially in the case of repeated fracturing or infill well fracturing, the original stress field changes greatly, and previous hydraulic-fracturing calculation models cannot take into account the complex and uneven distribution characteristics of the stress field [33]. This model takes into account the influence of the residual width under proppant support on the stress field and overcomes the shortcomings of previous numerical simulations. It can analyze the stress field distribution in different periods by simulating the historical fracturing process of all the wells on site.
• The spatiotemporal evolution patterns of the induced stress field during the propagation of multiple fractures were analyzed, and based on this, the interaction between the fractures was studied. Previous studies on the stress field during hydraulic-fracturing processes have been limited to specific time points and have not considered the dynamic evolution of the stress field over time, resulting in an incomplete understanding of the fracture interactions during the fracturing process. Especially in reservoirs with multiple wells, multiple fractures, and different construction techniques, it becomes crucial to analyze the fracture interactions from the perspective of the spatiotemporal evolution of the induced stress field.

4. Conclusions

In this study, a mathematical model has been developed to simulate the propagation of fractures within a hydraulic-fracturing process where multiple wells and fractures are present. The space–time evolution law of the stress field and the fracture propagation trajectory characteristics are studied.

• Based on the boundary element method, a computational model for the propagation of multiple fracture was established using MATLAB programming. The model takes into account the supporting effect of the support agents on crack surfaces, the interaction between multiple cracks, and the dynamic perturbation of the local stress field.
• The induced stress field during the hydraulic-fracturing process exhibits dynamic fluctuation characteristics. In the analysis of the stress field and the interaction between fractures, it is necessary to study the initiation and propagation of fractures by considering the temporal and spatial variation characteristics of the stress field.
• In the process of synchronous fracturing, no matter whether the fractures are arranged head-to-head or staggered between adjacent wells, the fracture tips tend to intersect when they approach each other. The deflection angle of the fracture tip decreases with the increase in the horizontal stress difference.
• It is easy for the zipper fracturing to cause the inter-well fractures to intersect in the process of approaching each other, and the problem of inter-well fracture interference cannot be solved simply by increasing the well spacing. The fracture propagation trajectory can be controlled to some extent by adjusting the fracturing sequence in zipper fracturing.