The abundance of valuable geological deposits has led to the development of deep-earth resources has become the new normal. It has become an effective breakthrough in the rigid demand for energy (Zecheng et al. 2022; Zhang et al. 2022; Wan et al. 2024; Fuxi et al. 2024). It is predicted that China's deep and ultra-deep oil and gas production will exceed 40% by 2030(Yang et al. 2020; Nie et al. 2023), which is crucial to reach a balance between energy supply and demand. Deep reservoir fracturing refers to establishing an oil and gas route several kilometers below the surface. This process is considerably more challenging than shallow reservoir fracturing. At high temperatures, confining pressure, and geostress differences, the mechanical rock deformation and failure characteristics enhances plasticity, which changes the hydraulic fracture propagation behavior. Engineering design suitable for shallow reservoirs does not yield the desired outcome when applied to deep reservoirs (Guo et al. 2023a; Xiao et al. 2022; Qun et al. 2021; Ma et al. 2021; Haikuan et al. 2022; Jia Yunzhong, et al. 2021). Understanding the mechanism behind elastoplastic hydraulic fracture propagation in deep reservoirs is crucial for enhancing the efficacy of fracturing.

Significant recent development has been evident in CO₂ capture and storage (Liu et al. 2024; Benn et al. 2024; Hao et al. 2019), geothermal resource utilization (Guo et al. 2023b; He et al. 2023; Zeng and Zhang 2023), and deep oil and gas exploitation (Zhai et al. 2023; Zhonghua et al. 2020; Lingling et al. 2023), while the influence of in-situ conditions on rock deformation and failure behavior has attracted increasing research attention. At different temperatures (within 800 ℃), rock strength generally shows a nonlinear downward trend as the temperature increases (Yide et al. 2020; Zhang et al. 2016; Masri et al. 2014), with most degradation thresholds around 400℃. The rock strength in most reservoirs increases slightly in a temperature range within 200℃ (Zheng et al. 2023; Vishal et al. 2022a; Lubo et al. 2011). However, at higher confining pressures, the rock is mainly manifested as mechanical strengthening and brittleness weakening from the macro perspective. Based on observation methods such as CT and acoustic emission, it is mainly manifested as inhibition of fracture generation and expansion from the microscopic perspective (Lei et al. 2021; Guanghui et al. 2019; Wang et al. 2020; Qiu et al. 2023; Sheng-Qi et al. 2019; Lyu et al. 2021; Zhu et al. 2022). However, to better reflect the real-world scenario, mechanical research started to consider the combined effect of temperature and confining pressure (Li et al. 2021; Zhao et al. 2023; Ye et al. 2022). At present, it is generally believed that the linear mechanics theory applicable to shallow reservoirs cannot be applied to deep reservoirs. Establishing an elastoplastic constitutive model of deep reservoirs simulating the in-situ environment and effectively characterizing the deformation and failure behavior of deep rocks are essential to examine the deep reservoir fracturing mechanism.

Rock fracture behavior is closely related to the fracturing mechanism. Due to the plastic enhancement of rock in deep environment, fracture initiation and propagation have new characteristics. The competition between the growth rate of micro-fractures in the rock and the thermal expansion deformation of the rock leads to different change rules of fracture toughness (Zuo et al. 2018; Feng et al. 2018, 2017; Li et al. 2023): Fracture toughness has little change below 100℃, but increases in the temperature range between 10 and 200℃, and begins to decrease around 400℃. Scanning Electron Microscope (SEM) images show that obvious heat-induced micro-cracks are formed at 400℃, which is also the fundamental reason for the deterioration threshold of rocks (Vishal et al. 2022b). Due to the limitations of fracture experiments, relatively few studies have considered the effect of confining pressure on the fracture properties of rocks. The current consensus is that the fracture toughness of rocks under high confining pressure exceeds those under environmental pressure. When cracks undergo significant plastic deformation in a wide range of temperature and pressure conditions, the traditional fracture propagation theory that relies on linear elastic fracture mechanics is no longer valid. This necessitates the development of new criteria for fracture propagation based on elastoplastic fracture mechanics (Funatsu et al. 2014; Yang et al. 2021).

Current studies have extensively examined the law of hydraulic fracture propagation using the elastic reservoir hypothesis (Tan et al. 2017; Bing et al. 2018; Haddad and Sepehrnoori 2015; Zhang et al. 2019). The linear elastic constitutive model is used to deal with reservoir rock deformation, while the brittle fracture criterion based on the stress intensity factor is used to reasonably predict the hydraulic fracturing process of shallow brittle reservoirs. However, researchs found that for reservoirs with obvious plasticity, the fracture pressure is much higher than the predicted value of the linear elastic hypothesis (Wang 2016). Some important studies on elastic–plastic reservoir fracturing models have been carried out (Wang et al. 2016; Wenzheng et al. 2021), and the influence of rock mass plastic deformation has been effectively considered in the framework of fracturing mathematics. However, the linear cohesive fracture model is still adopted in the fracture propagation part, and the nonlinear fracture in the rock mass under the influence of plastic deformation is not taken into account. Therefore, it is necessary to consider the complex stress-seepage coupled elastoplastic deformation and nonlinear fracture propagation under the influence of plastic enhancement in deep reservoir hydraulic fracturing process.

