Hydraulic fracturing is an efficient technology to extract hydrocarbon within natural caves. However, these caves can markedly affect the fracture propagation behavior. This paper proposes a novel hydraulic fracturing model to simulate the fracture propagation in poroelastic media containing the natural cave, utilizing the strengths of the phasefield method. By coupling the Reynolds flow with cubic law in fracture domain, free flow in cave domain, and lowpermeability Darcy flow in reservoir domain, the fracturecavereservoir flow governing equations are established. The Biot poroelasticity theory and fracture width are the links of hydromechanical coupling. The smooth phasefield is introduced to diffuse not only the sharp fracture but also the sharp cave edge. The fully coupling model is solved by a staggered scheme, which independently solves the pressure field and displacement field in inner cycle, and then independently solves the phase field in outer cycle. The proposed model is verified by comparing with the Khristianovic–Geertsma–de Klerk (KGD) model and Cheng's hydraulic fracturing model. Then, the interaction between hydraulic fracture and natural cave is investigated through several twodimensional and threedimensional cases. The result shows that the cave effect can make the hydraulic fracture deflect and raise its propagation velocity. Increasing the fracturecave distance, injection rate, and in situ stress difference can all decline the cave effect. The displayed cases also substantiate the capability and efficiency of the proposed model.
I. INTRODUCTION
Saturated geomaterials, such as carbonate rocks, often contain the natural caves, which results in multiscale pore structures. The geometry of pore structures can influence the mechanical properties of the geomaterials and also cause complex fracture growth path.^{1–3} For example, fracturecaved carbonate reservoirs are a common kind of carbonate reservoirs.^{4,5} They are composed of fractures, caves, and tight matrix, which leads to strong heterogeneity.^{6,7} Oil/gas is primarily stored within the natural caves. Hydraulic fracturing is an efficient technology to stimulate productivity of fracturecaved carbonate reservoirs, because it can connect the wellbore and natural caves by fluiddriven fracture.^{8,9}
Hydraulic fracture (HF) propagation is an extremely complicated process, including branching, deflection, merging, and so on. The existence of natural cave can cause strong stress disturbance, which will play a noticeable effect on the HF propagation behavior.^{10–12} However, the interaction mechanism between HF and natural cave is still vague. Moreover, according to reports in recent years, the hydraulic fracturing may create unfavorable transport channel accidentally. This circumstance can increase the possibility of groundwater being contaminated.^{13,14} Therefore, accurately predicting the HF propagation behavior in the existence of natural cave is essential and challenging, which contributes to the economic benefits and environmental protection. The primary objective of this paper is to develop a reliable numerical model to address this challenge.
Even now, various numerical techniques are developed to simulate the fracture initiation and propagation. In general, these techniques can be divided into discrete and continuous approaches. In the discrete approaches, the displacement discontinuity is directly or indirectly introduced into the computation. The extended finite element method (XFEM) is a widely used discrete approach.^{15–18} It adds the enrichment functions to depict the displacement discontinuity of fracture. In addition, the popular discrete approaches include cohesive element method,^{19} discrete element method (DEM),^{20,21} displacement discontinuity method (DDM),^{22,23} and so on. In the continuous approaches, the intact region and fractured region are viewed as a whole, and the crack is smeared at a certain area.^{24,25} In recent years, the phasefield method (PFM) has become the most common continuous approach. PFM originated in the late 1990s. Francfort and Marigo^{26} developed a variational model for the quasistatic cracks evolution. After that, Bourdin et al.^{27,28} established the early phasefield model based on the variational framework. Then, Miehe et al.^{29} established a solid–fluid coupling hydraulic fracturing model in porous media by PFM. So far, PFM has been utilized by many academic groups to establish hydraulic fracturing model.^{30–34} The highlight of PFM is that it applies a continuous scalar field to deal with the sharp fracture. Hence, the complex fracture evolution is governed by the scalar field. This treatment can eliminate the laborious tracing procedure of the fracture topological structure, which helps reduce the complexity of numerical implementation. PFM is constructed on the thermodynamic framework, and it is easily expanded to the multiphysical field problem.^{35–37} The strengths of PFM include the following two aspects: first, it can naturally detect and simulate the complex fracture evolution behavior, such as deflection, branching, and coalescence. The additional treatment for stress intensity factor is not needed. Second, PFM is conducted on the fixed topology, which can avoid remeshing.
In the past few years, some academic groups have studies HF propagation behavior with cave. Youn^{38} established a hydromechanical coupling model to study HF propagation in the existence of cave by XFEM. He used the levelset method to characterize the fracture and cave. The result indicated that HF will deflect and spread along the cave when it approaches the cave. After that, Cheng et al.^{39} improved Youn's model by adding a ramp function in blending nodes. They studied the influence of cave location and stress difference for the interaction mechanism. However, in Youn's and Cheng's models, the fluid flow in the cave is not considered. Luo et al.^{40} proposed that the cave is not empty and there is fluid or cement in it. They established a free flowseepagemechanical coupling model by XFEM. The simulation result revealed that the fluid pressure in cave can affect the HF propagation path. Liu et al.^{41} studied HF propagation with natural cave by PFM. They changed Young's modulus of cave to distinguish between cave and matrix. However, the matrix and cave are both impermeable in their model. Azarov et al.^{42} investigated the HF growth in a poroelastic media containing a cave. The result showed the stress arising near the cave can change the HF growth path. Yet, their model did not couple the fluid in the cave. Liu et al.^{43} simulated the HF evolution with natural fractures and natural cave by the element partition method. The cave is represented by an ellipse equation. Overall, the fluid flow in the cave is not considered in most studies. Therefore, it is necessary to couple the cave fluid in the hydraulic fracturing model. In addition, the HF propagation behavior is very complex near the cave, which results in great difficult to treat the stress intensity factor. If expand to the threedimensional issue, it becomes more complicated, and there are a few studies now. Hence, we employ the PFM to establish our hydraulic fracturing model, due to its strengths in dealing with the complex fracture propagation behavior.
In this study, we establish a novel hydraulic fracturing model for the fluiddriven fracture propagation in poroelastic media containing the natural cave. First, the fracturecavereservoir flow equations are derived. In fracture domain, the fluid flow obeys Reynolds flow with cubic law; in cave domain, it conforms to free flow, while in reservoir domain, it follows lowpermeability Darcy flow. Then, the PFM is used to diffuse the sharp cave edge and fracture. We reshape the total energy functional by adding the work done by fracturing fluid pressure and considering the pore pressure effect in elastic strain energy term. The coupled equations are discretized by the finite element method and solved by a staggered scheme. After that, the interaction mechanism between HF and natural cave is investigated in detail. The strong capability of the proposed model in dealing with 2D, 3D, and multiple cave issues is presented.
This paper is outlined as follows. The mathematical models and governing equations are presented in Sec. II. Then, the discrete process and numerical implementation scheme are shown in Sec. III. Next, the proposed model is verified in Sec. IV. After that, the interaction mechanism is studied by some numerical examples in Sec. V. The performance of the proposed model in multiple cave issues is displayed in Sec. VI. Finally, the concluding remarks of this work are drawn in Sec. VII.
II. MATHEMATICAL MODELS
For the mechanical model in Sec. II D, a specified displacement $ \u2009 u *$ is imposed on the Dirichlet boundary $ \u2202 D \Omega $. Meanwhile, a specified traction force $ \u2009 \tau $ is imposed on the Neumann boundary $ \u2202 N \Omega $.
Before obtaining the governing equations, some assumptions are made as follows:

The poroelastic domain is linearly elastic and isotropic.

There is only singlephase fluid in the poroelastic domain.

The fluid is viscous and slightly compressible.
A. The reservoir flow model
In this paper, we use the finite element method (FEM) to discretize the governing equations. Hence, Eq. (7) should be changed into weak formulation.
The notation $ \u301a \u301b$ denotes the “discontinuity” of the vector $ \u2207 p$ at the fracture face.
B. The fracture flow model
The second term on the right side of above formula only appears when fracture is connected with cave. Where $ \Gamma c$ is the boundary between fracture and cave, ds is the line integral.
C. The cave flow model
Caves form through a longterm geological movement, such as sedimentation, tectogenesis, and dissolution.^{52,53} Scholars classify caves into two categories: one is unfilled cave, and the other is filled cave.^{45} In this paper, we only consider the unfilled cave. Hence, there are no any loose substances inside the cave, and the fluid flow conforms to free flow. The Stokes equation is a very good appropriate choice to describe the fluid flow in unfilled cave. Many scholars combined Darcy–Stokes equations to establish flow model in porous domain and cave domain.^{54–58} Their simulation results show the pressure spreads very quickly in cave, and there is almost no pressure drop inside cave. Hence, some scholars proposed the equipotential body model in cave domain. This model assumes the pressure inside the cave is equal everywhere. Many literatures have proven the correctness of this model.^{59–65} In this paper, we choose the equipotential body model to describe the fluid flow in cave domain.
Due to the hydraulic communication, the fluid pressure is continuous. Therefore, a single pressure variable $ \u2009 p$ can represent $ p r$ and $ p f$. Meanwhile, the leakoff flux from fracture to reservoir is eliminated. So, there is no leakoff term in Eq. (23).
In Sec. II D, we will establish the mechanical model.
D. Energy functional
In the above formula, $D$ is the elastic matrix, and $\epsilon $ is the strain.
After establishing the total energy functional $F$, we will use the phasefield method (PFM) to regularize it.
E. Regularization by phasefield method
The phasefield representation of fracture and cave is shown in Fig. 2. The sharp fracture and cave edge are diffused by phase field.
Next, the work done by the fracturing fluid pressure needs to be regularized.
F. The decomposition of elastic strain energy
The crack is expected to open in the tensile regions and prohibited from interpenetration in compressed regions. Therefore, Amor et al.^{70} decomposed the trace of strain tensor into positive and negative part: $ t r + ( \epsilon ) = max ( t r ( \epsilon ) , 0 )$, $ t r \u2212 ( \epsilon ) = max ( \u2212 t r ( \epsilon ) , 0 )$.
G. Minimization of energy functional
III. NUMERICAL IMPLEMENTATION
This section will delineate the numerical implement of the proposed multifield coupling model. First, the FEM is used to discretize the spatial domain and the finite difference method is used to discretize the time domain. After that, a staggered scheme is designed to solve the discretized equations. The details are as follows.
A. Governing equation discretization
A structured eightnode hexahedral element is used to discretize the spatial domain in our work. Hence, the numerical algorithm can simulate the fracture propagation in threedimensional space. If we assume the plane strain conditions, the numerical algorithm is also suitable for simulating the fracture propagation in twodimensional space. The code of this numerical algorithm is completed by C language.
Equation (56) involves the fracture width $ \u2009 w$ and the volumetric strain $ \u2009 \epsilon vol .$ The fracture width is the difference of displacement on both sides of the fracture face: $ w = \u301a u \u301b \u22c5 n \Gamma $. Since the fracture width and the volumetric strain are the function of displacement, Eq. (56) contains the pressure field and the displacement field.
The fixed stress split method^{72} is employed to make, Eq. (56), decoupling from deformation. It is an unconditional stable algorithm, which is appropriate for the multifield coupling problem. Specifically, the fracture width and the volumetric stress are assessed using the value of previous iteration step when solves Eq. (56).
Section III B provides a detailed explanation of the implementation scheme.
Because the location of the fracture is variable during hydraulic fracturing, it is difficult to directly conduct numerical integration on the fracture domain. Hence, a regularized method is used to transform the integral in fracture domain $\Gamma $. Correspondingly, the integral in reservoir domain $ \Omega r$ is also transformed.
Please see Appendix for the discretizing process.
B. Numerical implementation scheme
In Sec. III A, we have obtained three discrete governing equations: Eqs. (53), (54), and (60) The process of solving the three equations is called module_u, module_v, and module_p, respectively. The three modules are interrelated closely. Module_u solves the displacement field $u$. Based on the solution of the displacement field, relevant parameters such as the stress, strain, and fracture width can be calculated. Module_u supplies the fracture width and volumetric strain to module_p. In addition, it supplies strain energy to module_v. Module_p is designed to solve the pressure field $p$. It supplies fluid pressure to module_u and module_v. Module_v is designed to solve the phase field $v$. It provides indicator function for module_p and phasefield variable for module_u. The Newton–Raphson method is applied in the three modules to solve the respective equations.
There are two ways to solve the proposed multifield coupling model, which are monolithic approach and staggered scheme. The monolithic approach is capable of simulating the real physical process. Yet, it needs more precise mathematical techniques to obtain the convergence, because of the nonconvexity of energy functional.^{73,74} The staggered scheme is robust, but it needs more costs of computing.^{75,76} In this paper, we choose the staggered scheme to solve the multifield coupling model.
The detail of the staggered scheme is shown in Fig. 3. First, we construct the initial $ v$ field for the cave domain and the perforated fracture domain. Next, the boundary conditions and initial conditions for the multifield coupling model are established. After that, the solving cycle of the staggered scheme begins. The staggered scheme involves three cycles: external cycle, middle cycle, and inner cycle. In the middle and inner cycles, the time is separated and the physical field iterates at each time step. Further explained, the displacement field $u$ and pressure field $p$ are solved by keeping the phasefield $v$ fixed in the inner cycle. Then, the phasefield $v$ is solved by keeping $u$ and $ \u2009 p$ fixed in the middle cycle.
In the inner cycle, the fixed stress split method is applied to solve $u$ and $p$. First, module_p is executed to obtain the fluid pressure $ \u2009 p$ through retaining the volumetric stress fixed. Then, module_u is executed to obtain the displacement $u$ by the updated $p$ from module_p. After that, it iterates between module_p and module_u until the error $ \epsilon inn$ is less than the predetermined tolerance $ \epsilon t$. In the middle cycle, the phasefield $v$ is solved using the convergent solutions of $ u$ and $p$ from the inner cycle. The irreversible condition is set here to ensure that the fracture does not recover to the unbroken state. If the error $ \epsilon mid$ is less than the tolerance $ \epsilon t$, the computation in new time step begins. In addition, the final convergent solutions of $ u , \u2009 p$, and $ \u2009 v$ in current time step are used as the initial values in next time step.
IV. MODEL VALIDATION
This section offers three numerical cases to verify the proposed model. The proposed model couples three systems including fracture, cave, and reservoir. However, as far as we know, there is no analytical solution of the hydraulic fracturing model when natural cave exists. Hence, the basic validation idea is as follows: first, we compare our model with the Khristianovic–Geertsma–de Klerk (KGD) model, which is a classic case for verifying the hydraulic fracturing. Through this comparison, the correctness can be verified when cave is not considered in our model. Next, in order to verify the correctness when cave is considered, we compare our model with the model of Cheng et al.^{39} through two numerical cases.
A. Comparison with KGD model
The KGD model was first established in the study of Khristianovic, Geertsma, and de Klerk.^{77,78} In the KGD model, a single fracture propagates under the constant fluid driven. It is a planestrain issue and the fracture height is invariant. The reservoir is assumed impermeable, so the fluid leakage from fracture is not existent. In addition, the reservoir is linear elastic, isotropic, and homogeneous. There are many analytical solutions and semianalytical solutions in the KGD model. In this paper, we choose the analytical solutions of Spence and Sharp^{79} to verify our model. The characteristic of Spence's model is injecting viscous fluid into an impermeable reservoir. Hence, their model is in viscositycontrolled regime.
The geometry and boundary conditions for simulating KGD issue are shown in Fig. 4. In the KGD model, a preexisting fracture propagates in an infinite domain. In our modeling, the computation domain is set as $ \u2009 200 \xd7 200 \u2009 m 2 ,$ and its center is placed an initial fracture, with the halflength equaling 2.5 m. The four boundaries are restricted for the displacement field. The above setting is to simulate the infinite domain approximately. An injection point is set at the center of fracture, with the rate $ \u2009 Q = 3 \xd7 10 \u2212 4 m 2 / s .$ The total injection time is 12 s. To simulate the impermeable reservoir, the reservoir permeability is set very small with the value of $ 2 \xd7 10 \u2212 22 m 2 .$ The initial pore pressure is $ \u2009 1 \xd7 10 \u2212 6 MPa .$ Other important parameters are as follows: critical energy release rate $ G c = 100 \u2009 N / m ,$ fluid viscosity $ \mu = 1 \xd7 10 \u2212 3 \u2009 Pa \u22c5 s ,$ Young's modulus $ E = 15$ GPa, Poisson's ratio $ v = 0.2 ,$ and time interval $ \u2009 \Delta t = 0.05 \u2009 s .$
We focus on the evolution of fluid pressure, fracture halflength, and fracture width vs time. Figure 5 shows the comparison results. We can find the fluid pressure at injection point rises nonlinearly prior to the fracture initiation in Fig. 5(a). Then, it gradually decreases with the fracture propagation. The analytical solution and numerical solution match well in terms of fluid pressure after the fracture initiation. The numerical solutions of fracture halflength and fracture width are also matched with the analytical solutions, which are shown in Figs. 5(b) and 5(c). Although these two matches are not as good as the fluid pressure match, these results are acceptable.
There are some reasons that cause the matching difference. First, the phasefield method makes several assumptions for regularizing the hydraulic fracturing model. Second, the computational domain cannot be infinite. When the fracture length becomes larger, the fixed boundary will influence the numerical result.
To sum up, the comparison results between our model and the KGD model are favorable. They verify the correctness of our model when cave is not considered.
B. Comparison with Cheng's model
In this verification example, the natural cave is considered in the reservoir. We compare the proposed model with the model of Cheng et al..^{39} They established a numerical model to simulate the fracture propagation when natural cave exists by XFEM. In their study, four numerical cases were set to study the impact of cave location. Here, we choose to compare with the first and second numerical cases to verify our model.
Figure 6(a) shows the geometric structure of the first numerical case in Cheng's model. The reservoir size is $ \u2009 50 \xd7 50 \u2009 m 2 .$ The right and bottom boundaries of reservoir are fixed. The maximum and minimum horizontal principal stresses are applied to the left boundary and top boundary, respectively. The reservoir center places a natural cave, with the radius of 3 m. In their model, the cave is empty and there is no fluid in the cave. An initial fracture is set at the left boundary, whose length is 1.8 m. The main input parameters in our model and Cheng's model are shown in Table I.
Parameters .  Our model .  Cheng' s model . 

