Effect of Vugs on Hydraulic Fracture Propagation with Phase Field Method

Department of Mechanics, College of Energy and Mining Engineering, Shandong University of Science and Technology, Qingdao 266590, China State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation, Southwest Petroleum University, Chengdu 610500, China Research Center of Multiphase Flow in Porous Media, China University of Petroleum (East China), Qingdao 266580, China Bohai Oilfield Research Institute, Tianjin Branch, CNOOC China Limited, Tianjin 300459, China


Introduction
Carbonate rock is characterized by low porosity, low permeability, and complicated pore throat structure. However, natural fractures and vugs are abundant and dispersed inside the rock. Through hydraulic fracturing, these natural fractures and vugs can be connected together to form a favorable channel for oil and gas flow [1]. There are lots of researches that have been done to investigate the interaction between hydraulic fracture and preexisting natural fractures [2][3][4]. However, few studies are focusing on the interaction between hydraulic fracture and vugs. It still remains unclear how hydraulic fracture propagates when it encounters the vugs.
Experimental and numerical studies have been done on hydraulic fracture propagation with vugs. Liu et al. [5] conducted experiments to investigate the effect of vugs on hydraulic fracture propagation by using true triaxial apparatus and identified three main interaction modes; that is, the hydraulic fracture may cross, bypass, or be arrested by the vug. Guo et al. [6] carried out a numerical simulation to examine the effect of dissolved cavern on hydraulic fracture by using the seepage stress damage coupling model, and some influencing factors including cavern property, reservoir parameter, and treatment parameters were analyzed. Luo et al. [7,8] presented numerical simulation to investigate the effect of caverns on hydraulic fracture propagation direction under triaxial stress condition based on the extended finite element method and concluded that fractures tend to deflect from the maximum horizontal principal stress direction due to caverns. The effect of holes on crack propagation is also studied in other research areas [9,10].
Regarding hydraulic fracturing modeling, there are many different kinds of numerical methods that have been proposed. The commonly used methods include the finite element method, boundary element method, and discrete element method [11][12][13][14][15]. The finite element method is usually combined with the cohesive zone model for the hydraulic fracturing problem [16]. Fractures propagate along element boundaries, and thus, mesh refinement around fractures is needed for accuracy due to stress singularity of fracture tips. Instead, the extended finite element method is a good alternative, which allows fracture propagation across elements without remeshing [17,18]. However, it becomes cumbersome when dealing with intersection of multiple fractures [19]. The boundary method only discretizes the fracture surface boundary, which results in much less elements and lower computational cost. The displacement discontinuity method is a kind of indirect boundary element method, which has been extensively employed for simulation of complex fracture propagation in unconventional reservoirs [20,21]. Although the boundary element method has the advantage of fast computational velocity, it is difficult to consider the rock anisotropy and heterogeneity. The discrete element method treats rock as an assembly of discrete particles. The movement of particles obeys Newton's Second Law. It is easy to simulate complex fracture propagation by using the discrete element method, but the computational cost is high [22,23]. Moreover, much progress has been made in considering multiphysics coupling effect and rock plastic strain on fracture propagation [24][25][26].
The above methods represent fractures explicitly, and there is strong displacement discontinuity between fracture top and bottom surfaces. In recent years, the phase field method has been introduced to simulate hydraulic fracture propagation, which was initially applied for microstructure evolution, grain growth, and two-phase flow problems [27][28][29][30]. Without intention to review all models proposed so far, only some representative ones are mentioned. Wheeler et al. proposed an augmented-Lagrangian method for pressurized fracture problem with the phase field approach [31] and extended it to fluid-filled fractures in poroelastic media [32]. Miehe and Mauthe [33] presented a modular model to link the crack phase field evolution with the hydro-poroelastic response of the porous media for hydraulic fracturing and put forward a robust finite element implementation. Zhou et al. also conducted phase field modeling of fracture propagation in poroelastic media [34] and investigated the dynamic cracking phenomenon [35]. It has been well proved that it becomes very convenient to handle fracture initiation, branching, and joining with the phase field method in contrast with other methods for these crack behaviors are automatically captured without other additional criterions. Considering the complex interaction modes between hydraulic fracture and vugs, the phase field method is an excellent approach for the problem. Thus, a phase field model is established to simulate hydraulic fracture propagation with preset vugs in this study for the purpose of getting a more comprehensive understanding of vug effects in hydraulic fracturing.
The structure of the paper is organized as follows. Section 2 presents the mathematical model, including phase field representation of crack propagation and fluid flow equation considering crack and vug with phase field. Section 3 describes the numerical solution, including calculation of driving force of phase field, and finite element approximation as well as iterative scheme. Section 4 validates the proposed model and analyzes the effect of vug property on fracture propagation. Section 5 gives some conclusions.