From the perspective of deep reservoir plasticity enhancement, this paper establishes a model for elastoplastic hydraulic fracture propagation in deep reservoirs that considers elastoplastic rock deformation and fracture nonlinearity. It examines the hydraulic fracture propagation mechanism in deep reservoirs and explores the correlation between geological factors, engineering factors, plasticity degree and fracture propagation. The conditions promoting fracture propagation are analyzed according to the propagation mechanism, which provides theoretical support for deep reservoir fracturing design.

2.1 Research areas and experimental scheme

Due to the continuous process in extracting deep shale gas in the Sichuan Basin, significant amounts of industrial gas have been successfully collected from several deep shale gas wells in Luzhou, Yongchuan, and other areas. Therefore, this paper selected the shale from the Longmaxi Formation in the Luzhou area of the Sichuan Basin (Fig. 1).

An overview of the shale gas reserves in the Sichuan Basin and its surrounding areas (Fangzheng 2019)

The initial depositional center of the Longmaxi Formation in southern Sichuan is located in the Luzhou region, where deep-water shelf sedimentary conditions have led to the development of black carbonaceous and siliceous shales, which is highly conducive to shale gas formation (Zhankun et al. 2019). Deep shale is characterized by low porosity ranging from 2.5% to 4.4% and low permeability ranging from 0.006 to 0.01 mD. X-ray diffraction analysis shows that it is mainly composed of 33% and 56% quartz and 36% and 61% clay minerals, with a small amount of feldspar and carbonate minerals (Shanshan et al. 2019).

The experimental scheme shown in Table 1 was used to analyze the combined effect of temperature and confining pressure on rocks at different buried depths. Geological calculations showed that the temperature gradient below the surface was generally about 20–30 ℃/km, while every 1℃ temperature change produced a 0.45 MPa-0.55 MPa alteration (Jun et al. 2018). This empirical criterion was used to propose the environmental conditions at different burial depths for the experiments. It was assumed that the surface displayed a normal temperature of 25 ℃ and standard atmospheric pressure (0.1013 MPa), which corresponded to the uniaxial loading state without confining pressure conditions, which will not be repeated in the subsequent text.

2.2 Characterization of the rock deformation behavior in deep reservoirs

The in-situ rock compression failure experiment in combined temperature and pressure conditions revealed the reservoir deformation characteristics, providing a basis for establishing the constitutive rock relationship and numerical model of deep reservoirs. As shown in Fig. 2, the specimen was prepared in strict accordance with the standards recommended by the International Society of Rock Mechanics (ISRM). Considering the layered characteristics of shale, the sampling direction was uniformly perpendicular to the shale bedding. The raw rock was slowly cut into standard cylindrical specimens with diameters of 25 mm and heights of 50 mm using a high-precision emery line. The weights, sizes, and wave velocities of the processed specimens were measured, and abnormal specimens were removed. The test was performed using a GCTS-RTR-1500 high-temperature, high-pressure rock mechanics experimental system. After installing axial and circumferential strain gauges, the specimens were pushed into the pressure chamber and were slowly applied to the target confining pressure at a compression rate of about 2 MPa/min. The target temperature was gradually increased at a heating rate of about 2℃/min for 2 h to ensure uniform specimen heating. Finally, the axial loading rate was controlled by the strain rate at 0.02%/min until the specimen showed damage.

The design and principle of the in-situ compression failure experiment

The rocks showed compaction elastic deformation, plastic yield deformation, and post-peak failure during the compression failure experiment (Fig. 2). The initial pores and cracks were gradually compacted and closed during the early loading stage. Elastic deformation occurred in the entire rock system after compaction, and the curve was approximately linear. As the load continued to rise to the critical value, the curve began to deviate from the straight line. The stress concentration at the tip of the primary crack inside the rock prompted its evolution and plastic deformation, while the curve showed a downward trend. The microcracks continued to expand and penetrate to form macroscopic cracks. The bearing capacity decreased rapidly, resulting in post-peak failure.