Young's modulus E  35 GPa  35 GPa 
Poisson's ratio v  0.23  0.23 
Reservoir permeability k_{r}  2 mD  2 mD 
Fluid injection rate Q  0.05 m^{2}/s  0.05 m^{2}/s 
Fluid viscosity μ  10 $ \u2009 mPa \u22c5 s$  10 $ \u2009 mPa \u22c5 s$ 
Horizontal stress difference $ \u2009 \sigma H \u2212 \sigma h$  5 MPa  5 MPa 
Biot coefficient $ \u2009 \alpha $  1  / 
Porosity ϕ  0.04  0.04 
Critical energy release rate G_{c}  200 N/m  / 
Parameters .  Our model .  Cheng' s model . 

Young's modulus E  35 GPa  35 GPa 
Poisson's ratio v  0.23  0.23 
Reservoir permeability k_{r}  2 mD  2 mD 
Fluid injection rate Q  0.05 m^{2}/s  0.05 m^{2}/s 
Fluid viscosity μ  10 $ \u2009 mPa \u22c5 s$  10 $ \u2009 mPa \u22c5 s$ 
Horizontal stress difference $ \u2009 \sigma H \u2212 \sigma h$  5 MPa  5 MPa 
Biot coefficient $ \u2009 \alpha $  1  / 
Porosity ϕ  0.04  0.04 
Critical energy release rate G_{c}  200 N/m  / 
It is worth noting that Cheng et al. did not establish the fluid flow model in cave domain. Hence, there is no pressure inside cave, while in our model, the injection fluid can enter the cave and the cave pressure will change vs time. To better match their model, our idea is keeping the cave pressure fixed. Under this setting, the cave pressure remains constant and is always equal to the initial pressure. The similarity is that there is no pressure change inside cave domain of the two models. In this way, the influence of cave pressure change on solid deformation can be neglected.
Figure 6(b) presents the numerical results of our model and Cheng's model. The fracture propagation path in Cheng's model is depicted by the white dotted line, while the fracture propagation path in our model is depicted by the phasefield evolution. We can see that the fracture propagates horizontally to the right side and finally is connected with the cave. The result of the two models matches well.
The second numerical case is presented in Fig. 7. The only variation is that the cave location is changed. The cave moves up, and the cave angle, which can be seen in Fig. 7(a), is $ \u2009 10 \xb0$. The input parameters remain unchanged. Figure 7(b) presents the numerical results. The fracture propagation path in Cheng's model is depicted by the white dotted line, while the fracture propagation path in our model is depicted by the phasefield evolution. We can find that the fracture first propagates to the right side and then deflects up. Finally, the fracture is connected with cave. In a word, the matching result of the two models is good. Yet, there is a little difference of the propagation path between the two models. This is because our model is built upon PFM, while their model is upon XFEM. There may be a little calculation difference between the two methods.
To sum up, the correctness of our model is fully verified through the comparison with KGD model and Cheng's model.
C. Sensitivity analysis of the regularization parameter
The regularization parameter $\u03f5$ can influence the numerical simulation. Making a reasonable choice for it is crucial in our model. Therefore, we will discuss the reasonable selection of $\u03f5$ here.
In our model, the total energy functional [Eq. (29)] is regularized to Eq. (39) by introducing the diffusive interface to replace the sharp interface. Equation (39) converges to Eq. (29) as the regularization parameter $ \u03f5 \u2192 0$. Hence, the error between regularized model and actual model decreases as $\u03f5$ decreases. The ideal thought is to choose the regularization parameter $\u03f5$ as small as possible. However, previous studies have shown that $\u03f5$ is related to the grid element size (h_{w}) in the FEM.^{27,28} When $ \u03f5 \u2264 h w$, the phasefield method will construct a sharp fracture profile. In this situation, the fracture surface energy is overrated. Hence, we need to make $ \u03f5 > h w$ in the numerical simulation. However, $\u03f5$ cannot be chosen too large because large $\u03f5$ will produce softening effects, which can result in the undervaluation of the strain energy. In the meanwhile, fracture and cave are overly smeared to the intact region, which can interfere with the fracture–cave interaction.
The simulation result is shown in Fig. 8(b). We can find when $\u03f5$ = 1.5, 2, and 3h_{w}, the COD under the three cases is approximately the same. However, when $\u03f5$ = 4h_{w}, the computed COD is obviously different from others. We can also find there is a sudden increase in COD near the fracture–cave connection position. It is because the interaction is affected by the excessive diffusive interfacial thickness. From the simulation result, it is concluded that the reasonable choice of $\u03f5$ is h_{w} < $\u03f5$ < 4h_{w}. In other simulation examples of this paper, we make $\u03f5$ = 2h_{w}.
V. NUMERICAL CASES AND DISCUSSION
The existence of natural cave can markedly impact the HF propagation. The interaction mechanism between HF and natural cave is very complex. In this section, we set the following numerical cases to investigate the interaction mechanism in depth.
A. Effect of the distance between fracture and cave
To investigate the effect of the distance between fracture and cave, four cases are conducted here. The geometric structure of these cases is shown in Fig. 9. The computation domain is 20 $\xd7$ 20 m in size, and an initial fracture is set at its center. The coordinates of point 0 are (0, 0). The injection point is in the fracture center, and its coordinates are (10, 10). The position of cave's center is (14, H_{c}), while its radius is 1.3 m. To simulate the different distances between fracture and cave, H_{c} is set as 10, 12, 13.2, and 15.5. It is obvious that when H_{c} increases, the distance between fracture and cave increases. The initial reservoir pressure is 0.1 MPa. The maximum and minimum principal stresses are applied to the four boundaries of the computation domain. Other parameters are shown in Table II.
Parameters .  Value .  Unit . 

