Introduction

Deep shale gas is taken as an important alternative field for unconventional oil and gas exploration and development, and how to form a large-scale production capacity is very important for China to guarantee energy security (Wei et al. 2022). Even if horizontal well staged fracturing technology is currently the main means to improve the recovery rate of deep shale gas, due to the three high characteristics of the deep shale formation environment (high temperature, high in-situ stress, and high stress difference), the hydraulic fracturing effect is not good and the productivity decreases rapidly after fracturing (He et al. 2023; Yijin 2014). Therefore, in order to find out the reasons for the poor effect of hydraulic fracturing in deep shale, it is necessary to understand the propagation rule of hydraulic fractures in deep shale.

However, conventional laboratory tests are difficult to meet the real shale reservoir environment (Qian et al. 2020), and most of the studies on numerical simulation focus on the fracturing technique and the interaction mechanism between hydraulic and natural fractures, both of which do not consider the temperature effect. The purpose of this paper is to study the propagation law of hydraulic fractures under temperature effect by means of numerical simulation, but some research results without considering temperature effect still have certain reference significance. Boone and Ingraffea (1990) first used the finite element method (FEM) to approximately simulate the fluid flow in the fracture and to determine the size of the fractures based on a KGD-like model. Hunsweck et al. (2013) used linear elastic fracture mechanics (LEFM) as a criterion for fracture propagation. Both of them finally simulate the propagation behavior of hydraulic fractures in impermeable elastic media. Carrier and Granet (2012) used the cohesive zone finite element model (CZM) to simulate the propagation behavior of hydraulic fractures in permeable media, and noted that the propagation path of hydraulic fractures has a strong correlation with the permeability of the matrix. Similarly, Yao et al. (2015) used the cohesive element to simulate the crack propagation process and noted that the brittleness index and fracture opening were inversely correlated with the elastic modulus. Settgast et al. (2017) combined FEM with finite volume method (FVM) to simulate the field-scale hydraulic fracture propagation behavior. Wang et al. (2012) developed a 3D hydraulic fracturing model and noted that the fracture height will decrease with high interfacial tensile strength, low elastic modulus, and in-situ stress. Based on the hydraulic fracturing parameters, Wang et al. (2017) optimized the perforation cluster distance, while Tang et al. (2023) and Zhu et al. (2023). adjusted the temporary plugging process parameters to maximize the fracturing volume. Some scholars have simulated the propagation process of hydraulic fractures based on extended finite element method (XFEM) (Gordeliy and Peirce 2013; Roth et al. 2020; Zeng et al. 2019), boundary element method (BEM) (Hattori et al. 2017; Wu and Olson 2014), discrete element method (DEM) (Hamidi and Mortazavi 2014; Zhang et al. 2016) and phase field method (PFM) (Mikelić et al. 2019; Wheeler et al. 2014; Yoshioka and Bourdin 2016), and some results with engineering guiding significance have been obtained.

The simulation of deep shale hydraulic fracturing process is actually a typical complex multi-field coupling problem involving fluid-rock mass heat transfer, rock deformation, rock damage and fluid flow. Most studies have shown that the change of temperature has a significant effect on the overall morphology of hydraulic fractures. Especially under high temperature conditions, the thermal shock effect of low temperature fracturing fluid will induce microfracture formation in openhole section (Guo et al. 2022; Wu et al. 2018), which not only hinders the formation of complex fractures, but also greatly reduces the conductivity of hydraulic fractures. In order to further study the causes of these differences and the propagation rule of hydraulic fractures, more dynamic evolution processes of hydraulic fractures and rock field variables need to be obtained. Unfortunately, on the one hand, the current test conditions are not able to monitor these important processes. On the other hand, compared to the typical hydraulic fracturing physical model test, the expense of the test under high temperature and pressure is significantly higher, and it does not seem feasible to repeat the test in large numbers. Therefore, with the help of numerical simulation technology, the initiation as well as the propagation of hydraulic fractures can be carried out in two steps (thermal shock and hydraulic fracture propagation). Based on the CZM, a numerical model of hydraulic fracture propagation considering the thermal shock effect is established. The effects of different temperatures, thermal shock rates, and times on the scale of thermal fractures are analyzed. The evolution rule of hydraulic fracture propagation behavior and characteristic parameters under thermal shock are revealed. This study can provide some inspiration for the development of deep shale by improving the complexity of hydraulic fractures.

Simulation method and model setup

In order to simulate the process of borehole surrounding rock cracking and subsequent hydraulic fracture propagation after thermal effect. A numerical calculation model is established based on the CZM theory. The model first calculates the damage zone of the borehole surrounding rock after the thermal shock, then reconstructs the damage zone and calculates the propagation path of the hydraulic fracture. In this model, the rock deformation is simulated by the physical property parameters given to the basic elements, and the rock cracking caused by thermal stress and fluid pressure is simulated by the cementation surface between the basic elements. Finally, the thermo-solid coupling, fluid–solid coupling and deformation fracture process in the overall model can be solved by the corresponding governing equations (Wu et al. 2023a, b).

