Fracture propagation laws of staged hydraulic fracture in fractured geothermal reservoir based on phase field model

Hydraulic fracturing is widely used in geothermal resource exploitation, and many natural fractures exist in hot dry rock reservoirs due to in-situ stress and faults. However, the influence of natural fractures on hydraulic fracture propagation is not considered in the current study. In this paper, based on the phase field model, a thermo-hydro-mechanical coupled hydraulic fracture propagation model was established to reveal the influence of injection time, fracturing method, injection flow rate, and natural fracture distribution on the fracture propagation mechanism. The results show that fracture complexity increases with an increase in injection time. The stress disturbance causes the fracture initiation pressure of the second cluster significantly higher than that of the first and third clusters. The zipper-type fracturing method can reduce the degree of stress disturbance and increase fracture complexity by 7.2% compared to simultaneous hydraulic fracturing. Both low and high injection flow rate lead to a decrease in fracture propagation time, which is not conducive to an increase in fracture complexity. An increase in the natural fracture angle leads to hydraulic fracture crossing natural fracture, but has a lesser effect on fracture complexity. In this paper, we analyzed the influence of different factors on initiation pressure and fracture complexity, providing valuable guidance for the exploitation of geothermal resources.


Introduction
With the rapid development of economic globalization, human demand for energy is increasing.The long-term use of conventional energy sources such as coal, oil and natural gas in large quantities has led to the emergence of environmental problems such as the greenhouse effect and land acidification, and it is urgent to find clean and renewable alternative energy sources.Geothermal energy, as a clean and renewable energy source, has good potential for development, and its development and exploitation scale are increasing in countries all over the world (Kepinska 2003;Trumpy et al. 2016;Vonsée et al. 2019).Geothermal resources are mainly divided into two major categories: hydrothermal type geothermal resources and hot dry rock type geothermal resources.Hydrothermal type geothermal resources are buried in 200-3000 m, mainly stored in high permeability pores or fractured media.Hot dry rock type geothermal resources are buried in 3000-10,000 m with a temperature higher than 180 °C.Currently, hydrothermal type geothermal resources that can be effectively developed account for only about 10% of the explored geothermal resources, and more geothermal energy is stored in the hot dry rock type geothermal resources (Xu et al. 2018).The research of hot dry rock started in the 1970s, and after decades of research and development, it has shown great application value and good development prospects (Bruhn 2002;Górecki et al. 2003;Ziabakhsh-Ganji et al. 2018;Leveni et al. 2019;Anderson and Rezaie 2019).Morton Smith from Los Alamos National Lab was the first to propose the idea and concept for the development of hot dry rock, and its development method is the Enhanced Geothermal System (EGS), the use of hydraulic fracturing and other methods to improve permeability in hot dry rock reservoirs, thereby increasing the volume of heat exchange between the fluid and the high-temperature reservoir.Hot dry rock reservoirs have natural fractures and high temperatures, how to adequately connect the natural fractures is the key to the efficient construction of EGS.
In recent years, scholars have studied the propagation mechanism of hydraulic fractures under natural fractures through experiments.It is generally believed that four methods after the hydraulic fracture meets the natural fracture, such as crossing, opening, arrest and offset (Zhang et al. 2021a).However, the defects of the current experimental research include: (1) the mechanical parameters of the experimental samples are different, which have impact on the experimental results, and it is difficult to get the ideal experimental law.(2) The experimental loading conditions cannot fully conform to the real formation requirements, resulting in the deviation of experimental results from the real formation.(3) The size difference between the experimental and the formation is huge, and the size effect is obvious, and the complex working conditions such as multiple wells and multiple fractures cannot be simulated due to the size limitation.
To better study the fracture propagation mechanism, numerical simulation methods are widely used.For fracture propagation problems, discrete models and continuous are widely used.The discrete models including the boundary element method (BEM) (Rezaei et al. 2019;Cong et al. 2021Cong et al. , 2022a)), the discrete element method (DEM) (Shimizu et al. 2011;Chamanzad et al. 2017;Cong et al. 2022b;Li et al. 2023) and the displacement discontinuity method (DDM) (Jiao et al. 2015).The discrete model can simulate the fracture initiation and propagation, but the matrix is assumed to be discontinuous, which is insufficient for the simulation of continuous media, and the fracture propagation results are affected by the grid size.The continuous model includes damage evolution law (Shojaei et al. 2014;Li et al. 2017Li et al. , 2019) ) and phase field methods (Miehe et al. 2010b;Mikelić et al. 2015;Miehe and Mauthe 2016).The damage evolution law that examines the stress condition of rock using failure criteria in relation to fracture initiation and propagation.The phase field model describes the fracture initiation and propagation by establishing a variational-based energy minimization framework, this is a better way than damage evolution law (Li et al. 2018) by solving the challenges of fracture tip discontinuity, or grid size sensitivity for fracture propagation, and hence the phase field model has been widely used for fracturing simulation in porous media.Bourdin et al (2012) used a phase field model for the first time to simulate hydraulic fracturing by assuming that the material of the model is a homogeneous, impermeable, and nonporous reservoir.Wheeler et al (2014) applied phase field model on porous media by the augmented Lagrangian approach and filled the gap for hydraulic fracturing simulation using the phase-field models.Xia et al (2017) and Zhou et al (2019) considered the impact of material heterogeneity on fracture propagation in 2D and 3D phase field model and realized hydraulic fracturing simulation.
Based on the phase field model, this paper establishes a fracture propagation model for hydraulic fracturing of hot dry rocks of natural fracture development.The influence of fracturing method, injection flow rate and natural fracture distribution on fracture propagation and morphology are analyzed respectively.The model adopts the simulation technique that realizes the study of fracture propagation and meets natural fractures in the hydraulic fracturing process of EGS and explores the fracture propagation mechanism of hydraulic fracturing.The research results can provide a theoretical basis for the optimal design of EGS.

