Hydraulic fracturing is a complex seepage-damage-stress coupling problem. To accurately simulate hydraulic fracture propagation, it is necessary to take into account the pore seepage of the rock, the deformation of the fractures, the fluid flow in the fractures, and the effect of natural fractures. The numerical method based on continuous damage theory is an important means to simulate hydraulic fracture propagation. The method uses the internal variable, damage factor, to describe the damage degree of an element, and associates the stiffness, strength, and seepage parameters of the element with the damage factor, so that the hydraulic fracture propagation process can be equated to damage evolution of the model. The main modeling theories are briefly introduced in the following.

The permeability of the carbonate reservoir matrix is extremely low. However, to simulate the propagation of hydraulic fracture by using the continuous seepage-damage-stress coupling theory, the matrix still needs to be assumed to be a porous elastic medium. The deformation of the matrix by external loads and the process of pore fluid seepage satisfy the differential equations of stress equilibrium and the equation of conservation of seepage mass (Detoumay and Cheng, 1993; Xu et al., 2010; Shan et al., 2018), which are denoted as: ∇ σ ′ + α p I + b = 0 (1) ∇ k γ w ∇ p − b w + m T ε ˙ = 0 (2)

Where σ ′ is the effectively stress tensor, tension is taken as positive; α is the Biot coefficient; p is the pore pressure; b is the body force vector of rock;k is the reservoir permeability; γ w is element weight of water; b w is the body force of water phase; ε ˙ is the strain tensor, I is the Kronecker delta; m = 1 , 1 , 1 , 0 , 0 , 0 .

The process of hydraulic fracturing treatment involves continuously injecting fluid into the reservoir and holding up high pressure. When the pressure exceeds the strength limit of the rock, hydraulic fractures occur. Thus, a rock strength model for predicting hydraulic fracture initiation is needed, and here a hyperbolic Drucker-Prager plasticity model (Drucker and prager. 1952) is used as a criterion for fracture initiation, which can take into account the prediction of tensile and compressive shear failure of the rock (Hibbitt et al., 2016): F = l 2 + q 2 − p ⁡ tan ⁡ θ − d = 0 (3)

Where F is the yield function; q is the equivalent stress; p is the effectively mean stress; and θ are the friction angle and cohesion in p ∼ q space, related to the angle of internal friction; d is the intercept of the yield surface on the q-axis in the p ∼ q stress spaces. Figure 2 shows the comparison of the linear form and hyperbolic form of the Drucker-Prager criterion.

When the rock element reaches its strength limit under external loads, damage occurs. And macroscopic fracture occurs when the damage develops to a certain degree. The stiffness and strength of the rock at this stage can be expressed as a function of damage factor (Tang et al., 2002), For example, the stiffness and strength change with damage as: E = 1 − D E 0 σ f = 1 − D σ f 0 (4)

Where D is the damage factor; E 0 and σ f 0 are the modulus of elasticity and strength of the rock before damage; E and σ f are the modulus of elasticity and strength of the rock after damage.

The brittle behavior of rock is characterized by damage-displacement response, which could be used to define the energy required to open an element area of fracture (Hillerborg et al., 1976), as shown in Figure 3. When the rock element is deformed in an elastic phase, the damage factor is zero. It is generally accepted that damage occurs with plastic deformation and changes dramatically with increasing plastic deformation. The relationship between damage factor and plastic deformation can be fitted based on rock mechanical experimental parameters (Shan et al., 2018) D = 1 e − 1 / a − 1 e − u ¯ p m / a − 1 e − 1 / a − 1 (5)

Where u ¯ p m is the normalized equivalent plastic deformation, which is the ratio of the current plastic deformation to the ultimate plastic deformation; a is a material parameter that reflects the rate of evolution of the damage with plastic deformation.