Governing equations

Stress balance equation

Considering the influence of temperature and pore pressure changes on rock deformation, there is the following Navier-type stress equilibrium equation (Zhang et al. 2021):

$$G{u}_{i,jj}+\frac{G}{1-2v}{u}_{j,ji}-\alpha {p}_{{\text{p}}}-K{\alpha }_{{\text{T}}}{T}_{i}+{f}_{i}=0$$
(1)

where \(G\) is the shear modulus, \(u\) is the displacement, \(\nu\) is the Poisson’s ratio, \(\alpha\) is the Biot’s coefficient, \({p}_{{\text{p}}}\) is the pore pressure, \(K\) is the volume modulus, \({\alpha }_{{\text{T}}}\) is the thermal expansion coefficient of rock, \(T\) is the absolute temperature, and \(f\) is the body force.

Heat transfer equation

The thermal shock problem of low temperature fluid can be considered as the problem of thermos-solid coupling in the process of rapid temperature reduction. In the model, it is assumed that the rock is in a local thermal equilibrium state during the rapid cooling process and only heat conduction is considered. Then there is the following heat transfer equation (Ogata et al. 2018):

$$(1-\varphi ){\rho }_\text{r}{C}_\text{r}\frac{\partial T}{\partial t}=\nabla \cdot \left((1-\varphi ){\lambda }_{{\text{r}}}\nabla T\right)+{Q}_{{\text{h}}}$$
(2)

where \(\varphi\) is the porosity of rock, \({\rho }_\text{r}\) is the density of rock, \({C}_\text{r}\) is the specific heat capacity of rock, \({\lambda }_{{\text{r}}}\) is the thermal conductivity of rock, and \({Q}_{{\text{h}}}\) is the heat source.

Criteria for fracture initiation and propagation

The finite element model used in this paper is to realize the initiation and propagation of fractures (thermal fractures and hydraulic fractures) by inserting global cohesive elements at the boundaries of all elements. The linear-elastic relationship of the cohesive element simulates the initial deformation of the fracture (Guo et al. 2015; Wu et al. 2023a, b; Wu et al. 2020):

$$\sigma =\left[\begin{array}{c}{\sigma }_\text{n}\\ {\sigma }_\text{s}\end{array}\right]$$
(3)

where \(\sigma\), \({\sigma }_{n}\) and \({\sigma }_{s}\) are stress matrix, normal stress and tangential stress. The above stress is obtained by stiffness and strain, and the corresponding strain is as follows:

$${\varepsilon }_{n}=\frac{{d}_{{\text{n}}}}{{T}_{0}},{\varepsilon }_{s}=\frac{{d}_{{\text{s}}}}{{T}_{0}}$$
(4)

where \({d}_{{\text{n}}},{d}_{{\text{s}}}\) and \({T}_{0}\) are normal displacement, tangential displacement and constitutive thickness.

The maximum stress criterion was used to simulate fracture initiation when the fracture started to expand:

$${\left\{\frac{\langle {\sigma }_\text{n}\rangle }{{\sigma }_\text{n}^{0}}\right\}}^{2}+{\left\{\frac{{\sigma }_\text{s}}{{\sigma }_\text{s}^{0}}\right\}}^{2}=\lambda$$
(5)

where \({\sigma }_{n}^{0}\) and \({\sigma }_{s}^{0}\) are the critical stresses in different directions; \(\langle \cdot \rangle\) is the tensile stress only, and \(\lambda\) is a constant between 1 and 1.05. The degree of weakening of the stiffness of the cohesive element when the fracture opens can be expressed by the following equation (Dahi Taleghani et al. 2018, 2016):

$$E=(1-d){E}^{0}$$
$$d=\frac{{\delta }_{m}^{f}\left({\delta }_{m}^{{\text{max}}}-{\delta }_{m}^{0}\right)}{{\delta }_{m}^{{\text{max}}}\left({\delta }_{m}^{f}-{\delta }_{m}^{0}\right)}$$
(6)

where \({E}^{0}\) is the initial elastic modulus, \(E\) is the elastic modulus after damage, \(d\) is the degree of stiffness degradation (SDEG), \({\delta }_{m}^{f}\) is the complete damage displacement, \({\delta }_{m}^{{\text{max}}}\) is the maximum load displacement, and \({\delta }_{m}^{0}\) is the initial damage displacement.

In addition, in the process of hydraulic fracture propagation, assuming that the fracturing medium is incompressible Newtonian fluid, the flow of fluid in the cohesive unit is as follows:

$$q=\frac{{w}^{3}}{12\mu }\nabla p$$
(7)

where \(w\) is the thickness of the element, \(\mu\) is the fluid viscosity, \(p\) is the fluid pressure. The process of fracture propagation is also accompanied by fluid leakage, and the rate of leakage can be expressed as follows:

$${q}_{{\text{l}}}={c}_{{\text{l}}}\left({p}_{{\text{f}}}-{p}_{{\text{p}}}\right)$$
(8)

