An XFEM Model for Hydraulic Fracturing in Partially Saturated Rocks

Hydraulic fracturing is a complex multi-physics phenomenon. Numerous analytical and numerical models of hydraulic fracturing processes have been proposed. Analytical solutions commonly are able to model the growth of a single hydraulic fracture into an initially intact, homogeneous rock mass. Numerical models are able to analyse complex problems such as multiple hydraulic fractures and fracturing in heterogeneous media. However, majority of available models are restricted to single-phase flow through fracture and permeable porous rock. This is not compatible with actual field conditions where the injected fluid does not have similar properties as the host fluid. In this study we present a fully coupled hydro-poroelastic model which incorporates two fluids i.e. fracturing fluid and host fluid. Flow through fracture is defined based on lubrication assumption, while flow through matrix is defined as Darcy flow. The fracture discontinuity in the mechanical model is captured using eXtended Finite Element Method (XFEM) while the fracture propagation criterion is defined through cohesive fracture model. The discontinuous matrix fluid velocity across fracture is modelled using leak-off loading which couples fracture flow and matrix flow. The proposed model has been discretised using standard Galerkin method, implemented in Matlab and verified against several published solutions. Multiple hydraulic fracturing simulations are performed to show the model robustness and to illustrate how problem parameters such as injection rate and rock permeability affect the hydraulic fracturing variables i.e. injection pressure, fracture aperture and fracture length. The results show the impact of partial saturation on leak-off and the fact that single-phase models may underestimate the leak-off.