When the rock is damaged under external load, the permeability of the rock is inevitably affected due to the formation of microfractures. Here, the permeability of damaged rocks is described in three stages based on the results of previous research (Zhou et al., 2010; Shan et al., 2018). For undamaged rock element, the permeability is a function of the mean stress; when the rock element is damaged but no macroscopic fractures occur, the permeability is considered to be mainly influenced by the mean stress and the degree of damage; after the formation of a macroscopic fractures, the permeability is defined using the cube law, which is expressed as: k = k 0 e − η σ m − α p D = 0 ξ D k 0 e − η σ m − α p 0 < D < 1 l l + b ⋅ ξ   u k 0 e − η σ m − α p + b l + b ⋅ b 2 12   D = 1   a n d   σ 3 > σ t 0 (6) where k 0 is the permeability of the undamaged rock; η is the coefficient of stress on rock permeability; σ m is the effectively mean stress; l is the characteristic length of the element; b is the width of hydraulic fracture, which is approximated here as the tensile plastic deformation of the element. ξ D denoted as the permeability damage multiplier as a function of the damage factor, which is assumed here to be linearly related to the damage variable D . ξ u is the permeability damage multiplier when the element is fully damaged.

The above are the main theoretical equations to realize the simulation of hydraulic fracture propagation. There is an obvious coupling relationship between seepage, stress, and damage. During the hydraulic fracturing process, as the fracturing fluid is injected, the rock near the wellbore is damaged under hydraulic pressure. With the damaged area increased, the rock permeability is also increased, leading more fracturing fluid infiltration into the rock to further reduce the effectively stress, which promotes the formation of macroscopic fracture under the tensile stress. This process repeatedly occurs at the tip of hydraulic fractures, which promotes the forward propagation of the hydraulic fracture. The specific coupling process can be seen in Figure 6. Due to the large number of physical quantities involved and the nonlinear deformation characteristics, the finite element method is used here to solve the above equations.

To improve simulation efficiency, commercial software, ABAQUS, is used for building the model geometry and solving formulas, and we develop a user subroutine code to associate the rock element parameters with the damage variable to achieve the coupling of stress, seepage and damage (Hibbitt et al., 2016). The model geometric is discretized by linear tetrahedral elements. The Formula (1–6) are solved by principle of minimum potential energy using implicit method. Before applying this simulation method to fracture-cavity reservoirs, the effectiveness of the simulation method should be verified first.

Due to the fact that the ability of continuous damage methods in simulating the free propagation of hydraulic fractures in three-dimensional space and interaction between natural fractures and hydraulic fractures has been proved (Shan et al., 2018 ; Shan et al., 2018), here we only demonstrate the ability of the continuous damage method to simulate the influence of cavities on the direction of hydraulic fracture propagation. We will refer to existing physical simulation results (Liu et al., 2019) and compare them with the numerical simulation results calculated by our model to verify the effectiveness of the numerical method proposed in this study. In Liu’s study, The experiments demonstrated the influence of cavities on the propagation path of hydraulic fractures under different conditions of geostress differences. Referring to his experimental process, we establish a geometric model of the same scale and assign the same material properties and boundary conditions. As shown in Figure 4 are the model geometry and the loading direction. It can been seen the cavity is pre-set on the potential expansion path of hydraulic fractures. Two loading schemes are designed to investigate the influence of geostress differences on fracture propagation paths refereeing Liu’s study, where in Case 1, the σ V / σ H ./ σ h are set to 16/14/7 MPa with a horizontal stress difference of 7 MPa; while in Case 2, the σ V / σ H ./ σ h are set to 16/14/9 MPa with a horizontal stress difference of 5 MPa. The two cases use the same injection scheme with a pump rate of 20 mL/min.

