Numerical investigation of poroelastic effects during hydraulic fracturing using XFEM combined with cohesive zone model

. Hydraulic fracturing is an effective technology used to stimulate hydrocarbon production from unconventional reservoirs. Numerical simulation of the hydraulic fracturing process plays a key role in the design of hydraulic fracturing. However, most of existing models are overly simplified by neglecting the poroelasticity of rock matrix, which may have a substantial effect on the growth of hydraulic fractures, especially in oil and/or water saturated formations. In this study, we developed a fully coupled finite element model for hydraulic fracture propagation in permeable formations by combining XFEM technique with cohesive zone model, and used it to investigate the effects of poroelasticity on the geometry evolution of single fracture which initiated from an injection point. Fluid flow within the fractures, Darcy flow within the rock matrix, hydraulic fractures propagation, and fluid leak-off into the formation are simultaneously taking into account in the model. Then, the model is validated by comparing the results with available analytical solutions. To understand and quantify the poroelastic effects on the propagation of hydraulic fracture, several cases with different matrix permeability, leak-off coefficients, and bulk modulus of pore fluid are performed. The simulation results show that the total volume of leakage is controlled by the combined action of matrix permeability and leak-off coefficient. The fracture aperture and length decrease with the increase of matrix permeability or leak-off coefficient, while as the bulk moduli of pore fluid increases, the fracture aperture and length tend to increase.


Introduction
Hydraulic fracturing is an effective technology used to stimulate hydrocarbon production from unconventional reservoirs, in which one or more fractures are driven by the internal pressurized fluid (Bunger and Cardella, 2015). The growth of hydraulic fracture in porous media is a complex, multi-physics problem, which can be generally considered as the following coupled processes (Adachi et al., 2007;Detournay, 2016): (i) rock deformation caused by the fluid pressure on the fracture surfaces, (ii) fluid flow within the fracture, (iii) fracture propagation, and (iv) fluid leak out of the fracture into the formation matrix. Accurately predicting the growth behavior of hydraulic fracture is necessary to design a suitable fracturing scheme, which requires a fully coupled numerical model that can simultaneous account for the aforementioned processes. In the early times, analytical and semi-analytical methods (Detournay, 2004;Nordgren, 1972;Perkins and Kern, 1961) have been developed to understand the evolution of variables of interests, such as fracture aperture, fracture length, and injection pressure. Due to the inherent limitation of analytical methods, these solutions, however, are usually limited to some simple situations such as constant injection rate and linearly elastic medium. To simulate the growth behavior of hydraulic fracture in more sophisticated situations, a large number of numerical approaches have been proposed in literatures, including the boundary element method (Gordeliy and Detournay, 2011;Long and Xu, 2017;Zhang and Jeffrey, 2008), the distinct element method (Damjanac and Cundall, 2016;Han et al., 2015), the finite element method (Carrier and Granet, 2012;Guo et al., 2015;Li et al., 2017a;Li et al., 2017b), the finite discrete element method (FDEM) (Lisjak et al., 2017;Yan et al., 2015), the discontinuous deformation method (DDA) (Choo et al., 2016), and the extended finite element method (XFEM) (Mohammadnejad and Khoei, 2013;Salimzadeh and Khalili, 2015). Although significant progress has been made, there are still some critical points that need to be elucidated. For example, most of existing models use Carter leak-off model to characterize the velocity of fluid leak out of fracture, but neglect the poroelasticity effect of rock matrix, which may have substantial impacts on the growth of hydraulic fracture, especially in oil and/or water saturated formations. To understand and quantify the poroelasticity effect on the evolution of hydraulic fracture, this paper develops a fully coupled XFEM-based cohesive model which is able to take into account several crucial physical aspects during hydraulic fracture propagation in porous media, including fracture propagation, fluid flow within fracture, fluid leak-off, pore fluid flow, and rock deformation. In the following, the governing equations used in the model are provided first. Next, this model is validated against analytical solutions. Finally, parametric sensitivity analyses are performed to investigate the impact of poroelastic effect on the growth of hydraulic fracture.

Deformation of porous medium
The rock formation is assumed to be a homogeneous, isotropic, and linear poroelastic medium. For quasi-static conditions, the stress equilibrium equation can be written as where F is the body force vector, and σ is the total stress tensor. Assuming that the rock matrix is fully saturated by single-phase fluid and obeys poroelasticity theory, the effective stress of rock matrix reads (Detournay and Cheng, 1993) p     σ σ I (2) where  σ is the effective stress, p is the pore pressure, I is the second-order unit matrix, and  is the Biot coefficient, which is defined as where K and s K are the bulk modulus of porous rock and of the rock matrix grain, respectively. The effective stress is related to the matrix strain as   σ Dε (4) where D is the elastic stiffness matrix, and ε is the strain tensor. The pore fluid flow obeys mass balance equation as where f  is pore fluid density,  is the rock porosity,  (6) in which s v is the solid phase (rock matrix) velocity, and d v is the Darcy's fluid velocity, can be defined as (7) in which m k is the permeability of the rock matrix, f  is the fluid viscosity, and g is the vector of gravitational acceleration.

Fracture fluid flow
Assuming that the hydraulic fracture is filled with a Newtonian fluid, the mass balance equation for the fluid flow within the fracture reads (Hibbitt et al., 2016) where is the fracture width, is the tangential fluid flow velocity within the fracture, � and � are the fluid leak-off into the formation across top and bottom fracture surfaces, respectively. The tangential fluid flow velocity obeys cubic law as in which � is the fluid pressure within the fracture. The normal fluid velocity per unit area ( � and � ) can be described as where � and � are the pore fluid pressures in the rock matrix abut to the top and bottom surfaces, respectively. � and � are the "leak-off coefficients" of top and bottom surfaces, respectively. Different with the classical Carter's model, � and � here characterize the resistance of fluid across the fracture surfaces. The leak-off velocity is controlled by the rock matrix permeability if � and � are large enough.

