A Numerical Study on the Propagation Mechanisms of Hydraulic Fractures in Fracture-Cavity Carbonate Reservoirs

Field data suggests that carbonate reservoirs contain abundant natural fractures and cavities. The propagation mechanisms of hydraulic fractures in fracture-cavity reservoirs are different from conventional reservoirs on account of the stress concentration surrounding cavities. In this paper, we develop a fully coupled numerical model using the extended finite element method (XFEM) to investigate the behaviors and propagation mechanisms of hydraulic fractures in fracture-cavity reservoirs. Simulation results show that a higher lateral stress coefficient can enhance the influence of the natural cavity, causing amore curved fracture path. However, lower confining stress or smaller in-situ stress difference can reduce this influence, and thus contributes to the penetration of the hydraulic fracture towards the cavity. Higher fluid viscosity and high fluid pumping rate are both able to attenuate the effect of the cavity. The frictional natural fracture connected to the cavity can significantly change the stress distribution around the cavity, thus dramatically deviates the hydraulic fracture from its original propagation direction. It is also found that the natural cavity existing between two adjacent fracturing stages will significantly influence the stress distribution between fractures and is more likely to result in irregular propagation paths compared to the case without a cavity.

During the past decades, based on different numerical approaches like the finite element method (FEM) [3,4], the displacement discontinues method (DDM) [5], the phase field method [6,7], the distinct element method (DEM) [8], the extended finite element method (XFEM) [9][10][11][12], and the proper generalized decomposition method (PGD) [13], a variety of numerical models have been established and applied by researchers to study hydraulic fracture propagation in consideration of different kinds of influence factors. Among these studies, Schrefler et al. [3] proposed an adaptive refinement technique to simulate hydraulic fracturing problems based on the generalized finite element formulation. Song et al. [4] performed a series of numerical simulations and examined the influence of some key factors on hydraulic fracturing propagation using RFPA2D-Flow. Zhou et al. [7] proposed a phase field model for fluid-driven dynamic crack propagation simulation in the poroelastic media. Wang et al. [13] numerically studied the hydraulic fracturing problems via the PGD method. Recently, the interaction mechanisms between hydraulic fractures (HFs) and pre-existing natural fractures (NFs) in reservoirs have increased markedly [11,14,15]. When an HF meets an NF, it can be deflected or arrested by the NF, active the NF and branch along one or two sides of the NF, or even bypass the NF [14], leading to very complicated patterns of HF-NF interactions. Readers are referred to [14,15] for a comprehensive review of this topic. However, despite these efforts, a systematic study of mechanisms and factors that govern the propagation behavior of hydraulic fractures around natural cavities has been seldomly reported in the literature. Hence, an attempt is made in this study to present a comprehensive investigation of the interaction mechanism of hydraulic fracture and naturally-formed caves and fractures by establishing a reliable numerical model.
The difference in the propagation mechanisms of hydraulic fractures between fracture-cavity reservoirs and conventional reservoirs is mainly caused by the stress concentration around cavities. As shown in Fig. 1, the analytical solution of stress in the tangential direction, σ θ , surrounding the cavity of radius r can be written as [16] where λ = σ h /σ H denotes the lateral stress coefficient. From Eq. (1), it can be easily concluded that if λ < 1/3, tensile stress zones exist at the left and right sides near the cavity, thus promoting the propagation of a nearby hydraulic fracture towards the cavity located ahead of the fracture tip. If 1/3 < λ < 1/2, although the tensile stress zone disappears, σ θ < σ H can still be satisfied around the left and right sides of the cavity, which is also favorable for hydraulic fracturing. However, if λ > 1/2, the hydraulic fracture will be arrested or deflected from its original propagation path due to the compressive stress concentration phenomenon. With this in mind, therefore, it can be conjectured that once the cavities are considered in reservoirs, quite different physical mechanisms and propagation-path scenarios of hydraulic fracturing can be obtained. For instance, Liu et al. [6] investigated the effects of full-filled cavities on hydraulic fracturing based on the phase field method. They found that the propagation paths are sensitive to Young's modulus ratio between rock formation and cavity, spatial location of the cavity as well as differential in-situ stress. Likewise, Luo et al. [17] established a coupled seepage-free flow-mechanical model and studied the influence of the fluid pressure inside a cavity as well as the perforation direction on the hydraulic fracture propagation paths. It was found that higher fluid pressure in the cavity is more likely to reduce fracture path deflection and the perforation direction has a profound influence on the fracture propagation path. Cheng et al. [11] built a numerical model and studied the influence of cavity location as well as the difference of in-situ stress on the propagation path of a hydraulic fracture. He et al. [18] proposed a coupled hydromechanical model and studied the effects of existing pressurized cavities on the propagation of fluid-driven fracture. Besides, a numerical simulator called FEMM-FracFlow was used by Wang et al. [19] to investigate the behaviors of hydraulic fractures and found that fracture paths are closely related to the cavities, pumping rate of fracturing fluid, and the in-situ stress level. In addition to hydraulic fracturing technology, acid fracturing is one of the major measures of reservoir stimulation, especially for carbonate reservoirs. Zhao et al. [20] established a numerical tool for acid fracturing simulation and studied the mechanism of hydraulic fracture propagation in cases involving a single cavity and a pre-existing natural fracture. Simulation results [21] suggest that a smaller approaching angle of the hydraulic fracture towards the natural fracture and longer natural fracture benefit the penetration of hydraulic fracture into the cavity. As another alternative, the explosive pulsed fracturing technique is receiving more attention from researchers [21]. Wang et al. [22] developed a numerical model to investigate the pulsed fracturing process in fracture-cavity reservoirs and found that natural fracture can potentially guide the pulsed fracturing path under the action of pulse pressure. Although some works aiming at explaining the mechanisms of interaction among HF, NF, and cavity have been done by researchers, it still lacks a systematic study and deeper understanding of crucial factors that influence the hydraulic fracturing efficiency in fracture-cavity reservoirs. On the other hand, existing studies mainly focus on simple geometrical configuration, and the study of the influence of cavity on hydraulic fractures in a wellbore is not available in the literature. In this study, the XFEM proposed by Belytschko et al. [23] is adopted to describe the displacement field around fractures and cavities, the intersection between an HF and an NF, and the penetration of a fracture into a cavity. Consequently, the tedious remeshing process can be avoided after the propagation of fractures [10]. The fluid flow within fractures is described by Reynold's equation [9] which is discretized using the FEM. Afterward, the displacements of all solid nodes including enriched nodes, and the pressure distribution of all fluid nodes along the hydro-fracture are obtained by solving the fully-coupled governing equations [10] using the Newton-Raphson (N-R) method. This paper is presented as follows. The description of the numerical model is given in Section 2. Implementation and validation of the proposed numerical model are presented in Section 3. Mechanisms and factors that govern the propagation behavior of hydraulic fractures around natural cavities are systematically studied in Section 4. Major conclusions are summarized in Section 5.