Figures 5–7 show the numerical results and their comparison with experimental results. The green area in numerical simulations is the rock damage area, and its dynamic evolution can be equivalent to the path of hydraulic fractures. The numerical simulations obtained results consistent with the physical simulations. That is under high horizontal stress difference condition, it is easy for hydraulic fractures to connect the cavity that located in the direction of maximum horizontal stress (Figures 5, 7A). However, when the horizontal stress difference is low, the hydraulic fractures are prone to bypass the cavity (Figure 6; Figure 7B) due to the stress concentration around the cavity, as expressed in Cheng and Kao’s work (Cheng et al., 2019; Kao et al., 2022). It can be seen the numerical simulation method perfectly reproduces the experimental results. Therefore, the new simulation method proposed is effective and can be used for simulating hydraulic fracture propagation in fracture-cavity reservoirs.

According to the typical heterogeneity characteristics of fracture-cavity carbonate reservoirs, it is not suitable to use plane assumption, therefore, a three-dimensional model is established for the simulation. The scale of the model is 100 m × 100 m × 60 m (Figure 8A), and a cluster of perforations is set up in the center of the model with perforation phase of 60° and perforation number of 6 (Figure 8B). Referring to on-site seismic logging data and literature parameters (Kao et al., 2022), as shown in Table 2, randomly distributed natural fractures and cavities are generated in pre-processing stage (Figures 8C, D), assuming the natural fractures are elliptical and the cavities are spherical. To ensure the accuracy of the simulation results and computational efficiency, the geometric model is meshed using a tetrahedral mesh, and the grid is encrypted at the locations of the wellbore, perforation holes, natural fractures, and cavities (Figure 9).

Table 3 shows the reservoir rock material parameters, including the parameters related to the damage-permeability evolution of the reservoir rock. For material parameters of natural fractures and cavities, the nonzero initial damage variables are set for the natural fractures and cavities so that the natural fractures and cavities have lower initial strength and stiffness and higher initial permeability compared with the rock matrix. Here the initial damage variable for the natural fracture elements is set to be 0.9 and the initial damage variable value for the natural cavity elements is set to be 1 (i.e., complete damage with very low strength and very high permeability).

It is assumed that the wellbore is along the direction of the minimum horizontal stress, and the reservoir is in the state of positive fault geostress. The formation pressure gradient is set to 1.1 MPa/100 m, the vertical geostress gradient is taken as 2.2 MPa/100 m, and the maximum horizontal geostress gradient and the minimum horizontal geostress gradient are taken as 2.07 MPa/100 m and 1.8 MPa/100 m, respectively. The constant pore pressure boundary condition is applied in our model considering the model geometry is large enough to avoid the boundary effect on hydraulic fracture propagation. The viscosity and pump rate of the fracturing fluid is set concerning the range of values commonly used in field design.

The above-mentioned is the basic theory and key parameters of the model, and the overall computational flow is shown in Figure 10. The geometric model is first established containing different sets, such as the rock body matrix, natural fractures, and cavities. Then different model sets are respectively assigned with different material parameters and initial damage variables. Then the whole model is meshed and the boundary conditions are applied before calculation. During the calculation stage, the elements’ damage state is updated in each analysis step to determine the elements’ mechanics and permeability parameters, which are used to obtain the new stress and deformation field of the model. As the coupling calculation continues, the dynamic propagation process of hydraulic fracture characterized by damage can be simulated.

Figure 11 shows the hydraulic fracture propagation process of a simulation case with reservoir burial depth of 3,500 m. Due to the fact that this study only considers the case of one perforation cluster, the selected pump rate value is relatively small. Here, fracturing fluid pump rate of 2 m³/min and viscosity of 3 mPa·s. It is worth noting that, although the legend in Figure 11 shows a damage range of 0∼1, the hydraulic fracturing is only marked in green. That is because we selectively expose one damage isosurface, which is damage of 0.5 marked in green, to display the three-dimensional morphology of hydraulic fractures more clearly. As we can see in Figure 11, at an injection time of 100 s, the hydraulic fracture starts from the injection hole location. When the injection time is 300 s, the hydraulic fractures initiated from each perforation converge to form a single fracture surface. With the continuous injection of fracturing fluid, the hydraulic fracture continued to propagate, connecting natural fractures and cavities.