Figure 3 shows the stress–strain curves in different experimental conditions. As the correlation between the temperature and confining pressure rose, the slope and peak values of the curve increased. An obvious plastic deformation stage appeared, showing a fundamental change in the plastic improvement characteristics. The volume strain of the rock changed from positive to negative, indicating volume expansion. On the basis of the main cracks in shale, some secondary cracks appeared. This was because the confining pressure had a restraining effect on the rock. The free thermal expansion of the shale was inhibited by higher temperatures, causing internal thermal stress and several microcracks after thermal fracture. The confining pressure inhibited the unstable propagation of the thermally induced microcracks, resulting in staggered slip, which enhanced the plasticity.

The deformation characteristics of the rock in different experimental conditions

The results showed that with the increase of burial depth, the elastic-brittle mechanical properties of the shale in the shallow reservoir environment gradually changed to elastoplastic mechanical characteristics in deep reservoirs. This paper used the Drucker-Prager yield criterion to establish a constitutive elastoplastic model for deep reservoir shale. The model considered the nonlinear hardening behavior of shale, specifically regarding its deformation properties under the combined effect of temperature and pressure. Additionally, the model accounts for the significant plastic flow observed during pre-peak nonlinear shale deformation.

The elastic stage of the rock before yielding is described by the linear elastic constitutive equation of generalized Hooke's law:

where \(\sigma_{ij}\) is the stress (MPa), \(\varepsilon_{kl}\) is the strain (dimensionless), and \(C_{ijkl}\) is the elastic stiffness tensor.

The plastic constitutive equation was used to describe the stress of the rock when it exceeded the yield point, while the Drucker-Prager yield criterion was used to evaluate the stress conditions required for the material to reach a plastic state:

where F < 0 indicates that the rock is in an elastic state, F = 0 indicates that the rock is in a plastic state, I₁ is the first stress tensor invariant, J₂ is the second deviatoric stress tensor invariant, and α and k are model parameters, which are related to rock cohesion c and the internal friction angle \(\phi\).

When the rock met the yield conditions, it entered the plastic deformation stage, and the stress–strain relationship showed a nonlinear trend. Flow criteria were introduced to characterize the nonlinear evolution process during the plastic stage, taking into account the assumption of the plastic strain direction. The plastic strain increment was determined using the non-associated flow rule, which was based on the constitutive framework of plastic mechanics and the first law of thermodynamics:

where \(G\) is the plastic flow potential and \({\text{d}}\lambda\) is the plastic factor. The direction and magnitude of the plastic strain increment were determined by the partial derivative of the plastic potential function and plastic factor, respectively. The mechanical behavior during the plastic deformation stage followed the consistency condition, with the stress state consistently located on the yield surface as the plasticity level changed. Consequently, the criteria for plastic loading and unloading were summarized as follows:

The plastic factor was determined according to the consistency condition. The geotechnical plastic mechanics showed that:

The incremental constitutive relation considering the elastic–plastic deformation characteristics was obtained as follows:

The equivalent plastic strain increment was obtained as follows:

The evolution process of the subsequent yield surface during plastic deformation was controlled by hardening criteria. Based on the isotropic hardening criterion, the equivalent plastic strain \(\varepsilon_{{{\text{eff}}}}^{{\text{p}}}\) was selected as the internal variable to reflect the change in the plastic loading history. A nonlinear isotropic hardening function k was constructed based on the first-order exponential decay function to control the hardening process and characterize the relationship between the subsequent yield stress and equivalent plastic strain. In this case, the Drucker-Prager loading condition was expressed as follows:

where A, B, and C represented hardening parameters. Hardening parameter A controlled the yield stress saturation value, the physical meaning of which was the critical yield strength during compression. Hardening parameter B controlled the initial yield behavior, the physical meaning of which was the initial yield strength during compression. Hardening parameter C controlled the saturation rate, the physical meaning of which was the rate at which the yield stress tended toward saturation during compression, that is, the ability to produce plastic flow.

The processed experimental data were used to obtain the variation laws of the yield and hardening parameters at different temperature and confining pressure combination and determine the constitutive model parameters. The relationship between the yield and hardening parameters and the burial depth was obtained by fitting the experimental data, as shown in Eq. (9):

where \(H\) is the reservoir burial depth (m), \(C_{\text{ini}}\) and \(C_{\text{fail}}\) are the initial and critical cohesion, (MPa), \(\phi_{\text{ini}}\) and \(\phi_{\text{fail}}\) are the initial and critical internal friction angles (°), and A, B, and C are the nonlinear hardening parameters (dimensionless).

The physical significance of the parameters in the nonlinear hardening function indicated that increased buried depth elevated the A and B nonlinear hardening parameters, which corresponded to a higher initial yield and peak strength of the shale exposed to combined high-temperature and high-pressure conditions. The decrease in the C nonlinear hardening parameter indicated that the shale approached the peak stress more slowly after the plastic yield, with a more obvious plastic flow ability. This analysis shows that the nonlinear hardening function established in this paper better reflected the physical phenomena during the experimental process in combined temperature and pressure conditions at different buried depths.