Description of the Numerical Model
A shown in Fig. 2, we consider the domain with a hydraulic fracture HF , a frictional [24] natural fracture NF , and a cavity inside. The boundary of the domain is with a normal vector n Γ .
is fixed at boundary Γ u and traction t is applied at boundary Γ t . s represents the coordinate system defined along the hydraulic fracture. The incompressible Newtonian fluid is pumped at a rate of Q inj into HF . To make the numerical model more tractable, the following assumptions are made. The fracture growth process is treated as quasi-static [25]. The fluid lag phenomenon [14] is ignored. Natural cavities usually exist in three forms, i.e., non-filled, partlyfilled, and fully-filled cavities [2]. Non-filled cavities are studied in this article. In addition, although natural cavities are of quite irregular shapes, they are treated as regular circles or ellipses in this article for the sake of modeling convenience.

Governing Equations
According to the linear elastic assumption of rock formation, the governing equations for quasi-static deformation can be written as [10] where σ denotes the Cauchy stress tensor, t cont represents the traction vector caused by contact behaviors between fracture surfaces.
The fluid flow along the hydraulic fracture can be described by Reynold's equation [25] ∂w where w, t, p, and μ represent fracture aperture, time, fluid pressure, and fluid viscosity, respectively. The fracture aperture w can be calculated by w = [[u]] · n Γ HF , where [[u]] denotes the displacement jump across the fracture. The term g L represents the fluid leak-off from the hydraulic fracture into the surrounding rock formation, which can be written as [9,26] where C L is the fluid leak-off coefficient, and t 0 (s) represents the time instant when fluid starts to leak for a point s along the hydraulic fracture.
After introducing the trial function u(x, t) and test function δu(x, t), we can obtain the weak form of Eq. (2): Besides, by introducing the test functionδp(s, t) and integrating by parts, the weak form of Eq. (3) can be written as: Spatial and time discretization of Eqs. (5) and (6) can be found in our previous work [10]. Afterward, the solid-fluid coupling problem is governed by Eqs. (2) and (3) can be iteratively solved by using the Newton-Raphson method [10].