Keeping the geostress gradient and pore pressure gradient unchanged, five burial depths are taken, which are 3,500 m, 4,500 m, 6,500 m, 8,500 m, and 10500 m. The pump rate of the fracturing fluid is set to 2 m³/min; the viscosity is 3 mPa·s; and the injection time is 2,000 s. Then the geostress parameters corresponding to different burial depths can be calculated by multiplying their respective gradients by depth and substituted into models for simulations, and the simulation results are compared afterward.

Figure 12 shows the influence of different burial depths on the propagation range of hydraulic fractures. To visually display the simulation results, part of the model surfaces are removed to observe the internal fracture morphology, and the hydraulic fractures are also extracted and displayed separately. It can be seen that under the influence of natural fractures and cavities, the hydraulic fractures show an obvious non-planar morphology and the propagation ranges decrease significantly with the increase of burial depth. The height of fracture propagation under each modeling condition is statistically presented in Figure 13, showing an approximate linear relationship with burial depth and decreasing by about 5% for every 1,000 m increased in reservoir burial depth.

To further examine the discrepancy of the influences of natural fractures and cavities on hydraulic fracture propagation under different burial depth conditions, a horizontal slice of each simulation result is displayed in Figure 14. By comparison, it is found that when the burial depth is relatively shallow (3,500 m, Figure 14A), the impact of geostress compression on fracture propagation is relatively small, and the hydraulic fracture expands sufficiently. The hydraulic fractures can both penetrate and deflect after meeting natural fractures, and then form a larger scale of complex fracture network. With the increase of burial depth (4,500 m–8500 m, Figures 14B–D), the restraint effect of geostress on hydraulic fracture propagation is enhanced, increasing the rock matrix breakage difficulty. The hydraulic fractures tend to deflect and propagate along natural fractures. As the burial depth continues to increase (10500 m, Figure 14E), on the one hand, the increase of geostress limits the propagation range of hydraulic fractures, and Natural fractures are compacted tightly enough to prevent the filtration of low viscosity fracturing fluid into them. Therefore, The trend of hydraulic fractures deflecting when encountering natural fractures will weaken, and the hydraulic fractures tends to penetrate the natural fractures, and the morphology of the hydraulic fractures tend to be simple and short, which is not favorable for connecting more cavities. When the burial depth is 3,500 m, 4,500 m, 6,500 m, 8,500 m, and 10500 m, the number of natural fractures and cavities connected by hydraulic fractures is 34, 30, 24, 17, and 15, respectively, as shown in Figure 15, showing a more obvious decreasing tendency with the increase of burial depth.

We summarize the impact of burial depth on hydraulic fracture propagation in Figure 16. It can be concluded that there are two types of cavities that can be relatively easily communicated by hydraulic fractures: cavities located on the main hydraulic fracture propagation path, and cavities closer to the edge of natural fractures have the potential to be connected by hydraulic fractures. These two types of cavities can be further divided into two categories based on their distance from the wellbore. Therefore, each block in Figure 16 is illustrated with four circles representing four categories of cavities. For reservoirs with shallow burial depth (3,500 m), hydraulic fractures have a wide range of expansion and complex morphology, which can communicate with these four types of cavities. As the burial depth increases (4,500 m–8,500 m), the increase in geostress have a negative impact on the expansion of hydraulic fractures. Hydraulic fractures are prone to deflecting rather than penetrating when encountering natural fractures. Therefore, due to the inability to communicate enough natural fractures, or the expansion path being altered by natural fractures, hydraulic fractures cannot communicate with cavities far from the wellbore. When the burial depth increases to 10500 m, hydraulic fractures tend to be short and flat, only able to communicate with hydraulic fractures and cavities closer to the wellbore. Therefore, for deeply buried fracture-cavity reservoirs, adjusting the viscosity of fracturing fluid is no longer meaningful. It is possible to consider increasing the amount of acid to expand the communication range of hydraulic fractures.