2.3 Fracture behavior characterization of deep reservoir rock

Rock fracture experiments are effective in revealing fracture initiation and expansion behavior. Since existing fracture mechanics experiments are typically performed in environmental conditions, it is difficult to transfer them to pressure vessels with high confining pressure for testing (Stoeckhert et al. 2016). Previous studies (Li et al. 2024) have conducted relevant fracture experiments and theoretical investigations involving deep reservoirs (Fig. 4). These studies designed and constructed a visual experimental system for rock fracture in combined in-situ temperature and pressure conditions. Digital image correlation technology is used as a monitoring method for fracture propagation. By calculating and analyzing the correlation and position changes of scattered spots in a series of images during the whole fracture process, the crack initiation and propagation behavior can be captured and collected, and the rock fracture process can be visualized. It provides a feasible method for studying rock fracture behavior in deep reservoir.

The design and principle of the in-situ fracture experiment (Li et al. 2024)

Since the fracture surface gradually opened and separated during the fracture process, the displacement change trend differed on either side (Fig. 5). The fracture opening displacement was determined by identifying the discontinuous zone of the surface displacement field during the fracturing. According to the cohesive fracture model proposed by Dugdale and Barenblatt (Barenblatt 1959), the mechanical fracturing behavior at the crack tip represented a closed cohesive fracture. The load–displacement curve and the displacement field determined via the DIC (Digital Image Correlation) method were used to quantitatively track the evolution of the fracture area during crack initiation and propagation, while the macro-global response was correlated with that of the local crack tip. For specific experimental analysis, see the research literature (Li et al. 2024). The plastic deformation of the rock mass was enhanced by higher combined temperature and pressure conditions, while the deformation openness increased and passivation occurred around the crack tip, making crack initiation more difficult. The curve of the relationship between cohesion and crack tip opening displacement tended toward a nonlinear convex trend. It gradually deviated from the linear constitutive cohesion relationship (black dotted line in the figure), while the fracture failure mode presented increased toughness. The relationship between fracture characteristic parameters and burial depth was obtained by fitting the experimental data (Li et al. 2024):

where H is reservoir depth (m), \(\delta_{\text{ini}}\) and \(\delta_{\text{fail}}\) are the critical and failure opening displacements of the cohesive cracks (mm), \(\sigma_{\text{c}}\) is the fracturing strength (MPa), and \(G_{\text{c}}\) is the fracture energy (N/m).

The rock fracture characteristics in different experimental conditions (Li et al. 2024)

The key to characterizing the fracture behavior of deep reservoirs using the cohesive fracture constitutive model is to establish the constitutive relationship between cohesive force and fracture surface opening displacement. The mechanism behind rock fracture crack formation essentially involves the progression from local damage accumulation at the crack tip to macro fracture. Therefore, based on the principle of damage mechanics, the relationship between cohesion and crack opening displacement was expressed in the form of damage factors. The traction and separation criterion proposed by Ren (Ren and Ru 2013) was used to characterize the constitutive nonlinear cohesive crack relationship:

where \(K\) represents elastic stiffness, \(T\) denotes the cohesive force with traction in the fracture plane (MPa) \(\varepsilon\) is the fracture surface deformation strain, and \(D\) represents the damage factor in scalar form. To describe the nonlinear damage evolution process of the rock after loading peak in deep reservoir conditions, the damage factors were expressed as:

where \(\delta_{0}\) and \(\delta_{{\text{f}}}\) represent the crack opening displacement at the beginning of cohesive crack damage and complete fracture (mm) and \(\alpha\) is the exponential coefficient, which is used to control the nonlinear degree of the damage evolution process. The relationship between the exponential coefficient and the burial depth was expressed as follows (Li et al. 2024):

where \(H\) is the buried depth of the reservoir (m) and \(\alpha\) is the exponential coefficient (dimensionless). Increased buried depth elevated the nonlinear degree of the post-peak softening stage, while the ductile fracture characteristics became more pronounced.

3.1 The hydraulic fracture propagation model

The physical hydraulic fracturing process involves the coupling of rock mass deformation, fracture dynamic expansion, and fluid seepage under hydraulic drive. It is assumed that the pores of the reservoir are completely saturated, the stress of the rock mass is shared by the solid skeleton and the pore fluid, and the rock skeleton is mainly controlled by the effective stress (Biot 1941):

where \(\overline{{{\varvec{\upsigma}}}}\) and \({{\varvec{\upsigma}}}\) represent the effective stress and total stress matrix (Pa), \(P_{\text{w}}\) denotes the pore pressure (Pa), and \({\mathbf{I}}\) is the identity matrix.

The rock mass stress balance equation was expressed by the virtual work principle. The virtual work of the rock mass at a certain time was equal to that generated by the force acting on the rock mass (Zienkiewicz and Taylor 2005):