Energy functional
According to Zhou et al. (2019) the total potential energy ψ of porous media rocks consists of kinetic energy ψ k , elastic energy ψ ε , fracture energy, external energy, fluid pressure, traction f and body force b.The total potential energy ψ is given by: To distinguish fracture and rock matrix, defining ϕ(x, t) ∈ [0,1], for fracture domain ϕ = 1, for matrix domain ϕ = 0. Length scale parameter l 0 is used to control ϕ change from 0 to 1, and the outer boundary is ∂Ω, Ω is an arbitrary bounded computational domain.as shown in Fig. 1.
The kinetic energy ψ k is: The elastic energy ψ ε is: (1) where λ and μ are Lame constants.

Phase field approximation
To divide the fractures and the rock matrix, ϕ(x, t) ∈ [0,1] is defined to divide the rock into fracture region (ϕ = 1) and matrix region (ϕ = 0), and the regularized fracture length l 0 is used to control the transition of ϕ from 0 to 1.The fracture surface density per unit volume of rock is (Miehe et al. 2010a): Fracture energy ψ k is: where, G c is the critical energy release rate, N/m.
The elastic energy ψ ε is decomposed into two components, tension and compression, with a strain tensor of the form ε = ε + + ε -, where ε + and ε -are the tensile and compressive tensors, respectively.
To avoid singularities in the elastic energy density, the model parameter k is defined and 0 < k << 1.Assuming that the phase field affects the stretching part of the elastic energy density, the elastic energy density is (Borden et al. 2012):

Governing equations for evolution of the phase field
By using Eqs.( 2), ( 5) and ( 6), the expression for the energy can be obtained as: Calculating the first-order variation of the energy functional L and setting its value to 0 to obtain: (3) where: üi =  2 u t 2 , σ ij are the stress components, which are calculated as: To prevent fractures already generated by the phase field model from reverting to an unopened state, the strain history field H needs to be introduced (Borden et al. 2012): + by H (x, t) in Eq. ( 10), the strong forms are rewritten as (Borden et al. 2012):

Governing equations for fluid flow
Define χ R and χ F as two linear functions, where the matrix region χ R = 1 and χ F = 0; the fracture region χ R = 0 and χ F = 1; the transition region expression is (Lee et al. 2016) (Fig. 2): where: c 1 and c 2 are the boundary of the transition region, the fracture region when ϕ ≤ c 1 , the matrix region when ϕ ≥ c 2 , and the transition region when c 1 < ϕ < c 2 .The fluid and solid properties of the transition domain are obtained by interpolating the matrix and fracture domains using the functions χ R and χ F .
In this paper, we assume that the fluid flow within a porous medium conforms to Darcy's law, and the fluid pressure control equation is (Yu et al. 2019): where, q m is fluid source term, kg/(m 3 / s); The Darcy velocity is expressed as (Xu et al. 2019): where where S is the storage coefficient, 1/Pa; compressibility coefficients c = c R χ R + c F χ F , c R and c F denote rock and fluid compressibility coefficients, respectively; φ R is the rock porosity; and K vol is the bulk modulus.
Therefore, the fluid pressure control equation is rewritten as:

Thermal stress calculation equations
The heat exchange between fluid and rock within the fracture generates thermal stress that effect fracture initiation and propagation.The thermal stress due to temperature changes are calculated (Zhang et al. 2021b).
where δ ij is the Kronecker symbol, uncaused; α T is the coefficient of linear thermal expansion, 1/K; T 0 is the initial temperature of the formation, K.
The equation for heat transfer between fluid and rock within a porous medium is (Li. 2018): where, c p,eff is equivalent isobaric specific heat capacity, J/ (m 3 K); λ eff is equivalent thermal conductivity, W/(m K); Q is heat source, W/m 3 ; c p,eff and λ eff are given by: where, c p,R and c p,F represent rock and fluid isobaric specific heat capacity, respectively, J/m 3 K; λ R and λ F represent rock and fluid thermal conductivity, W/(m K).

Validation
To verify the accuracy of the hydraulic fracturing model, we compare the simulation result with the experimental result carried out by Zhang et al (2017).The established physical model is shown in Fig. 3a, the model size is 0.2 m × 0.2 m, the stress of 10 MPa and 8 MPa are loaded at the top and right boundaries of the model, the left and bottom boundaries are supported by rolls, the injection fluid is water, the injection flow rate is 30 mL/s, the initial pressure is 0.1 MPa, and the specimen is heated at 333.15 K for 24 h before the experiment.Therefore, the fracture fluid and specimen temperature are assumed to be 333.15K.The radius of the circular hole is 0.0075 m.The specific simulation parameters are shown in Table 1.
The hydraulic fracturing experimental fracture morphology of Zhang et al (2017) and the fracture morphology of numerical simulation in this paper are shown in Fig. 3b and  c.Due to the stress difference is 2 MPa, the simulated results and the experimental results fractures all propagation along the direction of the maximum principal stress, and the fracture morphology and propagation paths are the same.By comparing with the experimental results, it is verified that the model in this paper has good accuracy.

Results and discussion
In this paper, a two-dimensional model is used to simulate the hydraulic fracture propagation problem during geothermal mining.The boundary loads of 2 MPa and 1 MPa are applied in the x-direction and y-direction, respectively, and other boundaries are fixed.The fracturing fluid is water, and the model has two wells with one injection and one production, and the fracturing method is zipper-type

Hydraulic fracture propagation morphology study
This section analyzes the fracture morphology during hydraulic fracturing.The model is fractured by zipper-type hydraulic fracture, and the fracturing sequence is the first stage of the first well, the first stage of the second well, the   second stage of the well, and the second stage of the second well, respectively.The injection flow rate is 8 kg/ (m 3 s), and the fracturing time of each stage is 4 s. Figure 6 shows the fracture propagation at different times.By comparing the fracture propagation at different times, with the increase of time, the first and second stages hydraulic fractures propagate and meet with the natural fractures, and then the hydraulic fractures open the natural fractures and propagate along the natural fracture direction.When the fracture propagates to the tip of the natural fracture, the hydraulic fracture deflects in the direction of the maximum principal stress due to the stress difference.When the third and fourth stages are fractured, the inter-stage and intercluster stress interference exist, resulting in difficulty in propagate the hydraulic fracture along the natural fracture direction when it meets the natural fracture, at which the hydraulic fracture cross the natural fracture.To compare the fracture complexity, the ratio of the number of fracture grids (ϕ > 0.95) to the total number of grids in the model is defined as the fracture area ratio according to Li et al (2022), and the fracture complexity can be reflected according to the fracture area ratio.Figures 7 and 8 reflect the fracture initiation pressure per cluster per stage and the fracture complexity, respectively.
From the fracture initiation pressure curve of Fig. 7, the fracture initiation pressure of second cluster is significantly higher than that of first and third clusters, which is caused by the inter-cluster interference during fracture propagation.With the increase of fracture stages, there is a small increase of fracture initiation pressure.Comparing the fracture complexity (Fig. 8), the fracture complexity gradually increases with the increase of fracturing time, but the growth of fracture complexity gradually decreases.Combined with the fracture morphology cloud diagram (Fig. 6), the fracture propagation is mainly concentrated in the first and second stages, and the fractures in the third and fourth stages are propagated much less than the first and second stages.In the second stage propagation, the fracture complexity increases the most significantly by 0.058, but in the third and fourth stages, the fracture complexity only increases 0.014 and 0.009, respectively.Analyzing the reasons for this, we can see that: (1) in the second stage of propagation, the pressure in the first stage still causes the fracture to continue to propagate; (2) the hydraulic fracture activates many natural fractures during the second stage propagation, and the above two reasons work together to increase the fracture complexity in the second stage.From third stage to fourth stage, from 8 to 16 s, the fracture complexity increases by only