where \({q}_{l}\) is the local leakage velocity per unit fracture surface area, \({p}_{p}\) is the pore pressure, \({c}_{l}\) is the leakage coefficient.

The fluid flow in reservoir rock follows Darcy 's law:

$${q}_{m}=-\frac{k}{\mu }\nabla {p}_{{\text{p}}}$$
(9)

where \(k\) is the permeability tensor, \({q}_{m}\) is the velocity vector of fluid flow.

Since the numerical modeling is divided into two stages of calculations: (1) the establishment and calculation of the thermal shock model of the borehole surrounding rock considering the temperature of the rock mass, the cooling rate and the number of thermal shock; (2) establishment and calculation of hydraulic fracturing model for reconstructing thermal shock damage zone (identify the damaged area → mark the damaged area → reconstruct the new model → eliminate the damaged area). In order to simplify the modeling and calculation process, this paper developed codes that automate modeling and computation and nested them into Abaqus software to realize batch calculation and processing. The modeling and calculation process of the whole model are shown in Fig. 1.

Fig. 1
figure 1

Modeling workflow of the fracture’s initiation and propagation

Full size image

Initial model setup

Firstly, in order to study the thermal shock effect of low temperature fluid, according to the geometric characteristics of laboratory test samples, the two-dimensional plane in the direction of horizontal wellbore is selected, as shown in Fig. 2. The model is a square with sides of 300 mm, and the diameter of the borehole is 30 mm. The whole model is in a constant temperature state. The upper and lower boundaries simulate the in-situ stress state of the model by applying loads. The inner wall of the borehole is set as the temperature boundary condition (cooling heat conduction process), where an initial temperature of the borehole is set as well as the temperature of the borehole after it has cooled down for a specified calculation time. Secondly, in the process of simulating hydraulic fracturing, the in-situ stress state is realized by the "in-situ stress balance" module in ABAQUS. The injection point of fracturing fluid is the boundary of the whole borehole wall, and the pumping displacement is constant. The rapid propagation of hydraulic fractures is driven by fluid pressure, and the whole simulation time is 100 s. It should be noted that the damage zone of the borehole surrounding rock in the first stage can be used as the initial damage point for the fracturing fluid injection in the second stage. It can not only eliminate the step that the traditional model needs to artificially define the initial damage point, but also further promote the convergence of the model.

Fig. 2
figure 2

Numerical model diagram

Full size image

Generally, the parameters required in the simulation are divided into formation parameters and engineering parameters, with the former generally including the stress state, formation temperature, and physical and mechanical properties of the rock, and the latter generally including the fracturing fluid pumping displacement and the viscosity of the fracturing fluid, etc. Where the process of initial stress equilibrium is solved by employing the initial equilibrium enhancement program in Abaqus/Standard (Wu et al. 2023a, b, 2021a, 2021b). In this paper, the simulation is conducted at the scale of laboratory tests, so the simulation parameters are mainly consistent with the laboratory test parameters. Considering the low porosity and permeability of shale reservoirs, the initial porosity ratio is set to 0.001, and the permeability is reduced from the set value to impermeable, as shown in Table 1.

Table 1 Main parameters of numerical model
Full size table

Numerical experiment results

Cracking characteristics of borehole surrounding rock under thermal shock

The thermal shock effect is obvious for the promotion of thermal fracture formation in rock. Therefore, this paper first analyzes the process of heat transfer between low temperature fluid and high temperature rock mass in the process of hydraulic fracturing of deep shale. The influencing factors mainly include rock temperature, thermal shock rate and thermal shock times, which correspond to the formation temperature, fracturing fluid type and fracturing technique of the actual construction site.

Effect of thermal shock temperature

The distribution characteristics of thermal fractures in borehole surrounding rock with shale temperature of 100, 200, and 300 °C, open hole section cooling to 20 °C, and cooling time of 25 s are shown in Fig. 3. According to Fig. 3, due to the rapid decrease in temperature of the borehole surrounding rock, some thermal fractures of different scales are formed under the action of thermal stress, which are evenly distributed along the circumference of the borehole. Under the condition of 100 ℃, the thermal fractures are sparser, and the opening degree is smaller, and the extension distance is shorter; under the condition of 200 ℃, the thermal fractures are gradually dense, and the opening degree is obviously increased, and they can be extended for a certain distance; under the condition of 300 ℃, the tendency to change the characteristics of the thermal fractures is more obvious. This shows that when the temperature of the thermal shock fluid is constant, the higher the shale temperature, the stronger the tensile stress generated on the surface of the borehole surrounding rock, the larger the thermal fracture occurs, and the more obvious the thermal shock effect.

Fig. 3
figure 3

Distribution of thermal fractures in borehole surrounding rock at different temperature

Full size image