Extended Finite Element Method
In this paper, we use the XFEM to describe the displacement field around fractures and cavities. By simply introducing additional degrees of freedoms (DOFs), i.e., the enriched DOFs, the remeshing and data mapping between old and new meshes can be avoided [27] and the demand on mesh refinement around fracture tip can be greatly alleviated owing to the fracture tip enrichment function, leading to a significant reduction of the problem complexity and computational time. This is a substantial improvement for numerical simulation of problems containing strong or weak discontinuities, especially for simulation cases where a large number of fractures and caves are involved. Since its introduction in 1999 [23,27], the XFEM has been adopted by many researchers to study various kinds of problems related to hydraulic fracturing. For problems where the intersection of an HF and an NF, as well as the intersection of a fracture and a cavity, are considered, the displacement of a point x inside the domain can be written as in which N I denotes the shape function of node I, S all denotes the set of all nodes whose displacement vector is denoted by u I .
where (r, θ) represents the local polar co-ordinates defined at fracture tips [27]. As depicted in Fig. 3, in addition to the Heaviside enrichment nodes, the junction enrichment nodes are necessary to describe the intersection of two fractures. The T-shaped junction considered in this article, J (x) takes different values of 1, −1, and 0 in different subdomains [10]. Besides the Heaviside enriched nodes and cavity enriched nodes, the penetration enriched nodes are presented to describe the penetration of a fracture into the cavity, just as illustrated in Fig. 4. The penetration enrichment function P (x) takes the following form which indicates that P (x) equals 0 if x lies inside the cavity, and otherwise equals 1 or −1 in accordance with the Heaviside enrichment function.

Frictional Behaviors of Natural Fractures
Frictional and cemented fractures are two common types of natural fractures in reservoirs [24]. In this article, the frictional natural fractures are considered due to their stronger influence on the creation of the fracture network in comparison with cemented natural fractures [24]. Thus, to avoid embedding of fracture surfaces under the effect of compressive stress, the frictional behaviors of natural fractures must be considered. The no-embedding conditions of fracture surfaces can be expressed as [28] where t cont T represents the tangential component of t cont , φ f represents the friction angle of the fracture surface, and c f is the cohesive strength between fracture surfaces. If F f < 0, then the fracture surfaces are in a bonding state. If F f = 0, then the fracture surfaces are in a slipping state. In this article, the penalty function method together with the Newton-Raphson method [10] is adopted to consider the frictional force and the contact status in the simulation. For each Newton-Raphson iteration Step i, the following linear system is solved and the global conventional DOFs vector u and enriched DOFs vector a are updated according to until the following convergence criterion is satisfied In Eq. (12), K uu , K ua , K au , K aa , 1 and 2 take the following forms, respectively In Eq. (15), N is the matrix of shape function, B is the matrix of shape function derivatives, D is the stress-strain matrix, and D cont is related to contact status of fractures and can be written as [10] where k N and k T are the penalty parameters in the normal and tangential directions, respectively, and m Γ NF is the unit vector in the tangential direction of the natural frictional fracture. In Eq. (15), the frictional force between fracture surfaces is considered by the last term F contact of 2 .

Fracture Propagation Criterion
The widely adopted maximum hoop-stress criterion [30] related to Mode-I (K I ) and Mode-II (K II ) stress intensity factors (SIFs) is used to check whether and along which direction the fracture will propagate where α = 2 arctan If K e exceeds the fracture toughness of reservoirs (K IC ), the fracture will propagate and deflect by the angle α compared to its original direction. The SIFs are calculated by using the interaction integral approach [31]. When an HF meets a frictional NF, whether the HF will propagate along the NF, cross the NF, or be arrested by the NF is determined using the extended Renshaw and Pollard criterion [32]. Implementation details of this criterion can be found in our previous article [24].