Study on the fracture propagation of different fracturing methods
This section focuses on the effect of zipper-type hydraulic fracturing and simultaneous hydraulic fracturing on the fracture propagation morphology.The injection flow rate is 8 kg/(m 3 s), and the injection time is 4 s per stage.The schematic diagrams of the two fracturing methods are shown in Figs. 4 and 9, and the fracturing fracture morphology is shown in Fig. 10.The fracture propagation morphology of the two fracturing methods are basically the same, both mainly along the direction of maximum principal stress, and there is connection between stages and clusters, the difference is that the fracture propagation of zipper-type hydraulic fracturing is more adequate, the fracture length is longer, and the fracture propagation of the second and third stages of zipper-type hydraulic fracturing is higher than that of simultaneous hydraulic fracturing.
To compare the differences between the two fracturing methods, we investigate the differences between the fracturing methods and fracture complexity, as shown in Figs.11 and 12.To visually compare the differences in fracture initiation pressure, we take the average of the three clusters of initiation pressure as the fracture initiation pressure of the stage and compare the fracture complexity at the end of the second fracturing stage and the end of the fourth fracturing stage of zipper-type hydraulic fracturing with simultaneous hydraulic fracturing.By comparing the fracture initiation pressure (Fig. 11), the simultaneous hydraulic fracturing has a higher fracture initiation pressure than zipper-type hydraulic fracturing, which in turn affects the fracture complexity less than zipper-type hydraulic fracturing.The third and fourth stages of simultaneous hydraulic fracturing have significantly higher pressure than the first and second stages due to inter-stage stress interference.

Study on fracture propagation by injection flow rate
The injection flow rate is increased from 4 to 10 kg/(m 3 s), and the fracturing method is zipper-type hydraulic fracturing.To ensure the same injection volume under different injection flows, each injection time is set to t = 8, 5.3, 4, and 3.2 s, respectively.The fracture morphology after hydraulic fracturing is shown in Fig. 13.Similarly, we choose the average value of three clusters of fracture initiation pressure as the fracture initiation pressure of this stage.The comparison of fracture initiation pressure and fracture complexity under different injection flow rate conditions are shown in Figs. 14 and 15, respectively.The comparison of fracture morphology at different injection flow rates shows (Fig. 14) that the difference in fracture morphology and orientation is small.Comparing with the injection flow rate of 4 kg/(m 3 s), it is easier to open the natural fractures and connect with each other when the injection flow rate is 6 kg/(m 3 s) and 8 kg/(m 3 s).When the injection flow rate is 10 kg/(m 3 s), the first cluster of the second stage cross the natural fractures and connects with the first cluster of first stage, therefore, the increased injection flow rate of hydraulic fractures easily leads to the propagation of hydraulic fractures cross the natural fractures.
The first and second stages fracture initiation pressure increases with higher injection flow rate due to the higher injection flow rate in the fracture will result in higher frictional, therefore, higher pressure at the injection point required to reach the initiation pressure at the fracture tip and higher injection flow rates require higher pressures.However, the inter-stage interference between the third and fourth stages leads to higher fracture initiation pressure at injection flow rates of 4 kg/(m 3 s) and 6 kg/(m 3 s) initiation pressure than at injection flow rate of 8 kg/(m 3 s).It takes shorter time to reach the initiation pressure at injection flow rate of 8 kg/(m 3 s), which in turn leads to the highest fracture complexity (Fig. 15).More fracture complexity can be obtained with a medium injection flow rate, which is also consistent with Tan's experimental results (Tan et al. 2017).Although the injection time is longer at an injection flow rate of 4 kg/ (m 3 s), the longer time to reach the fracture initiation pressure results in a shorter hydraulic fracture propagation time, therefore, slower fracture propagation leading to a lower fracture complexity.As the injection flow rate increases, the hydraulic fracture propagation time is longer, and the fracture complexity is higher.Due to the shorter simulation time