In order to further analyze the distribution characteristics of thermal fractures under different temperature conditions, the number, width and length of thermal fractures with SDEG (cohesive element stiffness damage value) of 0.98 ~ 1 were counted, as shown in Fig. 4. According to Fig. 4, at 100 °C, 29 thermal fractures were generated in the borehole surrounding rock; at 200 °C, 35 thermal fractures were generated. At 300 °C, 57 thermal fractures were generated. With the increase of temperature, the number of hot seams shows an upward trend. Interestingly, no matter what kind of temperature conditions, there are two partitions in the thermal fracture scale, which are defined as small-scale thermal fracture zone and large-scale thermal fracture zone, and the maximum fracture width thermal fracture is the optimal thermal fracture. At 100 °C, the scale zoning characteristics of thermal fractures are the most obvious. Small-scale thermal fractures are mainly distributed in the range of 2–8.5 μm in width and 2–4 μm in length. Large-scale thermal fractures ranged from 7 to 11 μm in width and 5.5–7 μm in length. The width and length of the optimal thermal fracture are 10.6 μm and 6.1 μm. At 200 °C, the scale zoning characteristics of thermal fractures are slightly weakened. Small-scale thermal fractures are mainly distributed in the range of 3–15 μm width and 6–13 μm length. Large-scale thermal fractures ranged from 15 to 30 μm in width and 13–19 μm in length. The width and length of the optimal thermal fracture are 29.9 μm and 18.9 μm. At 300 °C, thermal fractures gradually gather at the junction of the partition. Small-scale thermal fractures ranged from 2 to 15 μm in width and 8–20 μm in length. Large-scale thermal fractures ranged from 16 to 37 μm in width and 18–27 μm in length. The width and length of the optimal thermal fracture are 36.2 μm and 26.8 μm. With the increase of temperature, the width and length of thermal fractures become larger. The width and length of the optimal thermal fracture at 300 °C are about 3.5 times that of 100 °C. Under the driving of fluid pressure, hydraulic fractures are more likely to propagate directly along these large-scale thermal fractures.

Fig. 4
figure 4

Thermal fractures characteristic parameters at different temperature

Full size image

In general, the thermal fractures in the borehole will not be fully opened by hydraulic fractures. Due to the overall increase in the scale of the thermal fracture, the hydraulic fracture certainly tends to open the thermal fracture with a relatively large fracture width, so it can be considered that the large-scale thermal fracture is a fracture with opening potential. Figure 4d shows the total number of thermal fractures and the stacking histogram of thermal fractures with opening potential. According to Fig. 4d, 8, 19, and 28 thermal fractures with open potential were produced at 100 °C, 200 °C, and 300 °C, accounting for 27.5%, 54.2%, and 49.1%, respectively. It can be seen that the thermal shock with a large temperature difference is the material basis for the formation of complex fractures by hydraulic fracturing.

Effect of thermal shock rate

Due to the shale at 200 ℃ to form the highest percentage of thermal fractures with opening potential, this section sets the shale temperature to 200 °C, the open hole section to 20 °C, and the cooling time to 25 s, 50 s, 100 s. To clarify the influence of thermal shock rate on thermal fracture characteristic parameters. The distribution of thermal fractures in borehole surrounding rock under different thermal shock rates is shown in Fig. 5. According to Fig. 5, compared with the thermal shock duration of 25 s, when the duration is 50 s, the distance between the thermal fractures begins to increase, but the width and length begin to increase; when the duration is 100 s, the thermal fractures become sparse, and the width and length increase further. With the increase in thermal shock duration, the thermal shock rate slows down, the density of thermal fractures gradually decreases, but the scale gradually increases. This is mainly because the thermal shock time becomes longer, which means that the temperature of the borehole surrounding rock decreases slowly, and the heat exchange efficiency between the fluid and the rock becomes lower. The lower heat exchange efficiency will inhibit the initiation of rock, resulting in a decrease in the number of thermal fractures. However, the tensile stress caused by thermal shock will gradually increase, so the length and width will gradually increase.

Fig. 5
figure 5

Distribution of thermal fractures in borehole surrounding rock at different thermal shock rate

Full size image

The characteristic parameters of thermal fractures under different thermal shock rates are also counted, as shown in Fig. 6. According to Fig. 6, when the thermal shock duration is 50 s, 25 thermal fractures are generated in the borehole surrounding rock; when the duration is 100 s, the number of thermal fractures is reduced to 22. With the increase in thermal shock duration, the number of thermal fractures shows a slow downward trend. In addition, the zoning characteristics of the thermal fracture scale under different thermal shock rates are still obvious. At 50 s, small-scale thermal fractures ranged from 2 to 23 μm in width and 5–20 μm in length. Large-scale thermal fractures ranged from 27 to 46 μm in width and 17.5–32.5 μm in length. The width and length of the optimal thermal fracture are 44.3 μm and 26.5 μm. At 100 s, small-scale thermal fractures ranged from 3 to 28 μm in width and 6–26 μm in length, and large-scale thermal fractures ranged from 42 to 69 μm in width and 29–40 μm in length. The width and length of the optimal thermal fracture is 68.4 μm and 32.6 μm. With the decrease of thermal shock rate, the number of thermal fractures decreases slowly, but its width and length increase significantly. Under the condition of 200 °C, the width and length of the optimal thermal fracture 100 s are about 2.3 times that of 25 s.

Fig. 6
figure 6

Thermal fractures characteristic parameters at different thermal shock rate