In this simulation case, the reservoir condition of 8,500 m burial depth is taken, and the pore pressure of 93.5 MPa, vertical geostress of 187 MPa, and maximum horizontal geostress of 176 MPa can be calculated correspondingly. Keeping the injection fracturing fluid pump rate of 2 m³/min and viscosity of 3 mPa·s and the injection time unchanged, the minimum horizontal geostress are taken as 152 MPa, 158 MPa, 164 MPa, 170 MPa, and 176 MPa, to investigate the influence of different geostress differences on hydraulic fracture propagation. The simulation results are shown in Figure 17.

Figure 17 shows the influence of different horizontal stress differences on the hydraulic fracture propagation pattern, and it can be found that when the horizontal stress difference is large, the hydraulic fracture surface stays perpendicular to the direction of the wellbore (Figures 17A, B), and as the horizontal tress difference decreases, the hydraulic fracture surface no longer stays perpendicular to the wellbore during propagation (Figures 17C,D). When the horizontal stress difference is 0, it produces a fracture surface that propagates along the direction of the wellbore in addition to a fracture surface that propagates perpendicular to the wellbore (Figure 17E). Also, horizontal slices of simulation results are displayed to investigate the mechanism of effect of different stress differences, as shown in Figure 18. As can be seen from the figures, at a horizontal stress difference of 24 MPa, hydraulic fracture breaks the rock more easily due to the smaller minimum horizontal stress, and hydraulic fracture can easily penetrate the fracture without being induced by the natural fracture to deflect (Figure 18A). As the minimum horizontal stress increases and the horizontal stress difference decreases, it becomes more difficult to break the rock, and deflection occurs easily when natural fractures are encountered (Figures 13A–C). When natural fractures exist in the vicinity of the wellbore injection hole, it will induce hydraulic fractures to propagate along the wellbore direction, especially when the horizontal stress difference is 0, a complete fracture surface that propagates along the wellbore direction will be develop (Figure 17E; Figure 18E), increased the complexity of hydraulic fracture morphology. It can be seen from the simulation results that the influence of stress differences and natural fractures on the propagation path of hydraulic fractures is similar to that in fractured reservoirs such as shale (Zhou et al., 2010 ; Hou et al., 2018 ; Liu et al., 2018). This indirectly proves the effectiveness of the new method in simulating the influence of natural fractures and stress differences on the propagation morphology of hydraulic fractures.

Figure 19 shows the influence of different horizontal stress differences on the hydraulic fracture connecting natural fractures and cavities, and it can be seen that under certain conditions of fracturing fluid injection pump rate and viscosity, high stress differences are beneficial for hydraulic fractures to connect more natural fractures and cavities. The reason for that maybe high stress difference can promote the hydraulic fracture break through the natural fracture barrier, and then connect more natural fractures and cavities, while under conditions of low stress differences, hydraulic fractures are susceptible to natural fractures, and hydraulic energy is captured by natural fractures and cavities near the wellbore, causing limited reservoir stimulation. Another reason for this result may be that under low stress differential conditions, hydraulic fractures tend to bypass cavities, as presented by previous research findings (Cheng et al., 2019; Liu et al., 2019; Kao et al., 2022). Therefore, for fracture-cavity reservoirs, high stress differences are beneficial for improving hydraulic fracturing treatment effects, unlike shale fracturing which prefers low stress differences to increase fracture complexity.

Create a new set of five models with the same geostress data. For each model, the reservoir burial depth is 8,500 m; the pore pressure is 93.5 MPa; the vertical stress is 187 MPa; the maximum horizontal stress is 176 MPa; the minimum horizontal stress is 164 MPa. To investigate the effect of pump rate on hydraulic fracture morphology, these five models take five different pump rates, 2, 4, 6, 8, and 10 m³/min. And the viscosity of the fracturing fluid is set as 3 mPa·s for all five models.

