Numerical Simulation on Multi-Fractures Propagation Behavior Based on Hybrid Finite-Discrete Method in Deep Shale Formation

1. Introduction

Deep shale resources played an important role for the energy supply [1]. The primary technology for developing deep shale reservoirs was hydraulic fracturing, which involved high-pressure water injection into low-permeability rock formations at depths of several kilometers [2]^. This process initiated and propagated natural fracture networks, thereby creating high-permeability heat exchange reservoirs. Furthermore, numerical simulations and real-time monitoring were employed to optimize fracture geometry and connectivity [3]. In-depth investigations into the fracture propagation mechanisms and thermo-hydro-mechanical coupling processes not only improved thermal recovery efficiency and the operational lifespan of geothermal well fields, but also supported fracture network optimization and induced seismicity management [4]. These studies provided critical scientific support for the safe and controlled development of deep shale formation [5].

To clarify the dynamics of hydraulic fracture propagation in deep shale reservoirs during fracturing, researchers both domestically and internationally conducted experimental and numerical simulation studies on hydraulic fracture extension in these reservoirs [6]. The failure mechanisms of rocks under thermal stress were primarily categorized into two types: first, thermal expansion caused by increased temperature induced tensile stress within the rock, and when this stress exceeded the rock's tensile strength, thermal fracturing occurred; second, high-temperature rocks rapidly contracted upon cooling, concentrating internal thermal stress and inducing low-temperature cracking [7]. Tang et al [8]. demonstrated through experiments that the permeability of cores increased after high-temperature treatment, with a distinct permeability threshold for abrupt changes. They investigated the fracture threshold, thermal fracture evolution, and the factors influencing thermal fracturing in rocks of various lithologies. Guo et al [9]. evaluated the hydraulic fracturing behavior of large concrete specimens containing embedded cracks under high-temperature conditions through experiments and numerical simulations, revealing the effects of elevated temperature and crack networks on crack propagation. Rinaldi et al [10]. utilized the TOUGH2 simulator and found that thermal stress near the injection well was more significant than stress changes caused by water pressure, leading to the formation of new fractures or the opening of natural fractures. Zhu et al [11]. employed the discrete element method and concluded that thermal stress generated numerous microcracks perpendicular to the main fracture surface along the hydraulic fractures. The presence of these microcracks allowed fracturing fluid to penetrate the high-temperature rock mass along the fracture surface, leading to further microcrack generation. Colombero et al [12]. through microseismic monitoring in field experiments, observed that injecting low-temperature fluids into high-temperature deep shale created large-scale high-permeability zones. They concluded that thermal stress-induced fracturing was a key factor influencing the effectiveness of geothermal development.

In numerical simulation studies, Bujakowski et al [13]. applied the TOUGH2 software to evaluate heat extraction from thermal reservoirs and found that the volume and permeability of artificially fractured areas had a significant impact on the development effectiveness of EGS. Sun et al [14]. developed a THM coupling model to simulate the development effects of HDR, demonstrating that reducing flow resistance in EGS could improve heat extraction efficiency, although heat extraction decreased gradually over time. Li et al [15]. simulated the fracturing conditions and thermo-hydro coupling process of a well, showing that reservoir enhancement methods could enable the sustainable and stable development of geothermal energy. Lei et al [16]. conducted numerical simulations to analyze temperature and pressure variations in the fracture-porous media system during geothermal energy development.

The aforementioned studies identified shear failure as the primary failure mode in deep shale reservoirs with natural fractures, directly determining the fracture propagation patterns of hydraulic fractures. However, these studies neglected the influence of thermal stress generated by the injection of cold fluids during the modeling process on the formation of the fracture network [17,18]. Existing research primarily focused on the influence of natural fracture systems in deep shale reservoirs on fracture propagation and investigated the hydraulic fracture propagation patterns under circulation pump injection conditions. Nevertheless, few studies have examined the propagation patterns of multiple fractures during hydraulic fracturing. It is widely acknowledged that thermal stress-induced fracturing, caused by the injection of low-temperature fluids, significantly influences multi-fracture propagation. Understanding the multi-fracture propagation patterns in deep shale reservoirs is essential for enhancing the development effectiveness of deep shale [19].

Given the unclear understanding of the multi-fracture propagation patterns in deep shale reservoirs and the ambiguous influence of thermal stress-induced fracturing on these patterns, this study focused on deep shale reservoirs. By establishing a numerical simulation model, the study simulated and analyzed the multi-fracture propagation process within deep shale reservoirs, determining the multi-fracture propagation patterns under the influence of thermal stress fields. The technical approach is illustrated in Figure 1.

2. Method