Full size image

According to Fig. 6d, 19, 14 and 11 thermal fractures with opening potential were generated in 25 s, 50 s and 100 s, respectively, accounting for 54.2%, 56% and 50%, respectively. In addition, although the decrease of thermal shock rate increases the size of thermal fractures as a whole, the proportion of thermal fractures with opening potential has not been significantly improved due to the decrease of the number of thermal fractures, which indicates that weakening the heat exchange efficiency of fracturing fluid cannot significantly improve the possibility of complex fracture formation in hydraulic fracturing.

Effect of thermal shock times

According to the above results, the change of thermal shock rate has little contribution to the width and length of thermal fractures. Therefore, in this section, the shale temperature is set to 200 °C, the open hole section to 20 °C, the cooling time is 25 s, and the number of thermal shocks is 1,5 and 10 times, so as to clarify the influence of the number of thermal shock on the characteristic parameters of thermal fractures. The distribution of thermal fractures in borehole surrounding rock under different thermal shock times is shown in Fig. 7. According to Fig. 7, compared with a single thermal shock, when thermal shocks times is 5, the thermal fractures will not only widen but also extend forward obviously for a certain distance, and a large number of new small-scale thermal fractures will emerge. When the thermal shocks times is 10, only a small number of thermal fractures continue to extend forward, but new small-scale thermal cracks will continue to appear. At the same time, the thermal fracture will bifurcate during the extension process, which may be due to the local stress concentration at the tip of the thermal fracture, resulting in the secondary cracking of the fracture (Khalil and Emadi 2020). With the increase of the times of thermal shocks, the number of thermal fractures will increase, but the number of thermal fractures that continue to increase in width and length will begin to decrease. On the one hand, the cyclic thermal shock will continuously cause the cracking of the rock, so the formed thermal fractures have temporal characteristics (i.e., time-order character); on the other hand, the cyclic thermal shock will cause the accumulation of fatigue damage at the tip of the existing thermal fractures, so only the first formed thermal fractures have the most accumulated damage, and the trend of the increase in the size of their fractures is significantly better than the other accumulated thermal fractures with less damage. This difference in damage accumulation will increase with the continuous increase of the thermal shock times, thus showing a trend of decreasing the number of thermal fractures with increasing scale.

Fig. 7
figure 7

Distribution of thermal fractures in borehole surrounding rock at different thermal shock times

Full size image

The characteristic parameters of thermal fractures under different thermal shock times are counted, as shown in Fig. 8. According to Fig. 8, when the thermal shock times is 5, 44 thermal fractures are generated in the borehole surrounding rock; when the thermal shock times is 10, the number of thermal fractures continues to increase to 52. The higher the number of thermal shock, the number of thermal fractures shows a slow upward trend, which is mainly attributed to the generation of small thermal fractures under the action of cyclic thermal shock. The higher the number of thermal shock, the zoning characteristics of the thermal fracture scale become obvious. When the thermal shock times is 5, the small-scale thermal fractures are mainly distributed in the range of 2–24 μm in width and 5–28 μm in length, and the large-scale thermal fractures are mainly distributed in the range of 32–48 μm in width and 25–36 μm in length. The width and length of the optimal thermal fracture are 50.0 μm and 32.1 μm. When the thermal shock times is 10, the scale zoning characteristics are the most obvious. The small-scale thermal fractures are mainly distributed in the range of 2–35 μm width and 5–40 μm length. The large-scale thermal fractures are mainly distributed in the range of 58–75 μm width and 42–52 μm length. The width and length of the optimal thermal fracture are 73.7 μm and 50.4 μm. The optimal thermal fracture width and length of 10 thermal shock times are about 2.5 times that of single thermal shock. It is worth noting that the numerical simulation results of the thermal fracture width after thermal shock are consistent with the measured results by Guo et al. (2022), which indicates that it is reasonable and true to use the CZM to simulate the thermal cracking of borehole surrounding rock.

Fig. 8
figure 8

Thermal fractures characteristic parameters at different thermal shock times

Full size image

According to Fig. 8d, 19,7 and 5 thermal fractures with opening potential were generated after 1,5 and 10 thermal shocks, accounting for 54.2%, 15.9% and 9.6%, respectively. With the increase of thermal shock times, although the proportion of thermal fractures with opening potential decreases significantly, the width and length of these large-scale thermal fractures are significantly better than those of small-scale fractures, and the cumulative damage sum is also significantly greater than other thermal fractures. It can be inferred that if these large-scale thermal fractures can expand synchronously during hydraulic fracturing, the possibility of forming complex fractures will be greatly improved.

Hydraulic fracture propagation characteristics of shale after thermal shock

At present, hydraulic fracturing is generally divided into 5 stages: pipeline pressure test, pumping ahead liquid, pumping carrying fluid, pumping displacement fluid and pressure drop test. In the first two stages, the large amount of heat exchange process between low temperature fracturing fluid and high temperature rock mass can be regarded as a thermal shock. Based on the above concepts, this section reconstructs the damage zone of the borehole surrounding rock of, establishes a new numerical model of hydraulic fracture propagation and performs multi-factor calculation. The influencing factors include horizontal stress difference and pump displacement.