Study of natural fracture distribution on fracture propagation
This section mainly analyzes the effect of natural fracture distribution pattern on fracture propagation.The natural fracture angles are 0° and 90°, respectively.The schematic diagram of natural fracture distribution pattern is shown in Figs.16 and 17.
By comparing the natural fracture random distribution, natural fracture angles are 0° and natural fracture angles are 90° fracture patterns (Fig. 18a and c), it can be seen that: the natural fracture distribution pattern has a great influence on the fracture propagation.When the natural fracture angles are 0°, a large number of natural fractures are opened by the hydraulic fractures and propagated along the natural fracture direction because the natural fracture direction is the same as the maximum principal stress direction, while when the natural fracture angles are 90°, the hydraulic fractures cross the natural fracture obviously, and almost all the fractures in the second and fourth stages cross the natural fractures and propagate along the maximum principal stress direction.And when the natural fracture is random distributed, the randomness of hydraulic fracture propagation is significantly higher than that when the angles are 0° and 90° (Fig. 18).
The fracture initiation pressure in the second and third stages shows natural fracture angle 0° > random distribution > natural fracture angle 90°, which is caused by the increased stress disturbance due to the horizontal propagation of the fracture.When the fourth stage initiation fracturing, the natural fracture angle 0° < random distribution < natural fracture angle 90°, and it can also be seen from the fracture morphology that the fracture propagation is higher when the natural fractures are 0° than when the natural fractures are 90° (Fig. 19).
The comparison of the fracture complexity shows (Fig. 20) that the fracture complexity is higher than the random distribution when the natural fractures are 0° and 90°.This is because both 0° and 90° open the natural fractures to a higher degree than random fractures, resulting in a slightly higher final fracture complexity than the random distribution.The fracture complexity is basically the same when the natural fractures are 0° and 90°.This is because when the natural fractures are 0°, the fracture mainly propagates horizontally, and the fracture length is longer, but the fracture width is smaller.However, when the natural fractures are 90°, the first and third stages of fractures propagate along the natural fracture direction due to the influence of the natural fracture, and the continue to propagate along the natural fracture direction with difficulty due to the stress difference, and the fractures are deflected to the direction of the maximum principal stress, resulting in smaller fracture lengths but increased fracture widths.The fracture complexity is basically the same when the natural fractures are 0° and 90°.

Conclusions
In this paper, based on the phase field model, a hydraulic fracture propagation model is constructed to study the influences of fracturing time, fracturing method, injection flow rate and natural fracture distribution on fracture morphology, fracture initiation pressure and fracture complexity, and the following conclusions are obtained: (1) With the increase of fracturing time, the fracture initiation pressure of each cluster increases, but due to the inter-cluster interference problem, the fracture initiation pressure of the second cluster of each stage is higher than that of the first and third clusters; the fracture complexity increases with the increase of time, but the growth rate gradually decreases.(2) Compared with simultaneous fracturing, the use of zipper-type fracturing can reduce the fracture initiation pressure and obtain a higher fracture complexity.The average fracture initiation pressure is reduced by 2.9% and the fracture complexity is increased by 7.2% for the zipper-type fracturing.

Fig. 1
Fig. 1 Schematic diagram of matrix and fracture of the phase field model

Fig. 2 a
Fig. 2 a Linear indicator functions χ R and χ F ; b The reservoir and fracture domains and K F represent rock and fracture permeability respectively, m 2 ; μ is the fluid viscosity, Pa s.The storage coefficient S is calculated by the formula(Zhou et al. 2019): hydraulic fracturing of staged fracturing, each well is divided into two stages, each stage has three clusters, and the cluster spacing is 12 m.The natural fractures in the reservoir are randomly distributed, and the natural fracture length ranges from 8 to 12 m.The model size is 120 m × 120 m.The grid number of the model is 59254.The specific model parameters are shown in Table 2 (Figs. 4, 5).

Fig. 3
Fig. 3 Comparison of fracture morphology of hydraulic fracturing: a Model size and shape; b Zhang et al. experimental result; c numerical simulation result

Fig. 6
Fig. 6 Fracture propagation at different fracturing time

Fig. 13
Fig. 13 Fracture morphology changes with injection flow rate

Fig. 16
Fig. 16 Natural fracture angles are 0°F (3) The fracture complexity increases first and then decreases, reaching the maximum at the injection flow rate of 8 kg/(m 3 s).Both low and high injection flow rates result in short hydraulic fracture propagation time and lead to low fracture complexity.(4) When the natural fractures are horizontal to the direction of the maximum principal stress, the hydraulic fractures are more likely to open the natural fractures and propagate along the natural fractures, and when the natural fractures are perpendicular to the direction of the maximum principal stress, the hydraulic fractures are easy to propagate cross the natural fractures.The natural fracture distribution has less influence on the fracture complexity.

Fig. 19
Fig.19The influence of natural fracture distribution on initiation pressure