During the hydraulic fracturing of deep shale reservoirs, the fracturing process was categorized into three primary stages: (1) Rock Fracturing Initiation: This stage involved the injection of high-pressure fracturing fluid, where the rock experienced failure due to the combined effects of elevated fluid pressure and thermal stress. (2) Fracture Propagation: As fracturing fluid continued to be injected, the fractures extended under the combined influence of high-pressure fluid, thermal stress, and the induced stress field. (3) Fracturing Fluid Leak-Off: Due to the porous nature of the reservoir rock, the fracturing fluid infiltrated the rock matrix during injection, altering the local pore pressure and temperature fields [20]. The governing equations for these interrelated physical processes were presented below.

2.1. Control Equations of Hydraulic Fracturing in Deep Shale Reservoirs

The theoretical foundation of hydraulic fracturing in deep shale reservoirs was based on the coupled thermal-fluid-solid interaction, which was described by the solid deformation equation. The governing equations of solid deformation primarily comprised the equilibrium equation, geometric relations, and constitutive equation, expressed as follows:

$${\sigma }_{ij,j}+b_i-\rho u_{i,tt}-\alpha u_{i,t}=0$$

(1)

$${\epsilon }_{i,j}={\displaystyle \frac{1}{2}}(u_{i,j}+u_{j,i})$$

(2)

$$d\sigma =Kd\epsilon$$

(3)