Effect of horizontal stress difference

The model parameters keep the thermal shock times unchanged (1), the maximum horizontal principal stress unchanged (85 MPa), and the pumping displacement unchanged (0.0025 m3/s). By changing the minimum horizontal principal stress (70–77–81 MPa), a numerical calculation model with horizontal stress difference of 15 MPa, 8 MPa and 4 MPa was established to analyze the influence of horizontal stress difference on the propagation morphology and complexity of hydraulic fractures.

When processing and analyzing the numerical calculation results, the propagation morphology of hydraulic fractures can be identified by the distribution cloud diagram of pore pressure. When the horizontal stress difference is 15 MPa, the hydraulic fracture morphology of different pumping time is shown in Fig. 9. According to Fig. 9a, in the initial stage of hydraulic fracturing, the hydraulic fracture starts from two thermal fractures, and the propagation direction is roughly parallel to the direction of the maximum horizontal principal stress. In the middle stage of hydraulic fracturing, one of the hydraulic fractures propagated for a certain distance and then stagnated (above Fig. 9b, defined as the secondary fracture), and the other hydraulic fracture propagated for a long distance (defined as the main fracture), and the propagation direction of both did not change significantly. At the end of hydraulic fracturing, the main fracture continues to propagate toward the sample boundary, and the secondary fracture propagation remains stagnant. It can be seen that under the condition of horizontal stress difference of 15 MPa, whether it is the main hydraulic fracture or the secondary fracture, the direction of its propagation is difficult to deflect. Secondly, due to the limited fracture net pressure, even if multiple thermal fractures fracture synchronously, only one main hydraulic fracture can be formed in the end.

Fig. 9
figure 9

The propagation morphology of hydraulic fracture when the horizontal stress difference is 15 MPa

Full size image

When the horizontal stress difference is 8 MPa, the hydraulic fracture morphology of different pumping time is shown in Fig. 10. According to Fig. 10, in the initial stage of hydraulic fracturing, the hydraulic fracture mainly starts from 4 thermal fractures. The main hydraulic fracture propagated for a short distance while the other three secondary fractures intersect and propagated stagnation. In the middle stage of hydraulic fracturing, the main hydraulic fracture was deflected, firstly expanding to the west for a certain distance (defining the direction of maximum horizontal principal stress as due north), then a second deflection occurred, expanding in the northwest direction, and the state of propagation of the secondary fracture did not change. At the end of hydraulic fracturing, the primary hydraulic fracture continued to propagate along the northwestern direction to the boundary without being deflected again. This suggests that under of horizontal ground difference of 8 MPa, although the main hydraulic fracture is prone to deflection in favor of the formation of curved and complex fractures, the secondary fracture is prone to intersection or stagnation of propagation near the well zone, which may lead to an increase in the flow friction within the fracture, and make the propagation distance of the hydraulic fracture limited.

Fig. 10
figure 10

The propagation morphology of hydraulic fracture when the horizontal stress difference is 8 MPa

Full size image

When the horizontal stress difference is 4 MPa, the hydraulic fracture morphology of different pumping time is shown in Fig. 11. According to Fig. 11, the propagation paths and overall morphology of hydraulic fractures are similar to that of the horizontal stress difference of 8 MPa, with the main difference that the deflection angle of the main hydraulic fracture is larger at the initiation of the fracture, and the secondary fracture also appears to intersect in the near-well zone or the propagation is stagnant.

Fig. 11
figure 11

The propagation morphology of hydraulic fracture when the horizontal stress difference is 4 MPa

Full size image

In order to further clarify the opening difficulty and morphology characteristics of the hydraulic fractures, the breakdown pressures, fracture lengths, maximum fracture widths of the main hydraulic fractures, and deflection angles of the main hydraulic fractures (the angle between the direction of fracture propagation and the maximum horizontal principal stresses) were counted for different horizontal stress differences, as shown in Fig. 12. According to Fig. 12, the breakdown pressure gradually increases (83.6–92.7–106.1 MPa) as the horizontal stress difference decreases. As the pre-existing thermal fractures continue to propagate to form the main hydraulic fracture, the fluid pressure only needs to overcome the horizontal minimum principal stress. The minimum horizontal principal stress gradually increases, and the breakdown pressure rises. The maximum width of the main hydraulic fracture gradually increases (0.08–0.10–0.11 mm), as the horizontal stress difference decreases. This suggests that two hydraulic fractures were formed under high horizontal stress difference condition, one of which did not propagate sufficiently but still required fluid flow distribution to maintain the openness, to some extent making the main fracture width smaller than that under low horizontal stress difference condition (only one hydraulic fracture). The variation trend of the fracture length is not obvious (61.76 mm, 11.07 mm for secondary fracture-81.88–67.75 mm), but it is not difficult to see that, under the condition of low horizontal stress difference, due to the change of the initiation direction and propagation direction of the main fracture, the main fracture propagation is affected by the horizontal maximum principal stress in the process of the propagation (deflected to maximum horizontal principal stress), thus increasing the distance of fracture propagation to a certain extent, which is most obvious when the horizontal stress difference is 8 MPa. It is worth noting that the ratio of the propagation distance of the secondary fracture to the propagation distance of the primary fracture (FLR) is 0.16 for high horizontal stress difference condition, which is relatively large compared to the medium and low horizontal stress difference, implying that the secondary fracture has the potential to continue to propagate. The deflection angle gradually increases (0.45°–42.28°–50.00°), as the horizontal stress difference decreases. Under the condition of high horizontal stress difference, the propagation direction of the main hydraulic fracture is largely controlled by the maximum horizontal principal stress, and the control effect is weakened by the reduction of the stress difference, the fracture is easy to be deflected, and the angle of deflection is then increased. In summary, the high stress difference is unfavorable to fracture deflection, which makes the main hydraulic fracture internal friction decrease, which creates conditions for the propagation of secondary fractures, but the propagation distance of secondary fractures is limited, and the overall fracture width is lower; the low stress difference is favorable to fracture deflection, and the main hydraulic fracture friction is larger, and its propagation pressure increases accordingly, which is unfavorable to the propagation of secondary fractures, but the overall fracture width is increased.