Mathematical Model
2.1. Crack Representation by Phase Field Method. Following Miehe et al. [36], the sharp crack topology can be approximated by a diffusive one, which is described by an auxiliary field variable as shown in Figure 1. The field variable is characterized for d = 0 the intact state and for d = 1 the fully broken state in rock, denoted as phase field. The crack surface can be represented by a functional as follows: where ϒ ðd,∇dÞ is the crack surface density function defined as where l 0 is the length parameter which controls the spread of damage. As the length l 0 approaches 0, the sharp crack topology is recovered.
In the variational theory, the total potential energy is the sum of elastic strain energy and crack surface energy as follows: where ψðε e Þ is the elastic strain energy density and G c is the critical energy release rate of the rock. Furthermore, it needs to consider the effect of fluid pressure on the total potential energy in the porous media. There are several ways to add the pressure term [37,38]. The expression given by Lee et al. [38] is adopted as follows: where p is the fluid pressure and α is Biot's coefficient. Considering the diffusive crack topology, the crack surface energy can be rewritten in the form of domain integral with crack phase field: The involution of the crack phase field represents the damage of rock, which gives rise to stress degradation. Hence, the elastic strain energy density can be written as where σ is the total stress tensor and ε e is the elastic strain tensor. C is the stiffness tensor of the rock. k 0 is a parameter taken as small as possible while keeping the system of equations well conditioned. Substituting equations (5) and (6) into equation (4) gives the final expression of total potential energy: The variation of internal energy increment is written as follows: After some mathematical manipulations on equation (7), the variation of internal energy can be expressed as follows: The term σ ′ is the effective stress, and it is related to the total stress as follows: The term ψ 0 ðε e Þ represents the undamaged elastic strain energy density, which is given as By using the divergence theorem, the equation of energy variation is deduced to On the other hand, the variation of external work is expressed as where b is the body force and h is the boundary traction imposed on ∂Ω h . According to the variational theory, it is required that δW int = δW int holds for arbitrary values of δu and δd, which results in the strong form of governing equation for rock deformation and phase field evolution: And the corresponding boundary conditions are given as where u is the prescribed displacement on boundary ∂Ω u .

Fluid Flow in Porous Media with Vug.
Fracturing fluid is injected from the wellbore and used to initiate and extend fractures. Due to high pressure in the fracture, the fluid leaks into surrounding porous media. As the fracture propagates, the permeability of media changes. However, the fracture is not explicitly represented but by using phase field, which gives rise to the difficulty in capturing the effect of fracture on fluid flow. Moreover, the vug has also significant effect on fluid flow. Lee et al. [38] proposed a way to separate the domain into three parts including the unbroken reservoir domain, the fractured domain, and the transition domain by using phase field. Considering the existence of vug, the unbroken reservoir domain is further divided into the matrix domain and vug domain with different mechanical and seepage properties. Two thresholds parameters are given to separate the domain, denoted as c 1 and c 2 . If the phase field d is equal to or less than c 1 , the domain is considered the reservoir domain. If the phase field d is equal to or greater than c 2 , the domain is considered the fractured domain. The rest domain is the transition domain.
In the reservoir domain, the equation of mass conservation is written as where ρ is the fluid density, S r and α r are the composite compressibility coefficient and Biot's coefficient of the reservoir domain, and v r and Q r are the fluid velocity and source term in the reservoir domain. The composite compressibility coefficient S r is given by where φ is the porosity of the reservoir domain. c is the fluid compressibility, and K r is the bulk modulus of the reservoir domain.
The fluid velocity v r is related to the pressure gradient as follows: where μ is the fluid viscosity and k r is the permeability of the reservoir domain.
To capture the fluid flow in the vug, the reservoir domain is further divided into two domains, which are the matrix domain and vug domain with different mechanical and seepage properties. The vug domain is viewed as a domain with very large permeability k v (several orders of magnitude larger than matrix permeability k m ) and small modulus E v (several orders of magnitude smaller than matrix modulus).
In the fractured domain, the equation of mass conservation is written as where S f is the equivalent compressibility coefficient of the fracture domain, which is equal to the fluid compressibility. v f and Q f are the fluid velocity and source term in the fractured domain.
The fluid velocity v f in the fractured domain is given by Darcy's law as follows: where k f is the permeability of the fractured domain.
To describe fluid flow in the transition domain, two linear indicator functions are defined as presented by Lee et al. [38], which have shown good effectiveness and accuracy. The indicator functions are given as By using the indicator functions, the equation of mass conservation in the transition domain can be written as where S t and α t are the equivalent compressibility coefficient and Biot's coefficient of the transition domain and v t and Q t are the fluid velocity and source term in the transition domain.
Biot's coefficient α t is expressed as The equivalent compressibility coefficient S t in the transition domain is rewritten as where φ denotes the porosity of the transition domain, which is given as The fluid velocity v t is written as where k t is the effective permeability of the transition domain, which is expressed by So far, the governing equations for all domain parts can be generalized into the unified form (equations (25) and (29)) by using the phase field indicator functions. Through solving the unified governing equations, the fluid pressure can be obtained in all domains including the matrix domain, vug domain, fractured domain, and transition domain.