where V is the integral space (m³), S is the integral space surface (m²), \({\mathbf{\delta \dot{\varepsilon }}}\) is the virtual strain rate matrix, \({\mathbf{\delta v}}\) is the virtual velocity matrix (m/s) \({\mathbf{f}}\) is the force vector per unit volume (N/m³), and \({\mathbf{t}}\) is the unit surface force vector (N/m²).

The hydraulic fracturing mathematical model used the constitutive equation of the rock skeleton to calculate the deformation caused by combined fluid pressure and ground stress. The elastoplastic constitutive relationship, which considered nonlinear strengthening, was incorporated into the governing equation of solid deformation to describe the deformation behavior of deep reservoirs:

where \({\mathbf{D}}^{{{\mathbf{ep}}}}\) is the elastoplastic stiffness matrix,\({\mathbf{D}}\) is the elastic stiffness matrix,\({\mathbf{D}}^{{\mathbf{p}}}\) is the plastic stiffness matrix, \(\text{d}\sigma_{ij}\) is the stress increment, and \(\text{d}\varepsilon_{ij}\) is the strain increment.

The fluid continuity equation was established according to the fluid mass conservation equation and Darcy's law. The fluid mass change rate over time in the rock microelement V was (Malvern 1969):

where \(\rho_{\text{w}}\) is the fluid density (kg/m³) and \(n_{\text{w}}\) is the porosity (dimensionless).

When the fluid mass entering the cell V per unit time was equal to the increase in its own mass, the fluid seepage continuity equation was obtained according to the Gaussian formula as follows (Malvern 1969):

where \(J\) is the volume change rate of the porous medium (dimensionless) and \(v_{\text{w}}\) is the flow velocity of the fluid (m/s).

The porosity of the saturated porous media was expressed as follows(Guangming et al. 2009):

where P_w is the pore fluid pressure (Pa), n⁰ is the initial porosity of the porous media (dimensionless), and k_g is the bulk modulus of the rock skeleton (Pa).

Assuming that the fluid flow in the rock pores satisfied Darcy's law, the relationship between the fluid velocity and pore pressure gradient was determined as follows (Guangming et al. 2009):

where \(k\) is the rock permeability (m/s) and g is the gravity acceleration (m/s²).

The transient nonlinear incremental equation of the reservoir stress-seepage coupling was obtained by combining the fluid continuity equation of the porous media with the stress balance equation of the rock skeleton in the finite element format. The finite element equation of the entire model system was acquired by assembling all the elements. The Newton–Raphson iterative method was used to determine the distribution of the stress, displacement, and pore pressure of the reservoir pore media nodes.

The nonlinear cohesive fracture constitutive relationship was used to describe the hydraulic fracture initiation and propagation process. The fracture propagation unit mainly consisted of three layers of nodes (Fig. 6). The upper and lower fracture plane nodes were used to simulate the fracture initiation and propagation process and were mainly controlled by damage factor D (Eqs. (11)–(13)). The intermediate nodes were used to simulate the flow behavior of the fluid injected into the fracture and drive fracture propagation. The tangential and normal of the fracturing fluid in the parallel and vertical fracture propagation directions indicated that fracture proliferation was generated by the movement of fluids within the fracture, although seepage caused some fluid loss to the surrounding rock matrix.

A schematic diagram of the hydraulic fracture propagation simulated using the nonlinear cohesive fracture model

Assuming that the fracturing fluid was an incompressible Newtonian fluid, the tangential flow behavior of the fluid in the fracture was controlled by the lubrication equation based on Poiseuille's law (Guo et al. 2015):

where \(q\) is the fluid flow rate per crack unit length (m³/s), \(\mu\) is the fluid viscosity in the fracture (Pa·s), \(w\) is the crack opening width (m), and \(\nabla p\) is the pressure gradient of the tangential fluid flow in the fracture (Pa/m).

Assuming that the fluid flow in the fracture conformed to the laminar flow pattern, the fluid pressure gradient in the fracture was expressed as follows (Guo et al. 2015):

where \(P_{0}\) is the fluid pressure when the fracturing fluid enters the fracture (Pa), \(P_{\text{i}}\) is the fluid pressure as the fracturing fluid crosses the fracture extension length (Pa), \(h_{\text{f}}\) is the crack thickness (m), and q is the flow rate per crack unit length (m³/s).

Assuming that a fluid exchange occurred between the permeable crack surface and the surrounding porous rock and adhered to Darcy's law of seepage, the normal fluid flow on the upper and lower surfaces of the unit was calculated as follows(Hagoort et al. 1980):