Fracture propagation
XFEM technique, implemented in Abaqus, is used in this paper to simulate the hydraulic fracture propagation in porous media. The fracture propagation behavior is described by cohesive zone model with linear tractionseparation response, as shown in Fig. 1. Damage is assumed to initiate when a quadratic interaction function involving nominal stress ratios reaches 1.0, which can be expressed as (Hibbitt et al., 2016) where � , � , and � are the nominal traction in the normal and two local shear directions, respectively; � � , � � , and � � represent the cohesive strength in the normal and two local shear directions, respectively.
where � � ⁄ � ⁄ are the stresses in the corresponding directions predicted by linear traction-separation law without damage. is the scalar damage parameter.

Model validation
In this section, Terzaghi's problem and KGD problem are performed to verify the capability of the XFEM-based cohesive zone model on modeling the process of fluid diffusion and the growth of hydraulic fracture in rock matrix, respectively.

Terzaghi's problem
Consider a water-saturated, plane strain porous column, with the height of 100 m, as shown in . Fig. 2(a). The top surface is drained and other surfaces are undrained boundary conditions. The normal displacement of left, right, and bottom boundaries are fixed. A load � acts instantaneously on the top boundary at initial time. The relevant parameters of this case are shown in Table. 1 and the analytical solution are provided in Appendix A.

KGD problem
As shown in Fig. 3(a), a bi-wing fracture is created by injecting a Newtonian fluid into an impermeable medium. Because the typical growth of hydraulic fracture in an unconventional reservoir is under viscosity dominated, the parameters used in this case were tuned to this regime, as provided in Table. 2. The details of semi-analytical solutions are shown in Appendix B. Fig. 3 (b) and (c) show the temporal evolution of injection pressure and mouth aperture. The agreement between numerical results and analytical solutions is very good, which proves the adequacy of the XFEM-based cohesive zone model on modeling the growth of hydraulic fracture.

Poroelastic effect on hydraulic fracture propagation
This section aims to investigate the impact of poroelastic effect on the propagation behavior of hydraulic fracture.
To this end, a series of parametric analyses for KGD problem are performed in this section. As shown in Fig. 4, a 2D plane-strain fracture is propagated by injecting a viscos fluid from the initial crack. The dimension of the domain is 200×50 m, and the length of the initial crack is 0.4 m. The left boundary of the model is symmetric boundary condition, the normal displacements of other external boundaries are fixed, and a constant pore pressure boundary condition is also imposed on these boundaries. The input parameters for the base case of parametric analysis are listed in Table. 3. Additionally, in each parametric analysis, only the investigated parameter is changed while all other parameters are kept same with the values from the base case.

Effect of permeability
As an important parameter for porous media properties, permeability is expected to significantly influence the pore pressure distribution around the hydraulic fracture, and consequently affects the fracture evolution. Thus, we perform a sensitivity analysis by varying permeability in this section. Assume there is no filter cake developed on the fracture surfaces. Fig. 7 depicts the fracture profile and pore pressure distribution at the end of the injection for different rock permeability. It readily observed that fracture length increases monotonically as permeability decreases. Large permeability leads to high fluid loss volume, which indicates that the leak-off in a system with large leak-off coefficient is controlled by rock permeability. Fig. 7 Effect of permeability on the fracture profile and pore pressure distribution Fig. 8 shows the effect of permeability on the net pressure and aperture at the injection point as time elapse. It can be observed that the net pressure of high permeability case is much larger than the low permeability case, while the fracture aperture of high permeability case is much lower than the low permeability case. This is due to the large back-stress induced by poroelastic effect in high permeability case, which in turn reduces the aperture.

Effect of leak-off coefficient
Leak-off coefficient in this work characterizes the resistance of filter-cake during the fracture propagation. Fig. 9 shows that the impact of leak-off coefficient on the growth of hydraulic fracture is substantial. Fluid loss is controlled by permeability in base case while it is controlled by the leak-off coefficient in Fig. 9(d) and consequently more fluid is utilized to drive the fracture propagation. Fig. 9 Effect of leak-off coefficient on the fracture profile and pore pressure Fig. 10 shows the temporal evolution of net pressure and aperture at the injection point for different leak-off coefficient. The results indicate that lower leak-off coefficient results in lower propagation pressure and higher fracture aperture. This is because the back-stress induced by poroelastic effect is lower.

Effect of bulk modulus of pore fluid
Bulk moduli of pore fluid is one of factors determining the dissipation rate of pore pressure, which may influence the hydraulic fracture propagation. However, to the authors' knowledge, there is no previous studies have investigated its influence on the fracture propagation. Fig.  11 shows the fracture profiles for different bulk modulus of pore fluid. The results indicate that low bulk moduli of pore fluid tends to suppress the increase of the pore pressure around hydraulic fracture, which leads to high leak-off velocity and low fracturing fluid efficiency.

Conclusions
In this paper, a fully coupled finite element model is proposed to investigate the influences of poroelastic effect on the growth of hydraulic fracture, and the model was validated against available solutions. The simulation results show that the total volume of leakage is controlled by the combined action of matrix permeability and leakoff coefficient. The fracture aperture and length decrease with the increase of matrix permeability or leak-off coefficient, while as the bulk moduli of pore fluid increases, the fracture aperture and length tend to increase.