Numerical Solution
3.1. Driving Force of Phase Field. In the governing equation of crack phase field, the driving force of the damage is the elastic strain energy density. To avoid crack initiation in compression, the elastic strain energy density is divided into the tensile and compressive parts [39].
Firstly, the elastic strain is decomposed into the tensile and compressive parts as follows: where ε + e and ε − e are the tensile and compressive modes of the elastic strain. And they can be calculated as where m is the dimension number of the problem. ε i e and n i are the eigenvalue and eigenvector of the elastic strain.
The Macaulay bracket is defined as The decomposition of the elastic strain energy density can be written as where λ and ω are the Lame constants of rock. The operator tr is to access the trace of the tensor. Finally, the local history field of the maximum tensile reference energy density is adopted as the driving force of the crack phase field as follows: 3.2. Finite Element Approximation. There are three primary fields needed to be solved: displacement, fluid pressure, and phase field. The finite element method is adopted to discretize the governing equations of these fields. And then, the problem is coupled solved by a staggered iterative scheme. The weak form of the displacement field can be derived by the principle of virtual work as follows: ð where δu is the virtual displacement. The weak form of fluid pressure is given as The temporary discretization of the transient term in the above equation is approximated by using the implicit backward method. The discretized form is where n and n + 1 represent the current step and next step and Δt is the time interval between them. The superscript n + 1 is omitted for simplicity in the following derivation. By using the local history field, the weak form of phase field is written as Using the Voigt notation, the primary fields can be discretized at the element level as follows: where n e is the number of the element and N i is the shape function related to node i. N u i for the displacement vector is written as Substituting the above expressions for displacement approximation into the weak form (equation (36)) and invoking the arbitrariness of the test functions, the element stiffness matrix and the external force of an element at node i for the displacement field can be obtained as follows: Correspondingly, the element stiffness matrix and the external contribution at node i for the fluid pressure field can be obtained as follows: The element stiffness matrix and the driving force at node i for phase field can be obtained as follows: 3.3. Iterative Scheme. Miehe et al. [39] proposed a robust operator split scheme that successively updated the local history field, the phase field, and the displacement field for modeling crack propagation. When the fluid flow is considered, the scheme is extended by updating the pressure field in the time step successively. The iterative scheme is illustrated in detail as follows. To obtain the fluid pressure, displacement, and phase field for next time step n + 1, a staggered iterative scheme is presented.    Step 1. Initiate the primary fields and history field based on the values at current time step n: Step 2. Compute the pressure field by using fixed crack phase field d j : Step 3. Compute the displacement field using fixed crack phase field d j and pressure field p j+1 : Step 4. Compute the local history field using fixed displacement field u j+1 : Step 5. Compute the crack phase field using fixed history field H j+1 : Step 6. Compute the relative error of primary fields between successive iterative steps: If the relative error e r is smaller than the allowable error limit, then it goes to the next time step. Otherwise, go to Step 2 and reiterate until convergence. The allowable error limit is set to equal to 1 × 10 −3 in this study, and the staggered scheme has shown good convergence capability in previous studies.

Result and Analysis
The proposed model to simulate hydraulic fracture propagation is validated against the classic KGD model. And then, the effect of vug on fracture propagation is analyzed in detail in terms of vug radius, distance from fracture to vug, and vug modulus. The effect of the injection rate on its propagation is also investigated.

Validation of Proposed Model.
To validate the proposed model for hydraulic fracturing, the classic KGD fracture model is solved numerically, and then, the numerical solution is compared with the analytical solution [40]. The input parameters are given in Table 1. Based on numerical solution, the phase field evolutions at different times are shown in Figure 2. The phase field evolution represents fracture propagation, and it shows that the fracture propagates along the horizontal direction from two sides. Besides, the comparison between numerical solution and analytical solution is shown in Figure 3. It can be seen that the numerical solution is generally consistent with the analytical solution. Some inconsistency may take its rise from separating the domain into different parts for fluid flow in the porous media. However, the relative error is still in an allowable range. This proves that the proposed model is capable of capturing hydraulic fracture propagation in porous media.