where \(q_{\text{t}}\) and \(q_{\text{b}}\) represent the fracturing fluid flow rate through the upper and lower surfaces of the fracture to the surrounding formation (m³/s), \(C_{\text{t}}\) and \(C_{\text{b}}\) represent the filtration coefficient of the fracturing fluid on the upper and lower surfaces of fractures (m/(Pa·s)), and P_i, P_t, and P_b represent the fluid pressure in the fracture and the pore pressure of the rock mass around the upper and lower surfaces (Pa).

The mass conservation continuity equation was used to determine the relationship between the flow state of the fluid in the hydraulic fractures at different times and the pumping capacity of the fracturing fluid, as follows(Peirce 2008):

where \(Q\left( t \right)\) is the pumping rate of fracturing fluid at any time (m³/s) and \(\delta \left( {x,y} \right)\) is the crack opening displacement at the x,y cell coordinates (m).

3.2 Model verification

The laboratory experiments showed that deep reservoir hydraulic fracturing must consider reservoir plasticity enhancement. The hydraulic fracturing model framework in this paper was subjected to subroutine secondary development to establish the deep reservoir hydraulic fracture propagation model. Figure 7 shows the overall calculation process.

The calculation of the deep reservoir hydraulic fracturing model

Finite element simulation software (ABAQUS)(ABAQUS 2016) was used to establish a hydraulic fracturing numerical model with a reservoir geometry size of 50 m × 50 m (Fig. 8). The planar strain four-node pore pressure elements were used for structured mesh division, the crack expansion area was locally encrypted by transition mesh, the global mesh size was 0.2 m × 0.2 m, and the transition coefficient was 10. A total of 16,250 elements were divided. A zero-thickness cohesive crack expansion element with pore pressure freedom was embedded along the expansion path perpendicular to the minimum principal stress direction. The injection point was located in the center of the model, while the initial fracture element was set around the injection point to represent the perforation hole. To simulate the confining effect of the surrounding rock in the reservoir, the normal displacement was constrained by the boundary around the model, and the initial effective stress and pore pressure fields were applied. The transient analysis method was used to calculate the fracturing process after balancing the ground stress. A fracturing fluid with a certain viscosity was introduced at the injection point according to the specified fracturing displacement.

The numerical hydraulic fracturing model

The hydraulic fracturing model in this paper was based on the plane strain hypothesis, which was consistent with that of the classic hydraulic fracturing KGD model. To verify the accuracy of the fluid–structure coupling effect and fracture propagation calculation in this model, the numerical solution was reduced to a fracture propagation problem in elastic conditions and compared with the analytical solution of the KGD model. Data from the shallow shale gas wells in the Longmaxi Formation in southern Sichuan were selected to validate the model. The borehole was drilled in the direction of the minimum horizontal in-situ stress. The maximum and minimum horizontal in-situ stress were 20 MPa and 15 MPa respectively, the elastic modulus was 15 GPa, the Poisson ratio was 0.25, the fracturing fluid viscosity was 1 mPa·s, the injection displacement was 6m³/min, and the tensile strength was 6 MPa. The hydraulic fracture width and half-length showed good agreement (Fig. 9).

A comparison between the numerical solution and the analytical solution

Since the cohesive crack unit describing the expansion behavior takes mesh as the basic unit, the volume in the crack will increase with each expansion unit, and the liquid cannot be replenalized in time when the crack tip breaks, resulting in a temporary decrease in the pressure in the crack. As the fracturing fluid is added to the newly formed fracture tip, the pressure will rise, gather enough energy and reach the initiation pressure before fracturing again, and the above process will be repeated, resulting in the fluctuation of the injection pressure curve, which is an important factor affecting the fracture expansion behavior. In order to analyze the influence of model Settings on the accuracy of hydraulic fracture propagation results, the influence of five mesh sizes on the injection pressure curves of intermediate fractures was analyzed, and the influence was reduced from the mesh optimization Settings to improve the simulation accuracy of fracture propagation. Under the condition of relatively large mesh size (0.4 m, 1 m, 2 m), the fluctuation of injection pressure curve is relatively severe, which directly affects the accuracy of the result, and even fails to converge, and the result of seam length decreases continuously. When the mesh size is reduced (0.2 m, 0.1 m), the fluctuation of hydraulic crack expansion pressure is reduced, and the result of crack length tends to be stable (red box in the Fig. 10), which is consistent with the analytical solution of the classical theory. However, too small mesh size will take a long time to calculate and waste computing resources. When the mesh size is 0.2 m, the injection pressure curve is less volatile, the accuracy is reliable and the calculation time is greatly reduced. So the model building process in this paper can ensure the calculation accuracy.

Optimization testing of model Settings

3.3 The elastoplastic fracture propagation mechanism in deep reservoirs