Model Implementation and Validation
The proposed numerical model is implemented in an in-house Fortran code called PhiPsi (http://phipsi.top/): A general-purpose XFEM-based program. The level set method [33] is used to describe fractures and cavities, and to track the growth of fractures. The compressed sparse row (CSR) format is used to store global system matrices with optimized memory consumption. For the simulation cases where frictional natural fractures are considered (Section 4.1.6), two types of Newton-Raphson iteration loops are performed [10]. The outer loop solves the solid-fluid coupling equations (Eqs. (2) and (3)), and the inner loop determines the contact status of natural frictional fractures [10]. For each Newton-Raphson iteration step, an iterative solver called Lis [34] is adopted to solve the linear system. The flowchart of the simulation algorithm [10] is shown in Fig. 5.
The validation of the presented numerical model without considering cavities has been sufficiently conducted in our previous articles [9,10,24,35]. In this section, the propagation of an initial fracture towards a cavity in a panel under the action of tensile load is studied to verify the capacity of the proposed model for predicting interactions between fractures and cavities. As illustrated in Fig. 6, the bottom of a plane stress panel is fixed and the top of the panel is subjected to tensile stress σ . The geometric parameters H, h 1 , h 2 , l, b 1 , b 2 , and a are taken as 3.0, 1.5, 0.1, 1.2, 0.7, and 0.2 m, respectively. Besides, the values of Young's modulus and Poisson's ratio are given as 30 GPa and 0.2, respectively. The tensile stress σ is 5 KPa. The discretized model consists of 3,400 quadrilateral elements in total. The fracture propagates by an increment of 0.045 m and penetrates the cavity after 15 propagation steps. The resulting vertical displacement field of the panel and the enriched nodes are presented in Fig. 7 where six penetration enriched nodes can be seen. This problem has been studied by Nguyen et al. [36] using the extended meshfree Galerkin method. The comparison of the calculated propagation path with [36] is visualized in Fig. 8. As expected, the propagation path predicted by the proposed method is well consistent with that predicted by other numerical method, indicating that the interactions between fracture and cavity can be well captured using the proposed model.      Fig. 10, from which good agreements can be seen. In order to make a more accurate comparison, the relative shear displacement (displacement in the tangential direction) between fracture surfaces obtained from both the XFEM and the FEM solutions are presented in Fig. 11. Again, good agreements can be observed, which means that the proposed scheme to simulate contact behaviors [37] of frictional fractures is accurate and reliable.

Interactions between Hydraulic Fracture and Cavity
In this section, interactions between hydraulic fracture and cavity will be studied via a simple model shown in Fig. 12. The influence of factors including in-situ stress, viscosity of the injected fluid, pumping rate of the fluid, size of the cavity, shape of the cavity, and natural fracture will be thoroughly investigated. The size of the model is 5 m × 10 m. The radius of the wellbore is 0.1 m. The length of the initial hydraulic fracture is 0.5 m. The radius of the cavity denoted by a is set to 0.  The strongly deflected fracture propagation path and the stress distribution in y-direction are shown in Fig. 13. It can be conjectured that if there is no cavity in this example, the hydraulic fracture will propagate horizontally, i.e., in the direction orthogonal to the minimum principal stress σ h . However, as can be seen in Fig. 13, because of the strong stress concertation surrounding the cavity, the hydraulic fracture propagates and bypasses the cavity. In the meantime, the injection pressure at the wellbore required to drive the hydraulic fracturing process changes in a significantly different manner from the scenario without the cavity, as shown in Fig. 14. Specifically, higher fluid pressure is required at the beginning of the simulation because of the intense compressive stress (in y-direction) at the left side of the cavity, and the required fluid pressure drops quickly once the fracture grows into the area of low compressive stress above the cavity.

Effects of In-situ Stress
The level of in-situ stress increases with the increase of the reservoir depth. On the other hand, as discussed in the introduction section, the stress field surrounding the cavity is strongly affected by the lateral stress coefficient λ. Hence, three different cases of in-situ stress are investigated in this section, as shown in Tab. 1. In Case 1, we choose a higher level of in-situ stress: σ h is 90 MPa and σ H is 100 MPa. In Case 2, we choose a lower lateral stress coefficient: σ h is 10 MPa and σ H is 50 MPa. In Case 3, we choose the same lateral stress coefficient as Case 2 but low in-situ stress difference: σ h is 1 MPa and σ H is 5 MPa. The simulated fracture propagation paths are presented in Fig. 15. It can be found that both the lateral stress coefficient and the level of confining stress (or in-situ stress difference) have a strong influence on propagation paths. Specifically speaking, a higher lateral stress coefficient can enhance the influence of the natural cavity, causing a more curved fracture path. However, lower confining stress or smaller in-situ stress difference can reduce the influence, and thus contributes to the penetration of the hydraulic fracture towards the cavity.

Effects of Fluid Viscosity
Fluid viscosity is an important factor in hydraulic fracturing treatments. In this section, a fluid of higher viscosity, 0.1 Pa · s, is injected into the wellbore. The simulation results are presented in Fig. 16, from which it can be concluded that increasing the fluid viscosity can reduce the impact of the cavity on the propagation path. The main reason for this is that the fluid pressure inside the hydraulic fracture is higher for higher viscosity fluid [38], and the higher fluid pressure can counteract the effect of the stress concertation caused by the cavity.

Effects of Fluid Pumping Rate
In this section, we investigate the effect of another key factor, the fluid pumping rate. Different from the base simulation case, a higher pumping rate, 0.01 m 2 /s, is applied. The simulated propagation path is also shown in Fig. 16 in which a smoother path can be seen compared to the base case. Therefore, it can be concluded that a high fluid pumping rate is able to attenuate the effect of the cavity, and thus benefits the propagation of the hydraulic fracture along its original direction. This phenomenon can be explained with the aid of the analytical solution of the KGD model [38] which reveals that a higher pumping rate results in increased fluid pressure.

Effects of Cavity Size
Cavities with greatly different scales are widely distributed in fracture-cavity carbonate reservoirs [2]. In this section, we consider a smaller cavity of radius 0.3 m to study the effect of cavity size. The simulated path is illustrated in Fig. 17. Just as expected, for a smaller cavity, its effect is significantly weakened in comparison with the base case. Besides, it is worth noting that the hydraulic fracture does not start to deflect until its tip comes into the zone near enough to the cavity.

Effects of Cavity Shape
Natural cavities are of quite irregular shapes in fracture-cavity reservoirs [2]. In this section, two elliptical cavities oriented in different directions (a horizontal elliptical cavity and a vertical elliptical cavity) are studied to investigate the effect of cavity shape. The elliptical and circular cavities have the same area. The eccentricity of both ellipses takes a value of √ 3/2. From Fig. 18 we can conclude that cavities of irregular shape have a stronger influence on the propagation path than regular cavity and the horizontal elliptical cavity leads to a more twisty propagation path. The stronger influence of the irregular cavity is caused by its higher stress concentrations compared to the regular circle cavity.

Effects of Natural Fractures
As shown in Fig. 19, the influence of natural fractures will be investigated through four different cases. All of the natural fractures have the same length of 0.6 m and are positioned along the normal direction of the cavity. The friction angle of the fracture surface and the cohesive strength of the natural fracture are π /6 and 0.5 MPa, respectively. k N and k T are taken as 10 4 GPa/m. In Figs. 19a and 19c, natural fractures are connected to the cavity; In Figs. 19b and 19d, however, natural fractures are disconnected from the cavity with a small distance of 0.1 m. The angle of natural fractures in Cases (a) and (b) is set to π /4 with respect to the x-axes, and in Cases (c) and (d) the angle is taken as −π /4. Simulation results are depicted in Figs. 20 and 21. It can be seen that the connectivity between existing natural fracture and cavity is a very important factor that influences the propagation behavior of the hydraulic fracture. The connected natural fracture dramatically deviates the hydraulic fracture from its original propagation direction. However, for the disconnected natural fracture, the influence is limited since no sliding occurs between the surfaces of the frictional fracture. In Case (d), a T-shaped junction is formed since the hydraulic fracture is arrested by the natural fracture.

Effects of Cavity on Hydraulic Fracturing in a Wellbore
This section is aimed to study the influence of cavity on hydraulic fracturing in a wellbore. As illustrated in Fig. 22, for the sake of simplicity, only the first two stages (Stages 1 and 2) are considered [9]. The size of the model is 100 m × 100 m. The spacing between stages is set to 10 m. The length of the initial hydraulic fracture is 2.5 m. The radius of the cavity denoted by a is set to 2.0 m. The vertical distance between the cavity and the fracture denoted by d is taken as 3 m. The in-situ stresses σ h and σ H are 20 and 30 MPa (the lateral stress coefficient λ equals 2/3), respectively. The Young's modulus E is 30 GPa, Poisson's ratio ν is 0.2, and fracture toughness K IC of the rock formation equals 1.0 MPa · m 1/2 . The viscosity μ and pumping rate Q inj of the injected fluid is 0.01 Pa · s and 0.001 m 2 /s, respectively. The fluid leak-off coefficient C L is taken as 10 −5 m/s 1/2 . The proppant volume concentration of the injected fluid [39] is set to 0.3. The distribution of proppant along the hydraulic fracture is obtained according to the method proposed in our previous papers [9,35]. Then, the width of the propped fracture (w p ) can be determined by [9,35] w p = w o c/η, where w o denotes the fracture opening before significant pressure drops, c is the proppant volume concentration, and η = 0.74 is the proppant packing density [40]. Afterwards, the penalty function method [9] is adopted to simulate the propped fracture, namely, the Stage-1 hydro-fracture in this example. The model is discretized into 10,281 quadrilateral elements. The fractures of both stages propagate by an increment of 0.85 m for each propagation step. It should be noted that simulation without a cavity is also performed for comparison purposes.
The propagation paths of Stages 1 and 2 fractures for both cases without and with cavity are illustrated in Fig. 23. The stress contours in x-direction for cases without and with cavity are shown in Figs. 24a and 24b, respectively. These simulation results reveal that the existing cavity between two sequential fracturing stages significantly influences the propagation paths. For the case without the cavity, the Stage-2 hydraulic fracture curves away from the straight Stage-1 fracture because of the stress shadow effects [41] caused by the propped Stage-1 fracture, and the final paths are symmetrical about the wellbore. For the case with cavity, the Stage-1 fracture deviates leftward for both the upper and lower fracture tips, and the fracture path below the wellbore (y < 0) is more curved than that above the wellbore (y > 0). For the subsequent Stage-2 fracture, its fracture tip above the wellbore (y > 0) grows away from the previous Stage-1 fracture; However, the fracture tip below the wellbore (y < 0) firstly grows towards and then away from Stage-1 fracture. Besides the propagation paths, the stress fields are also quite distinct. The maximum stress values in x-direction for the cases without and with cavity are 36.1 and 43.2 MPa, respectively. For the case shown in Fig. 24b, both the maximum stress and the minimum stress in x-direction occur around the cavity. The maximum stress occurs on the upper and lower sides of the cavity, and the minimum stress occurs on the left and right sides of the cavity. The distinct propagation paths shown in Fig. 23 are a direct result of the stress concentration caused by the cavity.

Conclusions
In this paper, we established a fully-coupled numerical model to investigate the mechanisms of hydraulic fractures in fracture-cavity reservoirs using the XFEM. The Heaviside, cavity, fracture tip, T-shaped junction, and penetration enrichment functions are proposed to describe the displacement jump across the fracture surface, displacement discontinuity over the cavity boundary, singular displacement field near the fracture tip, intersection between the hydro-fracture and the natural fracture, and penetration of a fracture into the cavity, respectively. Hence, tedious remeshing can be avoided. The fluid flow within fractures is described by Reynold's equation which is discretized using the FEM. Afterwards, the fully-coupled governing equations are solved iteratively using the Newton-Raphson method. After the validation of the proposed model in Section 3, several cases are simulated to investigate the effects of factors such as in-situ stress, fluid viscosity, fluid pumping rate, cavity size and shape, and natural fractures in Section 4.1. Besides, the effects of a cavity on the sequential hydraulic fracturing in a wellbore are studied in Section 4.2. According to the cases studied in this paper, the major conclusions can be reached as follows: (1) Both the lateral stress coefficient and the level of confining stress (or in-situ stress difference) have a strong influence on propagation paths of hydraulic fractures near cavities. A higher lateral stress coefficient can enhance the influence of the natural cavity, causing a more curved fracture path. However, lower confining stress or smaller in-situ stress difference can reduce this influence, and thus contributes to the penetration of the hydraulic fracture towards the cavity. (2) The fluid viscosity and fluid pumping rate are two dominant factors on the propagation path in hydraulic fracturing treatments when natural cavities are considered. Higher fluid viscosity and high fluid pumping rate are both able to attenuate the effect of the cavity, and thus benefit the propagation of the hydraulic fracture along its original direction.
(3) The influence of a cavity depends not only on its size but also on its shape. Cavities of irregular shape (ellipse, for example) have a stronger influence on the propagation path of hydraulic fracture than regular circle cavity. (4) The frictional natural fracture disconnected from the cavity, even with a very small distance between the fracture tip and the cavity, has limited influence on the stress field around the cavity. Nevertheless, a frictional natural fracture connected to the cavity can significantly change the stress distribution around the cavity, thus dramatically deviates the hydraulic fracture from its original propagation direction. (5) Natural cavity existing between two adjacent fracturing stages will significantly influence the stress distribution between fractures and is more likely to result in irregular propagation paths compared to the case without cavity.

Availability of Data and Materials:
The data that support the findings of this study are available from the corresponding author on reasonable request.