Figure 20 and Figure 21 show the hydraulic fracture propagation patterns under different pump rates. It can be seen that when the pump rate is increased from 2 m³/min to 4 m³/min, the propagation range of hydraulic fracture increases more obviously, but the hydraulic energy is still insufficient for hydraulic fractures to penetrate the natural fractures (Figure 21B). Thus the morphology of the hydraulic fracture does not change significantly (Figures 21A, B). When the pump rate is increased to 6 m³/min or more, hydraulic fractures not only propagate along natural fractures, but also penetrate them, forming a multi-branched fracture morphology (Figures 21C–E). Typical features of the influence of pump rate on fracture morphology can be illustrated as shown in Figure 22. Figure 23 summarizes the number of natural fractures and cavities connected by hydraulic fractures under different pump rates. As can be seen in, the number of hydraulic fractures connecting natural fractures and cavities increases significantly when the discharge rate is increased from 4 m³/min to 6 m³/min. Whereas the pump rate increases to greater than 6 m³/min, the increase in the number of hydraulic fractures connecting natural fractures and cavities is not obvious, which may be related to the scale of the model (Figure 23). In summary, under existing geological and engineering conditions set in this simulation case, in order to connect hydraulic fractures with as many natural fractures and cavities as possible, the fracturing fluid pump rate should be at least 6 m³/min per cluster.

In this simulation case, the burial depth is set as 8,500 m; pore pressure is 93.5 MPa; vertical stress is 187 MPa; maximum horizontal stress is 176 MPa; minimum horizontal stress is 164 MPa; the pump rate is 2 m³/min. And five different viscosities of fracturing fluid for five models are taken to exam the influence of fracturing fluid viscosity on hydraulic fracture morphology, they are 3, 30, 100, 150, and 200 mPa·s.

Figure 24 and Figure 25 show the influence of fracturing fluid viscosity on the hydraulic fracture propagation pattern. It can be found that when the viscosity of the fracturing fluid is very low (3 mPa·s), the fracturing fluid is heavily lost into the natural fractures when the hydraulic fractures meet the natural fractures, inducing the hydraulic fractures to deflect and propagate along the natural fractures (Figure 25A). As the viscosity of the fracturing fluid increases, the seepage coefficient of the fracturing fluid in the natural fracture becomes lower, and the tendency of the hydraulic fracture to deflect along the natural fracture is weakened, which is conducive to hydraulic fractures penetrating natural fractures. Overall as the viscosity of the fracturing fluid increases, the hydraulic fracture morphology becomes more and more flat and perpendicular to the wellbore (Figure 24). Especially when the viscosity of the fracturing fluid exceeds 150 mPa·s, the hydraulic fracture meets the natural fracture and penetrates directly, and the influence of the hydraulic fracture on the propagation path of the natural fracture is already very small (Figures 25D, E). Typical features of the influence of fracture fluid viscosity on hydraulic fracture morphology is summarized in Figure 26, form which, it can be seen that the hydraulic fracture morphology undergoes a transition from simple to complex, and then to simple with the change the fracturing fluid viscosity from low to high. Also, It is worth noting that through the statistics of the number of natural fractures and cavities connected by hydraulic fracture under different fracturing fluid viscosities, either too high or too low viscosity of the fracturing fluid is not conducive to the connection of more natural fractures and cavities by hydraulic fractures. But there exists an optimum viscosity interval, where the viscosity of fracturing fluid is not too high to ensure the expansion of hydraulic fractures along natural fractures, nor too low to ensure that hydraulic fractures can penetrate natural fractures. In this simulation case, the optimal viscosity of the fracturing fluid is 30 mPa·s. Under this viscosity condition, the hydraulic fractures can not only propagate along the natural fractures, but also penetrate of the natural fractures (Figure 26), which is beneficial to connect more natural fractures and the cavities (Figure 27).