Based on the conditions of deep reservoirs at buried depths of 5000 m, this paper compared the hydraulic fracturing numerical model with brittle fracture in conventional elastic conditions and analyzed the differences between the simulation results of the two models. Additionally, it discussed the applicability of the two models for deep reservoir fracturing simulation and revealed the related hydraulic fracture propagation mechanism. Table 2 shows the input parameters of the specific models.

As shown in Fig. 11, when using the linear elastic reservoir fracturing model for calculation, the rock mass around the fracture met the linear elastic deformation, and no plastic strain was evident. However, when using the deep reservoir fracturing model considering the effect of plastic enhancement, the rock mass around the hydraulic fracture tip and the fracture propagation path showed obvious plastic strain, reducing of the fracture reconstruction effect.

A comparison between the fracture morphology of the elastic and elastoplastic reservoir conditions

When the hydraulic fracture extended to the same length, the width of the elastoplastic fracture increased significantly (Fig. 12). At the same injection volume, the crack length of the elastoplastic hydraulic fracture decreased by 33.3%, while the crack width increased by 44.2% due to the plastic deformation, compared with the elastic hydraulic fracture, resulting in short, wide fracture formation. The expansion of the elastoplastic fracture to the same scale required more injection energy, while the crack propagation pressure increased by 12.9%. Crack propagation was hampered by energy absorption during plastic deformation, causing pressure suppression in the crack. In addition, most of the current research models considering the influence of plastic deformation on hydraulic fracturing only consider the plastic deformation of the reservoir, and modify the linear elastic model to the constitutive model based on Mohr—Coulomb, ‌Drucker-Prager or damage theory. However, the influence of plastic deformation on fracture criterion of crack tip has not been taken into account, and it is proved by experiments in this paper that plastic deformation has a significant influence on fracture behavior. The model established in this paper considers both the characteristics of rock mass plastic deformation and nonlinear fracture, and is compared with other elastic–plastic models (Figs. 12b, c, d). It is found that the crack width increased by 13%, the crack length decreased by 19.6%, and the crack propagation pressure increased by 5.9% when nonlinear fracture is considered under the influence of plasticity. The model established in this paper has certain reliability.

A comparison between elastic and elastoplastic hydraulic fracture growth

Figure 13 shows the stress and strain field distribution during hydraulic fracture propagation in elastic and elastoplastic conditions, respectively, illustrating the mechanism behind the significant enhancement effect of plastic deformation on hydraulic fracture propagation in deep reservoirs. The reservoir rock was in the compressive state of the elastic deformation stage due to the initial ground stress environment. Local tensile stress was generated in the extruded rock mass when the fracturing fluid flowed into the perforation hole, concentrating the stress near the fracture tip. The compressive stress decreased in the direction perpendicular to fracture expansion. The difference between the maximum and minimum principal stress increased, allowing easier plastic zone formation near the fracture tip. The strain near the elastoplastic hydraulic fracture increased noticeably, while the fracturing fluid injection energy was primarily used for the plastic deformation of the surrounding rock mass near the crack tip. Consequently, the expansion volume decreased when the injection energy remained constant. According to the fracture experiment conducted in this paper (Fig. 5), when the elastoplastic fracture occurs in deep reservoirs, the fracture constitutive curve tends to be "convex" and deviates from the linear cohesion constitutive relationship (black dashed line in the figure), and the excess part is the energy dissipated by the plastic deformation. Calculated from the area of the fracture constitutive curve, the fracture energy increases by 18.9%. In addition, by comparing the stress field at the fracture front under the two reservoir conditions, it can be seen that the compressive stress of the rock mass at the hydraulic fracture front is larger under the elastoplastic condition, and the higher normal stress of the fracture surface increases the tensile strength and shear strength, that is, the expansion resistance increases, and the expansion of the hydraulic fracture is more difficult. The difference in the compressive stress at the fracture front under different reservoir conditions is caused by the difference in stress concentration.

The effective stress and strain fields around the elastic and elastoplastic hydraulic fractures

The induced stress curves in the fracture propagation direction and the wellbore direction (Fig. 14) indicated that plastic deformation softened the rock mass around the fracture tip, which weakened the stress concentration and reduced the tensile stress level in the expansion direction by 25%. This provided an effective extended shielding effect, which was also the primary cause of unsatisfactory fracturing in deep reservoirs. The crack tip passivation and the pressure inside the crack led to short, wide fracture formation. Consequently, the compressive stress gradient in the wellbore direction around the fracture increased, leading to higher interference stress between the expansion of multiple fracture clusters.

A comparison between stress induced by elastic and elastoplastic hydraulic fractures

Rock deformation and fracture behavior display significant changes in deep reservoir environments. The description by the linear elastic model exhibits some errors that may lead to an underestimation of fracture width and an overestimation of fracture length in the prediction results. Therefore, it is imperative to consider the utilization of the nonlinear hardening elastoplastic and fracture constitutive models proposed in this study. These models accurately depicted the behavior of elastoplastic fracture propagation in deep reservoirs, which involved plastic energy dissipation after the fracture tip was passivated, while the expansion resistance increased due to a decline in the stress concentration.