Fig. 12
figure 12

Characteristic parameters of hydraulic fractures under different stress differences

Full size image

Effect of pumping displacement

The model parameters keep the thermal shock times unchanged (1), the maximum (minimum) horizontal principal stress unchanged (85 MPa and 70 MPa). By changing the pumping displacement (0.005 m3/s and 0.01 m3/s), a new numerical calculation model was established to analyze the influence of pumping displacement on the propagation morphology and complexity of hydraulic fractures. The hydraulic fracture propagation morphology for a pumping displacement of 0.0025 m3/s as a reference group has been analyzed and illustrated in the previous section.

When the pumping displacement is 0.005 m3/s, the hydraulic fracture morphology of different pumping time is shown in Fig. 13. According to Fig. 13, in the initial stage of hydraulic fracturing, the hydraulic fractures were initiated by 2 thermal fractures, which is consistent with the reference group, and since it was under the condition of high stress difference, the direction of fracture propagation was all parallel to the direction of the maximum horizontal principal stress, and the two hydraulic fractures propagated at a comparable distance. In the middle stage of hydraulic fracturing, due to the increase of the displacement and the increase of fluid pressure, neither hydraulic fracture showed stagnation of propagation. The main hydraulic fracture propagated over a long distance, and the secondary fracture had a tendency to be deflected to the east. The secondary fracture is not fully extended because the direction of propagation is perpendicular to the direction of the maximum horizontal principal stress. At the end of hydraulic fracturing, the main hydraulic fracture appeared to be stagnant after continuing to propagate towards the boundary for some distance. But the secondary fracture first extended a distance westward, and then extended a long distance to the boundary in the direction parallel to the maximum horizontal principal stress. The fracture width and fracture length changed significantly compared with the low displacement.

Fig. 13
figure 13

The propagation morphology of hydraulic fracture when the pumping displacement is 0.005 m3/s

Full size image

When the pumping displacement is further increased to 0.01 m3/s, he hydraulic fracture morphology of different pumping time is shown in Fig. 14. According to Fig. 14, in the initial stage of hydraulic fracturing, the hydraulic fracture is still initiated and propagated by 2 thermal fractures parallel to the maximum horizontal principal stress. Different from the other two displacement conditions, the main hydraulic fracture has expanded a large distance compared with the secondary fracture. In the middle stage of hydraulic fracturing, the main hydraulic fracture continued to propagate and the width of the fracture increased significantly, while the secondary fracture appeared to be deflected to the west after propagating for a short distance, and the propagation rate was significantly lower compared with that of the main hydraulic fracture. At the end of hydraulic fracturing, the main hydraulic fracture propagation stagnates, the secondary fracture propagation rate becomes larger, and its path gradually returns to the direction parallel to the maximum horizontal principal stress while propagating to the boundary, and the length and width of the secondary fracture are significantly increased compared with the previous two displacement conditions.

Fig. 14
figure 14

The propagation morphology of hydraulic fracture when the pumping displacement is 0.01 m3/s

Full size image