Effect of Vug Radius.
Vugs are of different shapes and sizes. Here, the vugs are assumed to be circles with different radii, and the physical model of hydraulic fracture propagation with preset vug is illustrated in Figure 4. To investigate the effect of vug radius on fracture propagation, three cases are considered by setting the radius to 2.5 mm, 3.5 mm, and 4.5 mm. The results of phase field distribution are shown in Figure 5, in which the white circle represents the vug. It shows that the fracture firstly propagates along the horizontal direction from both ends and then diverts to the vug. This indicates that fracture propagation direction changes due to the vug, and it shows an attraction effect of vug on fracture propagation. If the vug radius is 2.5 mm, fracture reinitiates from the vug after its approaching. However, when the vug radius increases to 3.5 mm and 4.5 mm, fracture propagation mainly occurs on the other side, and it is difficult to reinitiate from the vug surface. It indicates that a large vug would hinder fracture propagation.
The pressure in the domain changes due to fluid injection and fracture propagation. And the pressure evolutions at two observation points including the injection point and the bottom point of vug circle are shown in Figure 6. At   Figure 7. It shows that when the vertical distance from fracture and vug center is less than and equal to 5 mm, the fracture diverts to vug. Otherwise, when the vertical distance between them increases to 6.25 mm, the fracture propagation direction keeps unchanged, which indicates the vug has little impact on fracture propagation. The evolutions of fluid pressure at the two observation points are shown in Figure 8. At the injection point, when the fracture diverts to the vug, the fluid pressure evolution becomes complex at a local time interval, and the fluid pressure could rise above that pressure when the vug has no effect on fracture propagation and finally descends below that pressure. However, at the bottom point of the vug, when the fracture diverts to the vug, the fluid pressure is greater than that when the diversion of fracture does not occur.

Effect of Vug Modulus.
The vug is usually filled with some minerals, so it is of different equivalent Young's modulus. To investigate the effect of vug equivalent modulus on fracture propagation, three cases are conducted by setting the vug modulus to 10 -6 , 10 -3 , and 10 times matrix modulus. The distributions of phase field are shown in Figure 9. It shows that new fracture reinitiates from the vug without offset when the vug modulus is less than that of the matrix, while the fracture propagates around the vug circumference and reinitiates from a favorable point. The intrinsic mechanism is that the vug modulus determines the deformation and failure in the vug. If the vug modulus is much smaller than the matrix modulus, the vug is very vulnerable to failure, and hydraulic fracture is prone to crossing the vug. And as the vug modulus increases, the vug develops to be an obstruction and it becomes difficult for the hydraulic fracture crossing.

Effect of Injection
Rate. At last, the effect of the injection rate on fracture propagation is analyzed by setting it to 0.24 mm 2 /s, 0.28 mm 2 /s, and 0.32 mm 2 /s. The obtained phase field is shown in Figure 11. It can be seen that as the injection rate increases, the newly initiated fracture from the vug has longer length, which indicates that the injection rate could affect the propagation velocity and length of the fracture.
The evolutions of fluid pressure at the two observation points are shown in Figure 12. It shows that as the injection rate increases, the final stabilized pressure rises slightly at both the injection point and the bottom point of the vug. The final distributions of fluid pressure are also presented as shown in Figure 13.

Conclusion
A hydromechanical coupling model based on the phase field method is presented to simulate hydraulic fracture propagation with preset vugs. The numerical model is validated by comparison with analytical solution to the classic KGD fracture model. A series of numerical cases are done to study the effect of vug on fracture propagation. Based on the simulation results, some conclusions are summarized as follows.
(1) Due to the existence of vugs, the fracture deflects from its initial propagation direction and tends to divert to the vug. When the vug radius is small, the fracture can propagate across the vug. Otherwise, the fracture may be arrested by it. As the distance between fracture and vug increases, the effect becomes weaker until it disappears, and the fracture proceeds along its initial propagation direction (2) The modulus of vug fillings determines the fracture propagation pattern. When the vug filling modulus is small, the fracture crosses the vug without offset.
However, when the vug filling modulus is large, the fracture propagates around the vug circumference and reinitiates from the offset point. The increase of the injection rate is conducive to fracture crossing the vug, which results in longer fracture length (3) The vugs have great influence on hydraulic fracture propagation. Thus, to obtain clear characterization of vug properties and geometrical location is very important, which underlies the optimization design.
Besides, there are also many natural fractures in carbonate reservoirs; the combined effects of natural fracture and vugs on hydraulic fracturing are important and will be investigated in the future work

Data Availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest
The authors declare that they have no conflicts of interest. Lithosphere