Young's modulus E  30  GPa 
Poisson's ratio v  0.2  / 
Initial fracture length L_{f}  0.6  m 
Cave radius r_{c}  1.3  m 
Biot coefficient $ \u2009 \alpha $  1  / 
Fluid injection rate Q  1 $\xd7$ 10^{–4}  m^{2}/s 
Maximum principal stress $ \sigma H$  0.5  MPa 
Minimum principal stress $ \sigma h$  0.4  MPa 
Porosity ϕ  0.05  / 
Fluid viscosity μ  3 $\xd7$ 10^{–4}  Pa·s 
Initial reservoir pressure p_{r}  0.10  MPa 
Rock bulk modulus K_{s}  10  GPa 
Fluid bulk modulus K_{f}  0.6  GPa 
Reservoir permeability k_{r}  $ 1 \xd7 10 \u2212 16$  m^{2} 
Critical energy release rate G_{c}  250  N/m 
Parameters .  Value .  Unit . 

Young's modulus E  30  GPa 
Poisson's ratio v  0.2  / 
Initial fracture length L_{f}  0.6  m 
Cave radius r_{c}  1.3  m 
Biot coefficient $ \u2009 \alpha $  1  / 
Fluid injection rate Q  1 $\xd7$ 10^{–4}  m^{2}/s 
Maximum principal stress $ \sigma H$  0.5  MPa 
Minimum principal stress $ \sigma h$  0.4  MPa 
Porosity ϕ  0.05  / 
Fluid viscosity μ  3 $\xd7$ 10^{–4}  Pa·s 
Initial reservoir pressure p_{r}  0.10  MPa 
Rock bulk modulus K_{s}  10  GPa 
Fluid bulk modulus K_{f}  0.6  GPa 
Reservoir permeability k_{r}  $ 1 \xd7 10 \u2212 16$  m^{2} 
Critical energy release rate G_{c}  250  N/m 
Figure 10 shows the simulation result of the HF propagation path. We can find the propagation path changes with the different distance. When H_{c} is 10, the HF propagates horizontally along the X direction (the maximum principal stress direction) and finally is connected with the cave. When H_{c} is 12, the HF first propagates horizontally and then deflects up near the location of cave. Finally, the HF is connected with the cave. When H_{c} is 13.2, the HF also deflects up and is connected with the cave. Yet, the location where HF begins to deflect is farther away from the cave compared with H_{c} is 12. When H_{c} is 15.5, the HF propagates along the X direction throughout. The cave does not affect the HF propagation path.
Through the simulation result, it is found the stress concentration caused by cave can play a notable role in the HF propagation. It can induce the HF to propagate away from the direction of maximum principal stress. However, when the distance between fracture and cave is very large, this stress concentration is negligible to the HF propagation.
The pressure distribution with different distance is shown in Fig. 11. As seen from Fig. 11(d), the HF is not connected with the cave. The fluid pressure in the fracture domain is maximum. Due to the fluid leak off, the pressure increases in the reservoir area, which is close to the fracture domain. However, the fluid pressure in the cave is still equal to the initial reservoir pressure. It is because the cave is far away from the HF, so the injection fluid does not enter the cave. As seen from Figs. 11(a)–11(c), the pressure in the cave significantly increases. It is because the HF is connected with the cave, resulting in highpressure injection fluid into the cave.
The fluid pressure at injection point is a key parameter of hydraulic fracturing. It can be measured in engineering practice and has an important guiding significance for the optimization design of hydraulic fracturing. In this work, we analyze the pressure variation at injection point in detail. When H_{c} is 13.2, the fluid pressure at injection point vs time is shown in Fig. 12. The pressure rapidly increases before the fracture initiation. When the HF propagates, the pressure gradually declines. However, when time is 12.0 s, the fluid pressure steeply declines. According to the phasefield evolution contour, we find the HF is connected with cave at this time. After the connection, the injection fluid enters the cave. The fluid pressure at injection point begins to gradually increase.
Figure 13 presents the fluid pressure vs time at varied H_{c}. When H_{c} is 15.5, the HF is not connected with the cave. Hence, the pressure curve does not experience a sharp drop. Meanwhile, the gradual growth stage is not existent in the pressure curve. Comparing the curve of H_{c} is 13.2 and H_{c} is 12, we can find the pressure experiences a sharp drop earlier when H_{c} is 12. It is because the distance between fracture and cave is shorter when H_{c} is 12. Hence, the HF is connected with the cave earlier. Finally, the fluid pressure in the three cases all tends to be stable.
B. Effect of injection rate
Injection rate is a critical factor of hydraulic fracturing operation. Engineers can improve the efficacy of hydraulic fracturing by adjusting the injection rate. Here, our study pays attention to how the injection rate affects the interaction between HF and cave. The position of cave center is (14, 12.7). The injection rate Q is set as 1.0 $\xd7$ 10^{−4}, 1.4 $\xd7$ 10^{−4}, 1.8 $\xd7$ 10^{−4}, and 2.2 $\xd7$ 10^{−4} m^{2}/s. Other parameters conform to Table II.
Figure 14 presents the phasefield evolution contour. We can find the propagation path varies with different injection rate. The interaction pattern can be divided into two types. When Q is 1.0 $\xd7$ 10^{−4} and 1.4 $\xd7$ 10^{−4} m^{2}/s, the HF is finally connected with the cave. When Q is 1.8 $\xd7$ 10^{−4} and 2.2 $\xd7$ 10^{−4} m^{2}/s, the HF deflects along the cave, but not connected with it. For a clearer comparison of the four cases, we use v = 0.4 to display the HF profile on the right side of injection point, which is shown in Fig. 15. Here, we define the deflection point, which represents the point where the HF begins to deflect up. It can be seen the deflection point is farther away from injection point when Q increases. In addition, deflection angle of HF deviation from X direction decreases as Q increases. The result reveals that the cave effect on HF declines when the injection rate increases. Hence, the HF is more inclined to propagate along the maximum principal stress direction.
Figure 16 shows the fluid pressure at injection point with different injection rates. We can find the fracture initiation pressure raises as the injection rate increases. When Q is 1.0 $\xd7$ 10^{−4} and 1.4 $\xd7$ 10^{−4} m^{2}/s, the fluid pressure experiences a sharp drop. It is because the HF is connected with cave. However, when Q is 1.8 $\xd7$ 10^{−4} m^{2}/s, the fluid pressure does not experience a sharp drop since the HF is unable to connect with cave.
C. Effect of cave radius
The cave radius r_{c} can also markedly impact the HF–cave interaction. In above section, we have found the stress concentration caused by cave can make the HF deflect. In this section, we devote more attention to investigate the impact on the propagation velocity of HF.
We conduct three cases here. In case 1, the cave radius is set as 0, which indicates there is no cave in the reservoir. In case 2 and case 3, the cave radius is set as 0.5 and 0.9 m, respectively. The coordinates of the cave center are (15, 10). Other parameters conform to Table II.
Figure 17 presents the phasefield evolution contour. When time is 6.0 s, we can find the fracture length is different in the three cases. In case 1, the fracture length is equal on both sides of injection point. However, in case 3, it is obvious that the fracture length on the right side of injection point is longer than that on the left side. In case 2, the length difference between the two sides is very small, with the rightside fracture length of 1.93 m and the leftside fracture length of 1.86 m. When time is 9.5 s, the HF is connected with the cave in case 3. The rightside fracture length is much longer than the leftside fracture length. When time is 19.0 s, the rightside fracture length does not increase in case 3, since the HF has been connected with cave, while the leftside fracture length has increased a little. In case 2, the rightside fracture length is 4.5 m, and the leftside fracture length is 3.9 m. There is a noticeable length difference between the two sides.
Figure 18 presents the rightside fracture length vs time. The HF begins to propagate when time is 2.5 s in all three cases. Then, it can be found the HF propagation velocity is different among the three cases. The HF propagates faster as the cave radius increases. Contrasting case 1 and case 2, we can find that the existence of cave can make the HF propagate faster. Moreover, if there is no cave, it is apparent that the propagation velocity will gradually decrease with time. However, if the cave exists, the propagation velocity can suddenly increase near the cave.
Some significant findings can be drawn from the simulation result. The stress concentration caused by cave is similar to an attraction force for the HF. If the cave location is far away from the maximum principal stress direction, it can make the HF deflect. If the cave location is in the direction of maximum principal stress, it can make the HF propagate faster. The propagation velocity raises as the cave radius increases. In addition, it can cause the asymmetric propagation of HF. The propagation velocity is greater on the side with cave.
D. Effect of in situ stress in threedimensional example
The in situ stress can alter the stress concentration caused by cave, which further influences the HF–cave interaction. Here, we set a threedimensional example to investigate the interaction under different in situ stress.
As discussed in Sec. III A, we use a structured eightnode element to discretize the spatial domain. Hence, our model is applicable to the threedimensional simulation. The schematic of the threedimensional example is shown in Fig. 19. There is a cubic computation domain, measuring 5 $\xd7$ 5 $\xd7$ 5 m in magnitude. Its center is placed a pennyshaped fracture, with 0.1 m radius. The coordinates of point 0 are (0, 0, 0). Hence, the coordinates of fracture center point are (2.5, 2.5, and 2.5). In addition, there is natural cave at the upper right of the pennyshaped fracture. The cave center is located at (3.5, 2.5, and 3.3), and the cave radius is 0.4 m. The pennyshaped fracture is pumped with fracturing fluid at a steady rate of $ 5 \xd7 10 \u2212 4 m 3 / s$. To investigate the effect of in situ stress, we set three cases here: case 1: σ_{x} = 1, σ_{y} = 1, and σ_{z} = 0.8 MPa; case 2: σ_{x} = 1, σ_{y} = 1, and σ_{z} = 0.6 MPa; case 3: σ_{x} = 1, σ_{y} = 1, and σ_{z} = 0.2 MPa. Other parameters are identical to those in Table II. In this simulation example, the element size is $ 0.05 \xd7 0.05 \xd7 0.05 \u2009 m ,$ and there are 1 030 301 nodes in total. We use the Intel Xeon Platinum 9242 CPU with 96 cores for parallel computing. The case 1 is computed within 12 h for 4000 time steps.
The propagation path of the pennyshaped fracture is presented in Fig. 20. We use the phasefield v = 0.25 to display the threedimensional fracture shape and cave shape. In case 1, the in situ stress difference $ \u2009 \sigma d = \sigma x \u2212 \sigma z = 0.2 \u2009 MPa .$ In cases 2 and 3, the in situ stress difference is 0.4 and 0.8 MPa, respectively. It is generally believed that the pennyshaped fracture propagates in the direction orthogonal to the minimum in situ stress. However, in case 1, the pennyshaped fracture deflects up. It is because the stress concentration caused by cave is similar to an attraction force, which can change the fracture propagation path. In case 2, the pennyshaped fracture also deflects up. However, the deflection angle of pennyshaped fracture deviation from x–y plane decreases. In case 3, the pennyshaped fracture propagates horizontally, and it does not deflect. From the three cases, we can find the impact of cave stress concentration declines as the in situ stress difference increases. Hence, the pennyshaped fracture is more inclined to propagate in the direction orthogonal to the minimum in situ stress.
Next, we test the capacity of our model in handling threedimensional issues when multiple caves exist. We set two cases for evaluation. In case 3, there are two caves in the computation domain, with 0.4 m radius. Their locations are (3.5, 2.5, and 3.2) and (1.5, 2.5, and 1.8), respectively. In case 4, we add a new cave based on case 3. The new cave is located at (3.5, 2.5, and 1.8), and its radius is 0.25 m. The in situ stress of both case 3 and case 4 is the same as case 1. The fracture propagation path is presented in Fig. 21. When t = 100 s, the fracture predominantly propagates horizontally. Comparing Figs. 21(a) and 21(c), we can find the fracture propagates farther to the right side in case 4. It is because the two caves on the right side create larger attraction force in case 4. When t = 400 s, the fracture is connected with the left and right caves in case 3. However, in case 4, the fracture is not connected with the right cave. It only deflects up a little.
The above example substantiates the capacity in our model to simulate the HF propagation with natural cave in threedimensional issues. Moreover, we use the fixed mesh to characterize the natural cave and simulate the HF propagation. The proposed model does not need any special node treatment for characterizing the cave. In addition, it does not need any remeshing (the element size is smaller than $\u03f5$ in the whole computation domain). These advantages fully prove the conciseness and efficiency of the proposed model.
VI. MULTIPLE CAVE SIMULATIONS
Sometimes, there are multiple natural caves in the geomaterials. These caves are distributed randomly and vary in size. The following example is presented to check the performance of our model in multiple cave issues. The geometric structure of this example is shown in Fig. 22. The computation domain is 20 × 20 m^{2} in size, and its center is placed an initial fracture. The initial fracture length is 0.6 m, and the angle between the fracture and horizontal direction is 25°. The four outer boundaries of computation domain are fixed. The injection rate of fracturing fluid is 1 $\xd7$ 10^{−4} m^{2}/s. Within the computation domain, there are six natural caves and their radii range from 0.85 to 1.1 m. In addition, the distances between the initial fracture and caves are 3.9–5.5 m. Other parameters are identical to those in Table II.
Figure 23 shows the phasefield evolution contour. We can found the HF propagation path is significantly influenced by the multiple caves. When t is 2.5 s, the HF begins to extend. After 1.5 s, the HF has evidently grown. Because of the stress concentration caused by caves, the propagation path is not straight. When t is 8.1 s, the HF is connected with the lowerleft cave. Since the lower two caves are larger than the up two caves, the attraction force of the lower two caves is larger. Therefore, the HF propagates faster along the lowerleft direction. Then, the HF continues to extend and is connected with the upright cave when t is 20.3 s.
Figure 24 shows the fluid pressure variation at injection point over time. The pressure reaches its peak when fracture initiation. After that, the pressure gradually decreases. Subsequently, we can found the pressure experiences two sharp drops. The reason is that the HF is connected with the caves at these two time points. The finding is consistent with the study in Sec. V.
This example demonstrates that our model can simulate the HF propagation in the geomaterials containing multiple caves. It also proves the wide applicability of our model.
VII. CONCLUSIONS
In this study, we establish a novel hydraulic fracturing model for the fluiddriven fracture propagation in poroelastic media containing the natural cave by the phasefield method. First, we derive the fracturecavereservoir flow governing equations. In reservoir domain, the fluid flow adheres to lowpermeability Darcy flow; in cave domain, it follows free flow, while in fracture domain, it conforms to generalized Reynolds flow. The hydromechanical coupling is achieved by employing Biot poroelasticity theory and updating fracture width. The phase field is introduced to diffuse the sharp fracture and cave edge. Compared to existing study, our model does not need additional enrichment function to characterize the cave. The proposed model can naturally detect and simulate the hydraulic fracture (HF) propagation behavior in the existence of natural cave. The complicate calculation of stress intensity factor when HF approaches cave is avoided. The effectiveness of the proposed model in twodimensional issue is fully displayed. Then, the model is extended to the threedimensional issue and multiple cave issues.
The interaction mechanism between HF and natural cave is also investigated in accordance with the proposed model. Several findings are obtained: (1) The stress concentration caused by cave can affect the HF propagation. If the cave is far away from the maximum principal stress direction, it can make the HF deflect toward the cave. (2) The cave effect can change the HF propagation velocity. The fracture propagates faster on the side with cave. (3) Increasing the fracture–cave distance, injection rate, and in situ stress difference can all decline the cave effect. (4) The cave effect raises as the cave radius increases. (5) The fluid pressure at injection point will experience a sharp drop if the HF is connected with cave. After the connection, the fluid pressure will slowly increase.
ACKNOWLEDGMENTS
The authors gratefully acknowledge the funding support from Fourteenth FiveYear Major Science and Technology Project of CNPC (Grant No. 2021DJ1505).
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
Jie Jin: Data curation (lead); Investigation (lead); Validation (lead); Writing – original draft (lead); Writing – review & editing (lead). Xiaoqiang Wang: Formal analysis (lead); Methodology (lead); Software (lead). Xiaohua Liu: Funding acquisition (lead); Resources (equal). Yanmei Xu: Resources (equal); Supervision (lead). Detang Lu: Conceptualization (lead); Project administration (lead).
DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.
APPENDIX: THE DISCRETIZING PROCESS OF THE FLOW EQUATIONS
Here, we define a indicator function $ v c$ to distinguish between cave domain and noncave domain. Specifically, it has the same form as $v$ in Eq. (31). In noncave domain, $ v c = ( 0 , 1 ]$. However, in cave domain, $ v c = 0$. Let $ \Omega d = \Omega r \u222a \Gamma $, it represents the noncave domain (reservoir domain and fracture domain). First, the integral in $\Gamma $ and $ \Omega r$ needs to be transformed to the integral in $ \Omega d$.
Chukwudozie et al.^{34} proposed a regularized method to transform the integral in $\Gamma $. It is realized through the integral term in fracture domain $ \Gamma $ multiplied by $ \u2009  \u2207 v  .$ Following their method, the integral in $ \Omega r$ is transformed through the integral term in $ \Omega r$ multiplied by $ v 2$. To avoid the illconditioned system when $ v = 0$, they proposed that Biot modulus term can be used as a stabilizing term.