Characteristic parameters of hydraulic fractures under different stress differences, are shown in Fig. 15. According to Fig. 15, the breakdown pressure gradually increased (83.6–84.8–96.9 MPa) with the increase of pumping displacement, which is consistent with the results of most scholars, mainly due to the fact that the higher the displacement, the faster the fluid flow rate in the wellbore, the higher the friction, the net pressure decreases, and the fluid pressure reflected to the wellhead increases at a certain rock breakdown pressure. The maximum width of the main hydraulic fracture increases gradually (0.08–0.14–0.20 mm). When cracking starts from the same thermal fracture, an increase in displacement means more fluid feed to the fracture, and if the hydraulic fracture is of a certain length, then the width of the fracture increases. A clear trend of increasing length of the fractures (67.56 mm, 11.07 mm for secondary fracture-81.30 mm, 85.51 mm for secondary fracture-98.80 mm,119.28 mm for secondary fracture) was observed. As the fracture propagation rate is accelerated by the increased fluid pressure, when the propagation of the main hydraulic fracture is stagnant, the secondary fracture will also propagate rapidly, and the fracture length of the secondary fracture may even exceed that of the main hydraulic fracture. This is evident from the FLR, which also increases from 0.12 to 1.05 eventually reaching 1.20 as the pumped displacement is gradually increased to four times that of the reference group. The variation of the deflection angle is less pronounced (0.45°–0.51°–0.49°), due to the fact that the direction of fracture propagation is tightly controlled by the maximum horizontal principal stresses at high levels of stress difference. In summary, high pumping displacement leads to an increase in breakdown pressure, which increases the construction difficulty. However, the high in-fracture net pressure provided by the increased displacement is, on the one hand, conducive to the increase of the main hydraulic fracture width and length; on the other hand, it enhances the propagation distance of the secondary fracture, which contributes to the formation of the complex fracture.

Fig. 15
figure 15

Characteristic parameters of hydraulic fractures under different pumping displacement

Full size image

Discussion

In view of the fact that the increase in pumping displacement has a significant increase in the characteristic parameters of both the main and secondary fractures, is it possible to maximize the extent of hydraulic fracture propagation if the pumping displacement is increased after the cyclic thermal shock. Based on the above concept, the numerical model with maximum and minimum horizontal principal stresses of 85 MPa and 70 MPa, pumping displacement of 0.01 m3/s, and 10 times thermal shocks was established and calculated, and the propagation morphology of its hydraulic fractures, as shown in Fig. 16. In the initial stage of hydraulic fracturing, all sizes of thermal fractures on the borehole perimeter are open, and the hydraulic fracture continues to be initiated by the two largest sizes of thermal fractures, similar to the previous numerical modeling results. In the middle stage of hydraulic fracturing, the previously defined main hydraulic fracture expands synchronously with the secondary fracture, and the secondary fracture expands at a slightly faster rate than the main hydraulic fracture, and the fracture widths of both increase significantly. At the end of hydraulic fracturing, both the main hydraulic fracture and the secondary fracture expanded steadily to the boundary, and there was no stagnation or deflection in the propagation. From the statistical results of fracture characterization parameters, the rupture pressure is 94.35 MPa, which is not significantly improved compared with the results of low-displacement and multiple thermal impacts. The main hydraulic fracture is 119.95 mm long and 0.18 mm wide, the length of the fracture is increased by 77.5% and the width of the fracture is increased by nearly 130% compared with a single thermal shock. The secondary length is 110.80 mm and the width is 0.08 mm, which is an increase of nearly 9 times in length and 22.3% in width compared to a single thermal shock. The FLR is 0.92, which is a significant improvement. In summary, the pumping mode of cyclic thermal shock followed by large displacement helps to further expand the degree of hydraulic fracture. For fracturing construction, the fracturing mode of cyclic thermal shock with increased displacement undoubtedly enhances the fracture complexity. When the ultimate pumping pressure is limited, only cyclic thermal shock is a good choice to increase the degree of fracture propagation. If it is difficult to implement the cyclic thermal shock.

Fig. 16
figure 16

The propagation morphology at large displacements and high times of thermal shock

Full size image

Conclusions

Based on the cohesive zone method, a coupled simulation method for sequential fracturing of deep shale is proposed. The thermal fracture initiation and hydraulic fracture propagation models of borehole surrounding rock are established, and the variation rule of thermal fractures characteristic parameters and the propagation behavior of hydraulic fractures are analyzed. The main conclusions are as follows:

  1. (1)

    During thermal shock, the number of thermal fractures increased by up to 96.5% with increasing shock temperature difference; the optimal thermal fracture width and length increased by up to 350%; and the percentage of thermal fractures with opening potential increased by up to 26.7%. Thermal fractures with more accumulated damage and opening potential are more likely to propagate under the driving force of fluid pressure.

  2. (2)

    Hydraulic fractures are always based on the initiation and propagation of thermal fractures. The magnitude of the horizontal stress difference is the main controlling factor of the hydraulic fracture propagation direction. The magnitude of pumping displacement is the main control factor of hydraulic fracture scale, with the increase of pumping displacement, the main hydraulic fracture width increases by up to 150%, the length increases by up to 46.3%, and the fracture length ratio also increases by nearly 7.5 times. The thermal shock times was the main controlling factor for the variation of secondary fracture length, and as the number of thermal shocks increased from 1 to 10, the secondary fracture length increased by up to 148.7%, and the FLR also increased by almost 3.4 times.

  3. (3)

    Increasing the thermal shock times can exacerbate the changes in the scale of thermal fractures, and increasing the pumping displacement can effectively promote the simultaneous propagation of multiple thermal fractures. The cyclic thermal shock and large displacement pump injection method are beneficial to improve the fracture formation effect of complex fractures in shale hydraulic fracturing and enhance the recovery of deep shale gas.