This paper investigated the impact of variations in the rock mass plasticity on the fracture propagation behavior in different geological and engineering conditions. It analyzed the relationship between geological and engineering factors, plasticity, and fracture propagation, as well as the methods that enhance the scale of fracture proliferation.

Deep reservoirs display complex structures and significant in-situ stress differences, which can reach 25 MPa in some areas. The different in-situ stress difference levels were analyzed (Fig. 15a). At an in-situ stress difference of 5 MPa, the maximum plastic strain of the reservoir around the fracture increased by 48.4%, the fracture expansion pressure rose by 3.8 MPa, the fracture length decreased by 21.4%, and the fracture width increased by 17.1%. As the horizontal in-situ stress difference increased, the stress state was closer to the yield envelope, and rock mass more easily entered the plastic state during fracturing. The internal friction angle significantly affected the heterogeneous rock strength due to the combination of sedimentary, diagenetic, and tectonic factors. Analysis was performed at different internal friction angles (Fig. 15b). As the internal friction angle increased, indicating a higher reservoir rock strength, the maximum equivalent plastic strain in the nearby reservoir region along the hydraulic fracture expansion path decreased. Compared with an internal friction angle of 15°, the maximum plastic strain decreased by 72% when the internal friction angle increased to 35°, while the fracture propagation pressure decreased by 9 MPa, the fracture length increased by 32.1%, and the fracture width decreased by 28.9%.

The changes in the fracture propagation due to geological and engineering factors

Fracturing fluid viscosity and injection displacement are important engineering parameters during deep reservoir fracturing, influencing fracture initiation and propagation. Therefore, based on the assumption that these geological factors remained constant, different fracturing fluid viscosities and displacements were selected for analysis (Figs. 15c, d). Compared with a 5 MPa·s fracturing fluid viscosity, a value of 45 MPa·s increased the maximum plastic strain by 4.6%, the expansion pressure by 1.6 MPa, the fracture length by 10.5%, and the fracture width by 8.3%. A higher fracturing fluid viscosity decreased the fluid loss to the surrounding reservoir matrix and increased the pressure gradient in the fracture. In general, the fracturing fluid viscosity minimally influenced the degree of plastic deformation but affected hydraulic fracture propagation. The plastic deformation degree of the rock was minimally affected by the injection displacement, even when the same amount of fracturing fluid was injected. When the injection displacement was increased from 6 to 18 m3/min, the maximum plastic strain increased by 5.2%, the expansion pressure increased by 3.4 MPa, the fracture length increased by 11.5%, and the width increased by 8.9%.

In summary, the effect of the fracturing fluid injection rate and viscosity on the plastic deformation of deep reservoirs is much lower than that of geological factors. The degree of rock plastic deformation is mainly affected by the horizontal ground stress difference and rock strength properties. Modifying the stress field between multiple fracture clusters and reducing the local stress difference is a feasible method to promote the expansion of elastoplastic hydraulic fractures in deep reservoirs.

It should be noted in particular that the model in this paper has improved the rock mechanics constitutive and fracture criteria within the original hydraulic fracturing mathematical framework, and can solve and analyze the fracture propagation problem under the characteristics of plastic enhancement in deep reservoirs. However, there are still some parts that can be further optimized, such as the instantaneous temperature difference effect of fracturing fluid injected under high temperature conditions in deep reservoirs. However, with the occurrence of reservoir seepage and heat exchange, this process is extremely complicated. In addition, from the field point of view, when the fracturing fluid is injected from the wellhead to the bottom hole, whether the temperature difference still exists is a problem that needs to be further studied. For the future research prospect, it is necessary to revise, supplement and optimize the model through rock mechanics research under the action of temperature difference and field investigation of actual temperature difference changes to achieve better results.

1. (1)

   The rock deformation and fracture behavior change significantly in deep reservoir, it is necessary to consider the nonlinear hardening elastoplastic constitutive model and nonlinear fracture model adopted in this paper, which can effectively reflect the fracture propagation behavior in deep reservoirs.

2. (2)

   The plastic deformation in deep reservoirs causes shorter, wider hydraulic fracture expansion. The propagation mechanism of hydraulic fractures is that the energy dissipates plastically after the fracture tip is passivated, and the propagation resistance increases due to the weakening of stress concentration effect.

3. (3)

   The degree of plastic rock deformation is mainly affected by in-situ stress differences. Reducing the local in-situ stress difference via reservoir reconstruction can decrease the inhibitory effect of plastic deformation on hydraulic fracture expansion, which can promote deep reservoir fracture proliferation.