$${\displaystyle \frac{\partial {\sigma }_{ij}^{\text{'}}}{\partial x_j}}+b_i-{\rho }_s{\displaystyle \frac{d^2{u}_i}{d{t}^2}}-c{\displaystyle \frac{d{u}_i}{dt}}=0$$

(4)

During the hydraulic fracturing of deep shale reservoirs, both tensile and shear failures occurred. Therefore, the maximum tensile stress criterion and the Mohr–Coulomb criterion were employed to characterize tensile and shear failures, respectively. In this context, tensile failure was considered positive, while compressive failure was considered negative. The governing equations were as follows:

$$\left\{\begin{array}{l}{F}_1={\sigma }_1-f_{t0}\hfill \\ F_2=-{\sigma }_3+{\sigma }_1{\displaystyle \frac{1+\sin \theta }{1-\sin \theta }}-f_{c0}\hfill \end{array}\right.$$

(5)

When the rock experienced tensile failure (F₁ ≥ 0), the governing equation for the damage variable $Ds$ was expressed as follows:

$$D_{\mathrm{s}}=\left\{\begin{array}{l}0\epsilon <{\epsilon }_{\mathrm{t}0}\hfill \\ 1-{\displaystyle \frac{f_{\mathrm{tr}}}{\epsilon E_0}}{\epsilon }_{\mathrm{t}0}<\epsilon <{\epsilon }_{\mathrm{tu}}\hfill \\ 1{\epsilon }_{\mathrm{tu}}<\epsilon \hfill \end{array}\right.$$

(6)

Among ${\epsilon }_{\mathrm{tu}}=\eta {\epsilon }_{t0}$

When shear failure occurred in the rock (F₂ ≥ 0), the governing equation for the damage variable $Ds$ was expressed as follows:

$$D_{\mathrm{s}}=\left\{\begin{array}{l}0{\epsilon }_{\mathrm{c}0}<\epsilon \hfill \\ 1+{\displaystyle \frac{f_{\mathrm{cr}}}{\epsilon E_0}}\epsilon <{\epsilon }_{\mathrm{c}0}\hfill \end{array}\right.$$

(7)

The governing equation for the elastic modulus during fracturing was expressed as follows:

$$E=E_0\left(1-{\displaystyle \frac{n}{N}}\right)\left(1-D_{\mathrm{f}}\right)$$

(8)

During the fracturing process, the porosity of HDR varied with applied stress, while the permeability increased in response to the accumulation of rock damage. The governing equations were as follows:

$$\begin{array}{c}{\varphi }_{\mathrm{r}}=\left({\varphi }_{\mathrm{r}0}-{\varphi }_{\mathrm{re}}\right)\exp \left(-{\alpha }_{\varphi }{\sigma }_{\mathrm{eff}}\right)+{\varphi }_{\mathrm{re}}\\ K_{\mathrm{r}}=K_0{\left({\displaystyle \frac{{\varphi }_{\mathrm{r}}}{{\varphi }_{\mathrm{r}0}}}\right)}^3\exp \left({\alpha }_K{D}_{\mathrm{f}}\right)\end{array}$$

(10)

$$f\left(I_1,\sqrt{J_2}\right)=\sqrt{J_2}-\alpha I_1-k=0$$

(11)

If $\sigma 1$ , $\sigma 2$ and $\sigma 3$ represented the maximum, intermediate, and minimum principal stresses, respectively, then: $I1=\sigma 1+\sigma 2+\sigma 3$ , $J2=\left[{\left(\sigma 1-\sigma 2\right)}^2+{(\sigma 2-\sigma 3)}^2+{(\sigma 1-\sigma 3)}^2\right]/6$

The relationship between variables $\alpha$ and $\phi$ ,as well as between variables $k$ and $\phi$ ,was defined as follows:

$$\alpha ={\displaystyle \frac{\sin \phi }{\sqrt{9+3{\sin }^2\phi }}}$$

(11)

$$k={\displaystyle \frac{\sqrt{3}c\cos \phi }{\sqrt{3+{\sin }^2\phi }}}$$

(12)

The tensile and shear strains generated during the fracturing of deep shale reservoirs were as follows:

$${\epsilon }_{pt}={\displaystyle \frac{u_n-({\sigma }_t(t_0)A_c)/K_n}{L}}$$

(13)

$${\epsilon }_{ps}={\displaystyle \frac{u_s-(F_n\tan \varphi +c(t_0)A_c)/K_s}{L}}$$

(14)

Considering the impact of temperature variation on permeability, the temperature sensitivity coefficient was incorporated into the permeability evolution equation to achieve temperature-stress-seepage coupling. Consequently, the strain and permeability evolution equations under thermo-hydro-mechanical (THM) coupling conditions were derived as follows [21]:

$${\displaystyle \frac{k}{k_0}}={\displaystyle \frac{e^{M_T}\cdot T(1-{\phi }_0^{\gamma }\sqrt{-M_{\sigma }\Delta \sigma })}{2-{\phi }_0^{\gamma }\sqrt{-M_{\sigma }\Delta \sigma }-{\phi }_0}}{(1+{\displaystyle \frac{1-{\phi }_0}{{\phi }_0(1-{\phi }_0^{\gamma }\sqrt{-M_{\sigma }\Delta \sigma }}})}^3$$

(15)

2.2. Fluid Flow Equation in Cracks

The Navier–Stokes equations governing incompressible Newtonian fluids were restated as follows:

$$\rho {\displaystyle \frac{\mathrm{d}\overrightarrow{v}}{\mathrm{d}t}}=\rho \overrightarrow{F}+\mu {\nabla }^2\overrightarrow{v}-\nabla p$$

(16)

Neglecting inertia and body forces, the above equation was simplified as follows:

$$\nabla p=\mu {\nabla }^2\overrightarrow{v}$$

(18)

The governing equation for fluid flow within deep shale was obtained:

$$\nabla \cdot \left[{\displaystyle \frac{w^2}{12\mu }}\nabla p\right]={\displaystyle \frac{\partial w}{\partial t}}+{\displaystyle \frac{2C_{\mathrm{L}}}{\sqrt{t-\tau \left(x,y\right)}}}$$

(20)

$$D\left(w\right)={\displaystyle \frac{w^3}{12\mu }}$$

(21)

3. Numerical Simulation

3.1. Establishment of a Hydraulic Fracture Propagation Model for Deep Shale Reservoirs

To investigate the multi-fracture propagation behavior during the horizontal well fracturing of deep shale reservoirs, a two-dimensional hydraulic fracturing model was developed using ABAQUS software. The model represented a 150 m × 100 m reservoir with a single-stage, three-cluster configuration. The rock matrix was discretized using 48,000 CPE4P elements (pore pressure structural elements), while hydraulic fractures were modeled with 3,760 COH2D4P elements [24]. Detailed parameters are provided in Table 1.

3.2. Model Validity Verification

During the study, a 30 cm × 30 cm numerical simulation model was developed to replicate the experimental fracture morphology of the main plane. The model's validity was evaluated by comparing the variations in the hydraulic fracturing injection pressure curve. Based on the data reported in the literature, the maximum horizontal principal stress was set at 12.0 MPa, the minimum horizontal principal stress at 10.4 MPa, and the vertical stress at 15.0 MPa. The reservoir temperature was maintained at 100°C, and the injection flow rate was set to 25.0 mL/min.

Figure 3. Deep shale hydraulic fracturing model verification.

Figure 3. Deep shale hydraulic fracturing model verification.

By comparing the numerical simulation results of deep shale reservoirs with published experimental data, it was found that: first, hydraulic fractures primarily exhibited planar characteristics during initiation and propagation, with no significant bending observed. The fracture morphology obtained from the numerical simulation closely matched the experimental results. Second, the experiments indicated a rock fracture pressure of 16.3 MPa during hydraulic fracturing and an average fracture extension pressure of 11.2 MPa. The numerical simulation yielded a rock fracture pressure of 15.8 MPa and an average fracture extension pressure of 11.6 MPa. The errors between the simulated and experimental rock fracture pressures and fracture extension pressures were 3.07% and 3.57%, respectively, indicating minimal deviation between the model accuracy and the actual experimental results.

3.3. Analysis of Numerical Simulation Results

3.3.1. Effect of Reservoir Temperature on Hydraulic Fracture Expansion

Based on the hydraulic fracturing crack propagation model of deep shale reservoirs, reservoir temperatures were set to 120°C, 140°C, 160°C, 180°C, and 200°C, respectively. The dynamic propagation behavior of hydraulic fractures under single-fracture conditions was analyzed while maintaining a horizontal stress difference of 5.0 MPa. Figure 4 illustrates the fracture propagation morphology at different formation temperatures. The simulation results indicated that when the formation temperature was below 160°C, the hydraulic fracture consisted of one main fracture and one branch fracture [26]. However, at formation temperatures of 180°C and 200°C, the number of branch fractures increased significantly, forming a fracture network characterized by one main fracture and four branch fractures.

To further investigate the characteristics of hydraulic fracturing under varying formation temperatures, four injection pressure curves corresponding to different temperatures were extracted. By comparing the injection pressure variations, it was observed that as the formation temperature increased, the onset time of rock damage during hydraulic fracturing was significantly advanced. Specifically, rock damage occurred at approximately 320 s when the temperature was 140°C, advanced to 296 s at 160°C, further advanced to 261 s at 180°C, and reached 148 s at 200°C.

Figure 5. Comparison of injection pressure under different formation temperature conditions.

Figure 5. Comparison of injection pressure under different formation temperature conditions.

Table 2. Statistics of reservoir rock rupture time and pressure under different formation temperature conditions.

Formation Temperature/℃	Rock Fracture Time/s	Rock Fracture Pressure/MPa
140	320	65.8
160	396	63.4
180	261	57.4
200	148	51.2

Table 2. Statistics of reservoir rock rupture time and pressure under different formation temperature conditions.

Formation Temperature/℃	Rock Fracture Time/s	Rock Fracture Pressure/MPa
140	320	65.8
160	396	63.4
180	261	57.4
200	148	51.2

By comparing the fracture pressure and fracture extension pressure at different formation temperatures, it was found that when the formation temperature was 140°C, the reservoir rock fracture pressure was 65.8 MPa, and the average fracture extension pressure was 47.2 MPa. At 160°C, the fracture pressure decreased to 63.4 MPa, and the average extension pressure was 45.7 MPa. When the temperature rose to 180°C, the fracture pressure further declined to 57.4 MPa, with an average extension pressure of 43.1 MPa. At 200°C, the fracture pressure dropped to 51.2 MPa, and the average extension pressure was 41.1 MPa [27]. These results demonstrated that both fracture pressure and fracture extension pressure decreased to varying degrees as the formation temperature increased, with a significant reduction observed at 180°C. Using a commonly employed reservoir rock fracture pressure calculation formula, the theoretical fracture pressure under these conditions was approximately 66.4 MPa. The measured fracture pressures at 140°C and 160°C were close to this theoretical value, whereas a marked decline occurred at 180°C.

$${\mathrm{P}}_f=3{\sigma }_h-{\sigma }_H+S_t$$

(20)

The primary cause of this phenomenon was that, as the formation temperature increased, the injection of fracturing fluid into the formation caused rapid cooling of the rock, generating a low-temperature-induced thermal stress fracturing effect. This thermal stress fracturing led to a sharp decrease in rock fracture pressure. When the temperature difference between the reservoir rock and fracturing fluid exceeded a certain threshold, the low-temperature-induced thermal stress fracturing effect occurred. During the simulation, at a formation temperature of 180°C (corresponding to a temperature difference of 160°C between the rock and the fracturing fluid), a significant thermal stress fracturing effect was observed. As this temperature difference further increased, the injection pressure required for rock failure decreased, resulting in earlier rock fracture initiation. The thermal stress fracturing phenomenon, caused by the temperature difference between the reservoir and the fracturing fluid, also contributed to an increase in the number of branch fractures during hydraulic fracturing due to the reduction in rock fracture pressure. Therefore, when the formation temperature exceeded 180°C, the number of branch fractures increased significantly.

Figure 6. Statistical diagram of reservoir rock fracture time and pressure under different formation temperature conditions.

Figure 6. Statistical diagram of reservoir rock fracture time and pressure under different formation temperature conditions.

3.3.2. Effect of Horizontal Stress Difference on the Extension of Multiple Cracks

During the synchronous propagation of multiple hydraulic fractures, the "stress shadow" effect caused by the induced stress between fractures significantly influenced their simultaneous expansion. To date, both domestic and international researchers have extensively studied the "stress shadow" effect, analyzing how fracture spacing, horizontal stress difference, and injection parameters affect fracture propagation. However, these studies did not consider the impact of temperature field variations on fracture growth. In this section, twelve numerical simulation experiments were conducted to investigate the influence of horizontal stress difference under varying formation temperatures on the synchronous propagation of multiple hydraulic fractures. The formation temperatures were set to 140°C, 200°C, and 260°C, while the horizontal stress differences were set to 6.0 MPa, 8.0 MPa, 10.0 MPa, and 12.0 MPa. To eliminate the influence of changes in the minimum horizontal principal stress on fracture propagation, the minimum horizontal principal stress was held constant throughout the simulations; different horizontal stress differences were achieved by varying the maximum horizontal principal stress. A single-stage, three-cluster hydraulic fracturing case was used as an example, with a fracture cluster spacing of 15.0 m.

In the first set of numerical simulations, the morphology of hydraulic fracture propagation was analyzed at a formation temperature of 140°C with horizontal stress differences of 6.0 MPa, 8.0 MPa, 10.0 MPa, and 12.0 MPa, respectively. The results indicated that as the horizontal stress difference increased, the discrepancy between the propagation length of the middle fracture and that of the side fractures also increased. When the horizontal stress difference was 6.0 MPa, the three hydraulic fractures exhibited nearly free propagation, and the stress shadow effect exerted minimal inhibition on the middle fracture, whose propagation length was only 2.4% shorter than that of the side fractures. At 8.0 MPa, the propagation of the middle fracture was slightly constrained, with its length reduced by 8.2% relative to the side fractures [29]. When the horizontal stress difference increased to 10.0 MPa, the inhibition of the middle fracture propagation became significant, with its length reduced by 34.7% compared to the side fractures. At 12.0 MPa, the restriction on the middle fracture propagation intensified further, the stress shadow effect was enhanced, and the difference in fracture lengths increased to 51.6%.

Figure 7. Multi-fracture extension diagram under different horizontal stress differences at formation temperature of 140℃.

Figure 7. Multi-fracture extension diagram under different horizontal stress differences at formation temperature of 140℃.

Figure 8. Comparison of multiple crack extension lengths under different horizontal stress differences.

Figure 8. Comparison of multiple crack extension lengths under different horizontal stress differences.

Table 3. Comparison of multiple fracture extension lengths under different horizontal stress differences at formation temperature of 140°C.

Horizontal Stress Difference/MPa	Average Length of Edge Cracks/m	Length of Middle Cracks/m	Length Difference/%
6.0	84.6	82.6	2.4
8.0	87.3	80.1	8.2
10.0	89.1	58.2	34.7
12.0	89.7	43.4	51.6

Table 3. Comparison of multiple fracture extension lengths under different horizontal stress differences at formation temperature of 140°C.

Horizontal Stress Difference/MPa	Average Length of Edge Cracks/m	Length of Middle Cracks/m	Length Difference/%
6.0	84.6	82.6	2.4
8.0	87.3	80.1	8.2
10.0	89.1	58.2	34.7
12.0	89.7	43.4	51.6

When the formation temperature increased to 200°C, the hydraulic fracture propagation morphology under horizontal stress differences of 6.0 MPa, 8.0 MPa, 10.0 MPa, and 12.0 MPa is shown in Figure 9. As illustrated in Figure 9 and Table 4, the disparity in propagation lengths among the three hydraulic fractures was significantly reduced at 200°C. When the horizontal stress difference was 6.0 MPa, the propagation length of the edge fracture was 94.2 m, while that of the middle fracture was 92.8 m. The stress shadow effect exerted minimal inhibition on the middle fracture, whose length was 1.4% shorter than the edge fracture. At 8.0 MPa, the edge fracture propagated to 93.7 m, and the middle fracture reached 88.9 m, representing a 5.2% reduction relative to the edge fracture. When the horizontal stress difference increased to 10.0 MPa, the edge fracture extended to 92.3 m, whereas the middle fracture length decreased significantly to 72.6 m, showing a 21.3% reduction. At 12.0 MPa, the propagation of the middle fracture was further constrained; the edge fracture extended to 92.1 m, while the middle fracture shortened to 69.2 m. Consequently, the stress shadow effect intensified, and the difference in fracture lengths increased to 24.8%.

Figure 9. Multi-fracture extension diagram under different horizontal stress differences at formation temperature of 200℃.

Figure 9. Multi-fracture extension diagram under different horizontal stress differences at formation temperature of 200℃.

Figure 10. Comparison of the extension lengths of multiple cracks under different horizontal stress differences.

Figure 10. Comparison of the extension lengths of multiple cracks under different horizontal stress differences.

When the formation temperature increased to 260°C, the hydraulic fracture propagation morphology under horizontal stress differences of 6.0 MPa, 8.0 MPa, 10.0 MPa, and 12.0 MPa is presented in Figure 12. As shown in Figure 12 and Table 5, at this temperature, the difference in propagation lengths among the three hydraulic fractures was significantly reduced. When the horizontal stress difference was 6.0 MPa, the edge fracture propagated to 97.3 m, and the middle fracture extended to 96.8 m. The stress shadow effect exerted minimal inhibition on the middle fracture, whose length was 0.51% shorter than the edge fracture. At 8.0 MPa, the edge fracture reached 96.9 m, while the middle fracture extended to 95.4 m, representing a 1.44% reduction relative to the edge fracture. When the horizontal stress difference increased to 10.0 MPa, the edge fracture propagated to 94.3 m, and the middle fracture to 91.5 m. The degree of restriction on the middle fracture’s extension increased slightly, with the middle fracture length being 2.97% shorter than the edge fracture. At 12.0 MPa, the propagation of the middle fracture was further constrained; the edge fracture extended to 93.7 m, while the middle fracture shortened to 90.2 m. Consequently, the stress shadow effect intensified, and the difference in fracture lengths increased to 3.74% [30].

Figure 11. Comparison of extension lengths of multiple cracks under different horizontal stress differences.

Figure 11. Comparison of extension lengths of multiple cracks under different horizontal stress differences.

A comparison of fracture propagation morphology at different formation temperatures revealed that the extension length of hydraulic fractures increased with rising formation temperatures. At 140°C, as formation pressure increased, the stress shadow effect on the middle fracture intensified, leading to a gradual reduction in its extension length. However, at 200°C and 260°C, the disparity between the extension lengths of the middle and edge fractures was significantly reduced.

As shown in Figure 13, when the formation temperature was 140°C, the difference in extension lengths between edge and middle fractures increased significantly once the formation pressure exceeded 8 MPa. When the horizontal stress difference reached 12 MPa, this difference reached 51.6%. Increasing the formation temperature effectively mitigated the impact of the stress shadow effect on fracture propagation. Specifically, at formation temperatures of 200°C and 260°C with a horizontal stress difference of 12 MPa, the difference in extension lengths between edge and middle fractures was reduced to 24.8% and 3.74%, respectively.

To further elucidate the influence of formation temperature on fracture propagation, the injection pressure curves were extracted, providing insights into hydraulic fracturing pressures under various conditions. The results indicated that increasing the formation temperature significantly reduced the rock fracture pressure during hydraulic fracturing. A comparison of rock fracture pressures at 140°C and 260°C revealed a 22.9% reduction in fracture pressure with rising temperature. This reduction in fracture pressure also promoted the formation of branch fractures. The distribution of hydraulic fractures at different formation temperatures demonstrated that the number of branch fractures increased substantially with rising temperature. At 140°C, hydraulic fracture propagation was predominantly characterized by the expansion of primary fractures. In contrast, at 200°C and 260°C, the number of branch fractures increased markedly during propagation. These findings suggest that a greater temperature difference between the formation and the fracturing fluid enhanced the low-temperature-induced thermal stress fracturing effect. This effect not only reduced fracture pressure but also increased fracture complexity, diminished the stress shadow effect’s impact on fracture propagation, and promoted the development of complex fracture networks [31].

Figure 14. Comparison of hydraulic fracturing fracture pressure under different formation temperatures and different horizontal stress differences.

Figure 14. Comparison of hydraulic fracturing fracture pressure under different formation temperatures and different horizontal stress differences.

3.3.3. Effect of Cluster Spacing on Multi-Crack Extension

To further investigate the impact of cluster spacing on the propagation of multiple fractures, the formation temperatures were set to 180°C and 240°C, with a horizontal stress difference of 14.0 MPa. Numerical simulations were conducted to examine the dynamic propagation of multiple hydraulic fractures at cluster spacings of 5.0 m, 8.0 m, 10.0 m, and 13.0 m.

At a formation temperature of 180°C, as the cluster spacing decreased to 10.0 m, the propagation of the middle fracture became restricted, with a difference of 1.38 m in extension length between the edge and middle fractures, resulting in a 0.97% reduction in the middle fracture's extension length. When the cluster spacing was further reduced to 8.0 m, the difference in extension length between the middle and edge fractures increased, leading to a 1.67% reduction in the middle fracture's extension length. At a cluster spacing of 5.0 m, this difference further expanded to 3.47%.

Conversely, increasing the cluster spacing reduced the number of hydraulic fracture branches. Considering both the extension length of hydraulic fractures and the number of branch fractures, the optimal cluster spacing under these conditions was determined to be 8.0 m.

Figure 15. Hydraulic fracture propagation morphology at formation temperature of 180℃ and different cluster spacing.

Figure 15. Hydraulic fracture propagation morphology at formation temperature of 180℃ and different cluster spacing.

Figure 16. Comparison of hydraulic fracture extension lengths at different cluster spacings.

Figure 16. Comparison of hydraulic fracture extension lengths at different cluster spacings.

Table 6. Comparison of hydraulic fracture extension lengths at different cluster spacings at formation temperature of 180°C.

Cluster Spacing/m		Average Length of Edge Cracks/m	Length of Middle Cracks/m	Length Difference/%
5.0	96.3		92.9	3.47
8.0	95.8		94.2	1.67
10.0	94.6		93.7	0.97
13.0	94.1		93.8	0.31

Table 6. Comparison of hydraulic fracture extension lengths at different cluster spacings at formation temperature of 180°C.

Cluster Spacing/m		Average Length of Edge Cracks/m	Length of Middle Cracks/m	Length Difference/%
5.0	96.3		92.9	3.47
8.0	95.8		94.2	1.67
10.0	94.6		93.7	0.97
13.0	94.1		93.8	0.31

At a formation temperature of 240°C, the dynamic simulation results of multiple hydraulic fracture propagation under cluster spacings of 5.0 m, 8.0 m, 10.0 m, and 13.0 m are presented in Figure 17. The results indicated that hydraulic fractures generally exhibited unrestricted propagation under different cluster spacing conditions. Only when the cluster spacing was 5.0 m did the propagation of the middle fracture experience noticeable restriction.

As shown in Table 7, when the cluster spacing was 5.0 m, the extension length of the edge fracture increased by 1.22% compared to the middle fracture. At a cluster spacing of 8.0 m, the extension lengths of the edge and middle fractures were nearly identical, with a length difference of only 0.83%. Further increasing the cluster spacing to 13.0 m reduced the difference between the edge and middle fractures to 0.42%. Based on these results, the optimal cluster spacing under the given conditions was determined to be approximately 5.0 m.

Figure 18. Comparison of hydraulic fracture extension lengths with different cluster spacings when the formation temperature is 240°C.

Figure 18. Comparison of hydraulic fracture extension lengths with different cluster spacings when the formation temperature is 240°C.

By comparing the multi-fracture propagation patterns under various formation temperature conditions, it was observed that, under the same cluster spacing, higher formation temperatures led to the formation of more branch fractures and increased hydraulic fracture propagation lengths. Additionally, higher formation temperatures further reduced the optimal cluster spacing, facilitated closer spacing between fractures during hydraulic fracturing, and mitigated the impact of the stress shadow effect on multiple fracture propagation.

3.3.4. Effect of Rock Thermal Conductivity on the Extension of Multiple Fractures

To investigate the impact of rock thermal conductivity on the dynamic propagation of multiple hydraulic fractures, the temperature of the deep shale reservoir was set to 200°C, with a horizontal stress difference of 10.0 MPa and a cluster spacing of 4.0 m. The rock thermal conductivity was varied at 1.2 W/(m·K), 1.3 W/(m·K), 1.4 W/(m·K), 1.5 W/(m·K), and 1.6 W/(m·K). Fure 19 presents the simulation results for the propagation of multiple hydraulic fractures under these different thermal conductivity conditions.

Figure 19. Dynamic simulation of multiple hydraulic fracture extension under different rock thermal conductivity conditions.

Figure 19. Dynamic simulation of multiple hydraulic fracture extension under different rock thermal conductivity conditions.

The simulation results indicated that as the thermal conductivity of the rock increased, the difference in extension length between the middle crack and the edge crack gradually decreased. At a thermal conductivity of 1.2 W/(m·K), the hydraulic fracture exhibited a pronounced internal inclusion phenomenon during propagation. This phenomenon was significantly reduced when the thermal conductivity increased to 1.3 W/(m·K). Furthermore, at a thermal conductivity of 1.5 W/(m·K), the three fractures propagated uniformly, accompanied by a substantial increase in the number of branch fractures.

Figure 20. Comparison of hydraulic fracture extension length under different rock thermal conductivity.

Figure 20. Comparison of hydraulic fracture extension length under different rock thermal conductivity.

Table 8. Comparison of hydraulic fracture extension length under different rock thermal conductivity.

Thermal Conductivity of Rock/ W/(m.K)	Average Length of Edge Cracks/m	Length of Middle Cracks/m	Length Difference/%
1.2	96.3	76.1	20.9
1.3	96.1	83.4	13.2
1.4	95.4	90.7	4.9
1.5	95.1	92.4	2.8
1.6	94.8	94.1	0.7

Table 8. Comparison of hydraulic fracture extension length under different rock thermal conductivity.

Thermal Conductivity of Rock/ W/(m.K)	Average Length of Edge Cracks/m	Length of Middle Cracks/m	Length Difference/%
1.2	96.3	76.1	20.9
1.3	96.1	83.4	13.2
1.4	95.4	90.7	4.9
1.5	95.1	92.4	2.8
1.6	94.8	94.1	0.7

Figure 21 illustrates the variation in injection pressure under different rock thermal conductivities. The simulation results revealed that as the thermal conductivity of the rock increased, the hydraulic fracturing injection pressure gradually decreased, and the rock fracture occurred earlier. This observation indicated that higher rock thermal conductivity enhanced the thermal stress fracturing effect induced by the injection of low-temperature fluids while mitigating the "stress shadow" effect during the propagation of multiple fractures.

4. Discussion and Analysis

To further elucidate the variation in reservoir rock properties during the fracturing of HDR reservoirs, a functional relationship between the degree of reservoir rock damage and permeability was established based on damage mechanics theory. Specifically, Equation 42 provided the calculation formula for determining the degree of reservoir rock damage.

$$D=1-{\displaystyle \frac{E^{\text{'}}}{E_0}}$$

(21)

During the research, the physical properties of reservoir rocks under single-fracture expansion conditions were initially analyzed. Figure 22 and Figure 23 illustrated the variations in Young's modulus and reservoir permeability of the rocks during hydraulic fracturing at different formation temperatures. The analysis of Young's modulus degradation under single-fracture conditions revealed that as the formation temperature increased, the degradation area of Young's modulus expanded significantly. Notably, when the formation temperature exceeded 180°C (specifically reaching 200°C), the degradation zone of Young's modulus exhibited a marked increase. Furthermore, Figure 23 demonstrated that due to reservoir damage, the area of the permeability enhancement zone also expanded significantly.

Figure 24 illustrated that at a formation temperature of 140℃, the average reservoir permeability exhibited a distinct three-peak pattern. The permeability was highest at the main fracture due to its extension and expansion, with a maximum value of 27.4 mD. As the distance from the main fracture increased, the reservoir permeability gradually decreased, showing a sharp decline approximately 12 m away from the main fracture. This decline resulted primarily from compressive deformation of the reservoir rock caused by the fracture surface opening during the formation of the main fracture.

Beyond the compression zone, permeability slightly increased due to damage induced in the distal rock by the extension of the main fracture and its branches, before eventually returning to the initial permeability value. When the formation temperature increased to 160℃, the reduction in permeability caused by fracture surface compression was significantly mitigated. This mitigation was mainly attributed to low-temperature-induced thermal stress fracturing, where rock damage from thermal stress compensated for the permeability reduction caused by compression. At formation temperatures of 180℃ and 200℃, the reservoir permeability decreased gradually from the main fracture to the initial value, with permeability reductions caused by rock compression being fully offset by thermal stress fracturing-induced damage.

5. Conclusions

Thermal stress damage caused by low-temperature fluid injection significantly influenced the propagation of multiple fractures during hydraulic fracturing. Based on linear elastic fracture mechanics theory, this study employed finite element numerical simulations to analyze factors affecting the propagation of multiple fractures in deep shale reservoirs, leading to the following conclusions:

(1) Thermal stress damage induced by low-temperature fluid injection effectively mitigated the restriction imposed by the "stress shadow" effect on the propagation of intermediate fractures, facilitating balanced initiation and propagation of multiple fractures. Additionally, low-temperature-induced thermal stress fracturing significantly reduced both rock fracture pressure and fracture propagation pressure during hydraulic fracturing.

(2) Higher reservoir temperatures and greater rock thermal conductivity intensified the low-temperature-induced thermal stress fracturing effect during hydraulic fracturing, thereby diminishing the influence of horizontal stress differences and cluster spacing on multiple fracture propagation.

(3) With increasing formation temperature, the degradation zone of Young's modulus in the reservoir rock expanded markedly during hydraulic fracturing, indicating that the stress disturbance zone enlarged as formation temperature increased

Nomenclature		Dimensionless Quantities are Marked with(*)
${\sigma }_{ij}$	Cauchy tensor(N/m²)	${\alpha }_T$	expansion coefficient (°C−1)
$b_i$	body force(N/m³)	$\overrightarrow{F}$	force per unit mass(m/s²)
$\rho$	density(kg/m³)	$p$	fluid pressure(N/m²)
$\alpha$	damping coefficient(N·s/m)	$q_{\mathrm{L}}$	fluid loss rate (m2/s)
$u_i$	displacement(m)	$C_{\mathrm{L}}$	filter coefficient (m/s)
${\epsilon }_{ij}$	strain(*)	$D$	fracturing damage degree (*)
$w$	weight function(*)	$E^{\text{'}}$	residual Young's modulus (Pa)
$J$	jacobian value (*)	${\sigma }_{1,2,3}$	Maximum, intermediate, minimum principal stresses (Pa)
$\alpha$	drucker-Prager criterion(Pa)	$\phi$	internal friction angle and cohesion (N/m²)
${\rho }_b$	bulk density of stone(kg/m³)	$\alpha$$k$	drucker-prager criterion material constants(Pa)
$D_f$	fatigue damage (*)	${\rho }_b$	bulk density of stone(kg/m³)
$E$	elastic modulus(Pa)	$\mu$	fluid viscosity(Pa·s)
$f$	frequency(Hz)	$\nabla p$	pore pressure(N/m²)
$f_{c0}$	compressive strength(Pa)	$D_S$	damage variables(*)
$f_{t0}$	tensile strength(Pa)	$E_0$	initial elastic modulus(Pa)
$K_0$	initial permeability(m2)	$f_{cr}$	compression strength(Pa)
$\epsilon$	strain(*)	$f_{tr}$	tensile Strength(Pa)
${\epsilon }_{t0}$	tensile principal strain(*)	$F_1{F}_2$	state Function(Pa)
${\varphi }_r$	rock porosity(%)	$K_r$	rock permeability(m2)
${\varphi }_{re}$	residual porosity(%)	${\epsilon }_{c0}$	compressive principal strain(*)
${\alpha }_k$	permeability variation coefficient(*)	${\epsilon }_{tu}$	ultimate tensile strain(*)
${\varphi }_{r0}$	Initial porosity(%)	${\varphi }_{r0}$	Initial porosity(%)
${\alpha }_{\varphi }$	influence coefficient(Pa-1)	${\sigma }_{eff}$	effective stress(Pa)