Introduction
Flow through fractures in deforming porous media has been the subject of great interest in many engineering disciplines.Typical examples are: fractured oil reservoirs [1], hydraulic fracturing for enhanced hydrocarbon production, tight gas reservoirs [2], weakly consolidated offshore sediments [3], soft coal bed methane extraction [4], geothermal energy [5,6], isolation of hazardous waste [7], measurement of in situ stresses [8], fault reactivation [9], and remediation of soil and ground water aquifers [10] to name a few.
Hydraulic fracturing is a complex multi-physics problem.A fluid is injected into the host rock; the fracture initiates and propagates due to the induced hydraulic loading.The simulation of the process requires several components: (i) flow of the fracturing fluid within the fracture; (ii) the mechanical deformation of the porous medium due to the applied hydraulic load; (iii) the fracture propagation; (iv) the leak-off flow to the surrounding porous media (host rock); and (v) the fluid flow through the host medium.The first three components are widely addressed in the literature, while last two are considered only in few studies albeit with some restrictions.
Based on the energy-dissipation mechanism and the ability of the rock matrix to dissipate fracturing fluid, four combined asymptotic regimes are defined: storageviscosity, storage-toughness, leak-off-viscosity and leakoff-toughness [11].Asymptotic solutions provide fundamental understanding of the hydraulic fracturing process, and provide benchmarking cornerstones for numerical models.However, analytical solutions are restricted to simplified fracture geometries in homogenous media and are typically constrained to a set of fixed boundary conditions.Standard geometries are 2D plane-strain fractures (PKN fracture [12,13]; KGD fracture [14][15][16], and radial (penny-shaped) fractures [17][18][19].
In majority of the available models, the flow through rock matrix and the fluid exchange between fracture and rock matrix, is either ignored by assuming impermeable rock formations, e.g.[20], or simplified using onedimensional analytical leak-off model, e.g.[21].The impermeable matrix assumption is proven to be an unrealistic assumption by the substantial evidence of conductive rock matrix [1,4].One-dimensional leak-off model, on the other hand, is based on Carter's onedimensional leak-off model [22], which is represented as a sink term in the mass balance equation of fracture flow.This approach has several shortcomings including onedimensional assumption of flow, time-dependency of flow instead of pressure-dependency and more importantly ignoring the poroelastic effects.Flow dissipation through rock matrix is multi-dimensional unless the permeability in the direction normal to the fracture plane is significantly higher than in other directions [23].Also, the fluid exchange between fracture and matrix depends on the pressure gradient, while in Carter's leak-off model the fluid leak-off rate at a given location within the fracture depends only on the fracturing fluid pressure history at this position and does not depend on the pressure in the adjacent regions.As time elapses, the leak-off rate at each position decreases proportionally to 1/√t, where t is time; therefore, a scenario of fracture arrest is not possible [24].Finally, the seepage of the fracturing fluid into the rock formation increases fluid pressure in the matrix causing dilation of rock matrix.Dilated matrix applies stresses back on the fracture, known as 'back-stress' in hydraulic fracturing context, tending to close the fracture [19,25,26].These shortcomings have also affected the available semianalytical solutions for leak-off-dominated regimes that use simplified one-dimensional leak-off model in deriving solution equations.These solutions fail to predict hydraulic fracturing parameters in leak-offdominated regimes accurately as shown in [23], and in [27] for single-phase flow, and in [24] for two-phase flow in two-dimensional case.
In this paper, a fully coupled extended finite element model for hydraulic fracturing in porous media is proposed.Two fluids: fracturing fluid and host fluid are considered.The independent fluid flow for the porous host medium makes the present model capable of simulating multidimensional leak-off fluid flow through high/low permeable, homogeneous/inhomogeneous host medium with transient pressure boundary condition at the fracture-rock boundary.The fracture discontinuity in the deformation model is handled using XFEM technique, while cohesive crack model is used as fracture propagation criteria.The model is verified against several benchmark solutions from the literature.It is shown that leak-off is higher when partial-saturation of rock matrix is considered.

Governing equations
The model proposed comprises five phases and three models: rock skeleton (deformation model), wetting and non-wetting fluids in fracture (fracture flow model), and wetting and non-wetting fluids in matrix (matrix flow model).The differential equation describing poroelastic deformation of rock matrix can be written as [24] div(1/2 D( du+du )-β w dp pw I-β nw dp pnw I)+dF=0 (1) where u donates the displacement vector of the porous medium, β is the incremental effective stress parameter, p is the fluid pressure, I is the second-order identity tensor, and F is the body force per unit volume.
Following the cohesive fracture mechanics, a tractionseparation law governs the nonlinear behaviour of the fracturing medium in the cohesive zone.In addition to cohesive traction, in fluid saturated fracture discontinuities, the hydraulic loading is applied to the fracture planes.The hydraulic loadings can be either positive for pressurised fluid within the fracture or negative due to suction in the lag region.Therefore, the natural boundary conditions for loading at the fracture discontinuity is expressed as where σ is the total stress, t is the cohesive tractions, p f is the average fluid pressure within the fracture, and n c is the outward unit vector normal to the discontinuity (Figure 1).Notice that in the above equation Biot's coefficient is not applied to fracture loading, and that the fluid pressure in the fracture is treated identical to an external pressure applied to fracture planes.Such a pressure does not require scaling according to Biot's coefficient.
The flow through planar fracture is commonly modelled using lubrication theory [28] for an incompressible Newtonian fluid obeying cubic law.Lubrication equation is derivable from general Navier-Stoke's equation for the flow between two parallel plates [29].The fracturing fluid can have properties substantially different from those of the host fluid.Furthermore, depending on the in-situ stresses, the fracturing fluid front and the fracture front may not coincide, giving rise to the so-called fluid lag phenomenon [1,30,31].In the lag region, the host fluid enters the fracture under the developed suction.Therefore, a second fluid is required to complete the fracture flow model.
Independent fracture flow model is considered for fracture discontinuity in this research.Fracture pressures are calculated directly, which may be linked to multiple pore pressures within the porous medium.The objective is to obtain a more realistic representation of fracture flow compared to current models presented in the literature in which gradient of fracture flow is calculated entirely based on the enriched component of the fluid pressures [32,33,34].The extended lubrication equation considering compressibility of fluid constituent and fluid transfer between fracture and matrix (leak-off) can be written as [24] ∂/∂l [S fξ k frξ w 3 /(12μ ξ ) ∂p fξ /∂l] = S fξ wc ξ ∂p fξ /∂t + S fξ ∂w/∂t where w is the fracture aperture, μ is the fluid dynamic viscosity, S fξ is the degree of saturation of the fluid ξ, k fr is the relative permeability, c ξ is the fluid compressibility, p pξ is the fluid pressure in the adjacent matrix, , y is the distance between fracture and adjacent matrix, and l is the dimension along the fracture.The change in saturation is linked to the change in capillary pressure (suction) through available capillary pressuresuction curves.In this study the relationships proposed by van Genuchten-Mualem (VGM) model are utilised for capillary pressure-saturation and for relative permeability-saturation interpolations.Details are given in [24].Note that fracture pressure loading on the fracture walls (Eq.2) and fracture aperture in Eq. 3 provide where, n is the porosity and s p is the capillary pressure(suction) in the porous medium.The pressure field is assumed continuous across the fracture discontinuity, with the discontinuity in the Darcy flow velocity normal to the fracture taken into account through leak-off loading.Note that leak-off terms in Eqs. 3 and 4 provide symmetric coupling between two flow models.

Schematic representation of the problem with the
body Ω, boundary Γ, fracture boundary Γ c .

Numerical approximation
Spatial discretisation has been performed using the standard Galerkin method with displacements and fluids pressures defined as primary variables.Extended finite element method (XFEM) is used to model the fracture discontinuity within displacement field at element level.Arbitrary discontinuities can be handled using XFEM without the need for remeshing.In XFEM, the displacement field consists of two parts, one continuous and the other discontinuous.The continuous part is the standard displacement field corresponding to the situation without any cracks.The discontinuous displacement field (also known as enriched part) is based on local partitions of unity and allows the element to include a discontinuity.
Heaviside step function, also known as 'jump' function, is typically used to enrich the nodes in fully cracked elements if their support was cut by the crack into two disjoint pieces [35].To obtain numerical solution of the governing equations, the rate form of the discretised equations is integrated over the time domain.The time integration is performed using a generalized trapezoidal method.
The set of discretised equations are highly nonlinear.Nonlinearities arising from fluid pressure dependency of constitutive coefficients are handled using Picard method such that the coefficient matrices appearing in the stiffness matrix as nonlinear functions of unknown variables are updated in every iteration.The four-node quadrilateral elements together with two-point Gaussian quadrature rule are used for numerical approximation and numerical integration.Numerical integration of the elements bisected by the fracture discontinuity is done using common method of element partitioning in which numerical integration is done separately within each subdomain.The model has been implemented in a wellknown code (Matlab) and verified against several analytical solutions and published results [27,36].

Hydraulic fracturing in partially saturated medium
In these examples, a viscous fracturing fluid is injected into a fracture with initial length of 1 m.The host medium is assumed to be saturated with gas as nonwetting fluid at atmospheric pressure.The model dimensions and mesh configuration is shown in Figure 2. The model parameters used are: Young's modulus E=50 GPa, Poisson ratio υ=0.25, tensile strength f t =1000 kPa, and fracture energy G f =0.1 KN/m.The fracture is assumed to have plane-strain behaviour (KGD fracture) with fracture height of 100 m.
The simulated examples in this section include: i) the case with impermeable rock matrix (k p =0) to validate current model against available semi-analytical solutions, ii) the case with low matrix permeability (k p =1×10 -13 m 2 ) and iii) the case with high matrix permeability (k p =1×10 - 12 m 2 ).The medium porosity in all examples is n=0.01.Geertsma and de Klerk [14] obtained an approximate analytical solution for the case without leak-off and with negligible toughness (KGD model).Later, Spence and Sharp [15] extended the solution to include fracture toughness.
The results from present study for the net well pressure, the fracture aperture and the fracture half-length are shown in Figures 3 to 5, respectively.Included in these figures are the analytical solutions by Geertsma and de Klerk [14] and Spence and Sharp [15].In Figure 5 the solution by Carbonell et al. [37] for the limiting case of storage-viscosity regime and the solution by Adachi and E 2016 -Detournay [16] for leak-off dominated regime are also included.The leak-off dominated solution is calculated for the most permeable case with k p =1×10 -12 m 2 which represents the most suitable case for this solution.
Model used for hydraulic fracturing simulations.
For the impermeable case very good agreement is found between present model results and KGD analytical solutions for net well pressure and fracture aperture at the wellbore.KGD solution is for viscosity dominated regime (zero toughness), and the parameters used in the present simulation falls in to the viscosity-dominated regime.This is verified through computing the dimensionless viscosity and toughness: M=308 and K=0.239 [38].Viscosity-dominated regime is in accordance with actual hydraulic fracturing practice, where the energy dissipation associated with the flow of viscous fracturing fluid often dominates over the energy dissipation associated with the rock damage [19].Spence and Sharp [15] calculated higher net well pressure due to including finite toughness.The fracture half-length (distance from fracture origin to the fracture tip) is larger than the analytical solution due to the pressure lag ahead of the fracturing fluid front as shown in Figure 6.The size of the lag region can be affected by the in situ stresses [11,30].The in situ stresses are considered zero in these simulations to allow the formation of lag zone.Negative pressure (suction) is induced in the lag region.If the position of the fracturing fluid front is considered instead of the fracture tip, as shown in Figure 5, a perfect match is observed between present model and KGD results.Again, Spence and Sharp [15] calculated slightly lower fracture length, while Carbonell et al. [37] calculated higher fracture length for the limiting case of viscosity-storage regime.
When the leak-off is allowed by considering permeable host medium, a significant portion of the injected fracturing fluid is dissipated into the surrounding media.This causes smaller fracture aperture and lower fracture propagation.The net well pressure on the other hand remains higher and maintains a steady leak-off flow towards rock matrix which becomes equal to the injected flow rate as time elapses.The higher the medium permeability, the sooner the equilibrium.For the case of higher medium permeability, the equilibrium is reached within five minutes, consequently, the fracture aperture and half-length are stabilised and no fracture propagation is observed.This implies that all the injected fracturing fluid are dissipated into the host rock.The model results for the high permeable case are in good agreement with the analytical solution by Adachi and Detournay [16].The solution is for leak-off-dominated-regime based on the 1D Carter's leak-off model.However, in late time response the two models diverge as the current model approaches an equilibrium in which there is no further fracture growth and the entire injected fracturing fluid is dissipated into the host medium, while the analytical solution predicts fracture growth proportional to square root of time.Note that the Carter's leak-off model is for single-phase flow.
Well pressure versus elapsed time.
Fracture aperture at the well versus elapsed time.The invaded zone showed by contours of degree of saturation of hydraulic fracturing fluid within host medium, top: k p =1×10 -12 m 2 , bottom: k p =1×10 -13 m 2 .
The shape of the invaded zone, i.e. the zone within host medium saturated by fracturing fluid, also depends on the permeability of the host medium.Figure 7 shows the contours of the degree of saturation of fracturing fluid within host medium for the two examples with permeable rock matrices.For the case with higher permeability, lower pressure gradient is required to push fracturing fluid into the porous medium so the invaded zone shape is short and dispersed, whereas for the case with the lower permeability it is elongated and thin.In the latter case, higher pressure-gradient is required to push fracturing fluid into the porous medium, so the higher fracture pressure is developed which causes longer fracture propagation and the preferred path for the fracture flow is along the fracture.

Conclusions
A fully coupled XFEM model for hydraulic fracturing in partially saturated medium is presented.Five phases are identified: solid skeleton, wetting and non-wetting fluids in fracture, and wetting and non-wetting fluids in rock matrix.Simulation results show that fluid dissipation is higher than what is predicted by asymptotic solutions.One-dimensional analytical leak-off model cannot accurately capture multidimensional leak-off.Furthermore, leak-off is higher when the matrix partial saturation is accounted for.