## Abstract

In this paper, we present a two-dimensional numerical model for modelling of hydraulic fracturing in anisotropic media. The numerical model is based on extended finite element method. The saturated porous medium is modelled using Biot’s theory of poroelasticity. An enhanced local pressure model is used for modelling the pressure within the fracture, taking into account the external fluid injection and the leak-off. Directional dependence of all the rock properties, both elastic and flow related, is taken into account. A combination of the Tsai–Hill failure criterion and Camacho–Ortiz propagation criterion is proposed to determine the fracture propagation.

#### Authors

#### V. Valliappan^{1}, J. J. C. Remmers^{1}, A. Barnhoorn^{2}, D. M. J. Smeulders^{1}

^{1}Department of Mechanical Engineering, Eindhoven, University of Technology, Eindhoven, The Netherlands. ^{2}Faculty of Civil Engineering and Geosciences, Delft, University of Technology, Delft, The Netherlands

Received: 3 April 2017 / Accepted: 16 November 2017

© The Author(s) 2017

We study the impact on fracture propagation (in both magnitude and direction) due to anisotropies induced by various parameters, namely ultimate tensile strength, Young’s modulus, permeability and overburden pressure. The influence of several combinations of all these anisotropies along with different grain orientations and initial fracture directions on the fracture propagation direction is studied. Different regimes are identified where the fracture propagation direction is controlled by the degree of material anisotropy instead of the stress anisotropy.

## 1 Introduction

Hydraulic fracturing is a process of inducing fractures in rock structures by injecting fluid at high pressures. Interest in hydraulic fracturing has been increasing in recent years due to its applications in the oil and gas industry for enhanced oil recovery in conventional reservoirs, for heat recovery from geothermal reservoirs and for the profitable extraction of oil and gas from unconventional reservoirs. A better understanding of the fracture growth phenomena will enhance productivity and also reduce the environmental footprint as less fractures can be created in a much more efficient way.

Several models have been developed for this purpose, the earliest of which was by Khristianovic and Zheltov (1955) and later extended by Geertsma and Klerk (1969) to form the KGD model, which proposed an analytical solution to the problem by assuming plane strain conditions. Although this model has been accepted as a standard case for hydraulic fracturing due to its simplicity, it suffers from the assumption that there is no leak-off from fracture into formation and also the formation is assumed to be solid. Another analytical solution with a different geometrical assumption (fracture length ≫ fracture height) was given by Perkins and Kern (1961) and extended to include fluid loss by Nordgren (1972). Analytical asymptotic solutions were formulated (Adachi and Detournay 2002; Garagash and Detournay 2005; Detournay 2016), based on a parametric space to identify the significant parameters in a given regime. Later, these asymptotic solutions have been modified to take into account the leak-off, fracture toughness and fluid viscosity (Adachi and Detournay 2008; Kovalyshen 2010; Dontsov 2017).

The first numerical models for modelling hydraulic fracturing were developed by Boone and Ingraffea (1990) using finite elements for modelling the formation and using finite volume for modelling flow with cohesive zones along element edges describing the fractures. In recent times there have been several models developed based on different numerical techniques. A linear elastic fracture mechanics (LEFM)-based finite element model (FEM) was proposed by Hossain and Rahman (2008). To avoid the singularity problems at the crack tip in LEFM, FEM models with zero-thickness elements (describing the fractures using cohesive zones) were developed (Carrier and Granet 2012; Chen 2012).

The FEM approach requires re-meshing to capture the fracture propagation accurately, whereas extended finite element (XFEM)-based models allows for fracture propagation in arbitrary directions without the need for re-meshing (Mohammadnejad and Khoei 2013; Remij et al. 2015; Meschke and Leonhart 2015). A novel approach in which the asymptotic behaviour near the fracture tip was resolved with extended finite element method (Gordeliy and Peirce 2013a, b, 2015). An alternate approach based on phase field modelling which combines FEM with continuum damage mechanics has been developed (Mikelic et al. 2015; Miehe and Mauthe 2016) which provides a convenient way for modelling complex fracture interactions. But all the above hydraulic fracture models assume the rock formation to be isotropic in nature.

Most rocks (especially shales, which are the most common rock type to be hydraulically fractured) are highly anisotropic in nature (Jaeger et al. 2009; Barton 2007). Kaarsberg (1959) and Sayers (1994) observed that shales have a bedding plane along which the grains are oriented causing the properties along the grain direction to be vastly different from the properties perpendicular to the grain direction. This causes a special type of anisotropy, called transverse isotropy, where the material properties in any direction in the plane can be obtained by using the material properties along any two mutually perpendicular set of directions in that plane. Although there are several studies (Abousleiman et al. 2008; Zhubayev et al. 2015) experimentally obtaining the anisotropic parameters, there are few papers by Cheng (1997) analytically deriving the anisotropic poroelastic coefficients. Abousleiman et al. (1996) modelled the deformation and pressure in a transversely isotropic porous medium without any fracture. Porous material with a stress-driven fracture in an orthotropic medium was modelled by Remij et al. (2015). More recently, the influence of rock anisotropy on tensile fractures was studied experimentally by Mighani et al. (2016).

In this paper, we enhance the aforementioned XFEM model by Remij et al. (2015) to include the effects of anisotropy on hydraulic fracturing. Using the model we analyse the effect of anisotropic rock properties (Young’s modulus, ultimate tensile strength and permeability) on the fluid-driven fracture propagation and also the impact when combined with loading.

## 2 Mathematical Formulation

In order to solve a poroelastic problem, we need to solve for the solid deformation (u) and fluid pressure (p) at all points within the porous media. Biot’s theory of poroelasticity (Biot 1941) is used to describe the porous media, and fluid flow is described using Darcy’s law. The porous medium is assumed to be saturated.

In addition for the hydraulic fracture problem, we make use of an enhanced local pressure (ELP) as proposed by Remij et al. (2015) to model the pressure inside the fracture at all points along the fracture length. This additional degree of freedom (p_{d}) enables us to model the complicated phenomenon happening within the fracture, namely the injection of external fluid, moving boundaries of fracture surface and the leak-off. Leak-off from fracture to formation is described using the 1-D Terzaghi’s consolidation equation.

For solving the unknowns, a set of governing equations along with auxiliary equations are used. The governing equations used to describe the poroelastic problem are of two types: solid deformation-based momentum balance and fluid flow-based mass balance. We consider an additional equation to ensure the mass balance inside the fracture. The auxiliary equations are used for relating these governing equations with the unknowns and also for coupling them. A schematic flow chart of the mathematical formulation is represented in Fig. 1.

**Fig. 1** Solution procedure for poroelastic fracture problem

A schematic of a body Ω with a discontinuity Γ_{d}, which splits the body into two domains Ω^{+} and Ω^{−}, along with prescribed boundary conditions is represented in Fig. 2.

**Fig. 2** Schematic of a domain Ω with discontinuity Γd

## 3 Implementation

### 3.1 Discretisation

We need a numerical method in order to solve the set of coupled differential equations and also incorporate the discontinuous jump in various parameters due to the fracture. Hence, we make use of XFEM which models the discontinuous jump due to fractures by using additional degrees of freedom. XFEM allows for the fracture to propagate through the elements, thus ensuring accuracy in capturing the fracture even with a coarse mesh making it computationally very efficient. Similar to traditional FEM, the unknowns are obtained at certain control points (nodes) by solving the differential equation and the unknowns in the intermediate regions are obtained by interpolation using shape functions. The discretised form of the unknown displacement and pressure in a porous medium intersected by the fracture is expressed as:

where H_{Γd} is the Heaviside step function given as:

u^{^} and p^{^} represent the regular nodal degrees of freedom for displacement and pressure, respectively, while u^{~} and p^{~} represent the enhanced nodal degrees of freedom representing the discontinuous jump in displacement and pressure across the fracture. p^{^}_{d} represents the pressure inside the fracture at points where the fracture intersects the element edges (Fig. 3). N, L are two-dimensional interpolation or shape functions for the displacement and pressure fields, whereas V is a one-dimensional shape function for interpolation of the pressure along the fracture length.

**Fig. 3** Discretisation in XFEM

### 3.2 Solution

The governing equations are combined with the auxiliary equations as shown in Fig. 1. Weak form of this set of equations is obtained by integrating them along with a test function. By substituting the discretised unknowns given by Eqs. (1), (2) and (3) into the weak form, we convert the set of differential equations into a set of algebraic equations. In order to solve this set of equations simultaneously, we make use of the Newton–Raphson iterative solver in combination with Euler’s forward scheme for obtaining the time derivative, and Euler’s implicit scheme for time-independent parameters. A detailed description of the solution procedure is given in Remij et al. (2015).

The unknown

degrees of freedom are solved at each grid point (nodes) for every time step.

### 3.3 Propagation

A propagation criterion is needed in order to determine the propagation initiation and also the magnitude and direction of propagation. Hence, we make use of the Camacho–Ortiz criterion (Camacho and Ortiz 1996) which can be used for mixed mode fractures as well. This criterion states the condition for fracture propagation as:

where t_{eq} is the equivalent traction ratio in a specific direction, τ_{ult} and S_{ult} represent the ultimate tensile and shear strength of the porous media, respectively, t_{n} and t_{s} are the normal and shear tractions along that orientation which are obtained from average stresses as:

where Ψ is the coefficient of friction, and n and s are the unit normal and tangent vector to the direction. σ_{av }is the average stress, used to obtain a better approximation of the stress state in the vicinity of the fracture tip. The average stresses are obtained by assigning weight functions to the Gaussian integration points within a certain distance (generally three times the characteristic element length) from the fracture tip, as derived by Jirasek (1998). As an averaged stress is used, this may lead to a slight delay in the onset of fracture propagation.

The direction of propagation is taken to be the direction in which the equivalent traction is maximum. The fracture is assumed to propagate through the entire element length in a single time step in a straight line. Further details on the implementation of the solution are described by Remmers (2006) and Remmers et al. (2003).

## 4 Anisotropic Parameters

In this section, we highlight the parameters which need to be modified to incorporate the effect of anisotropy.

### 4.1 Constitutive Relation

The effective elastic stress (σ_{e}) in the solid grain is related to the elastic strain by means of the generalised Hooke’s law.

where C is the constitutive relationship matrix and ϵ is the elastic strain in the solid grains.

The coefficients of the constitutive matrix depend on the material type and geometrical assumptions. The constitutive relation matrix can be obtained as the inverse of the compliance matrix, C=S^{−1}. We assume a plane strain case with a transverse isotropic material for which the compliance matrix in the grain direction is given as:

where E_{∥} and E_{⊥} are the Young’s moduli parallel and perpendicular to the grain direction and ν_{in} is the in-plane Poisson’s ratio, representing compressive strain perpendicular to the grain direction due to a tensile stress parallel to the grain direction and ν_{out} represents the out-of-plane Poisson’s ratio.

The constitutive matrix at any arbitrary direction is obtained from the following expression:

where T is the transformation matrix given as a function of the angle (β) between the global direction and the grain orientation direction.

### 4.2 Ultimate Tensile Strength

The ultimate tensile strength (τ_{ult}) is an important parameter which determines the fracture propagation. The ultimate tensile strength is maximum in the grain orientation direction and minimum perpendicular to it. Some previous studies (Remij et al. 2015; Yu et al. 2002; Lee and Pietruszczak 2015) have assumed cosine functions to interpolate the ultimate tensile strength at arbitrary directions from the values along the grain direction and the perpendicular direction. Here we make use of the Tsai–Hill failure criterion (Jones 1998) to model the directional dependence of the ultimate tensile strength. The Tsai–Hill failure criterion for a 2-D transverse isotropic material is given as:

where σ_{||},τ_{ult||} and σ_{⊥}, τ_{ult}_{⊥} are the stresses, ultimate tensile strengths parallel and perpendicular to the grain direction, respectively, and σ_{s} and S_{ult }indicate the shear stress and ultimate shear strength in the plane.

Assuming a fracture propagating at an angle ϕ with respect to the x-axis, the ultimate tensile strength perpendicular to the fracture propagation direction τ_{frac} is given as the normal stress perpendicular to fracture propagation direction which can satisfy the Tsai–Hill failure criterion. To obtain τ_{frac}, this stress state (σ_{frac}=[0τ_{frac}0]^{T}) is rotated to the grain orientation direction by using the stress transformation relations and then substituted into Equation (12).

where θ=γ−ϕ,γ is the grain orientation angle with respect to the x-axis, and ϕ is the fracture propagation angle with respect to the x-axis (Fig. 4).

**Fig. 4** Schematic for rotation of parameters

Hence, for a randomly oriented fracture, we obtain the ultimate tensile strength in the direction perpendicular to the fracture propagation direction (ϕ+90^{∘}) as:

In order to validate the proposed Tsai–Hill-based theory, we make a comparison with experimental results. Mighani et al. (2016) conducted tensile fracture experiments on Lyon’s sandstone and pyrophyllite rocks and observed the ultimate tensile strength variation with θ which is the angle between the grain orientation direction and the fracture propagation direction. We make use of the experimental values for maximum and minimum ultimate tensile strengths and interpolate for various angles. As one can observe in Fig. 5, Tsai–Hill failure criterion provides a much better fit for the experimental values when compared to the previously assumed cosine functions (Remij et al. 2015; Lee and Pietruszczak 2015).

**Fig. 5** Comparison of analytical and experimental values of ultimate tensile strength at various angles with respect to grain orientation direction

### 4.3 Poroelastic Coefficients

Cheng (1997) derived the analytical expressions for the transverse isotropic poroelastic coefficients based on the constitutive relationship matrix.

In a fully saturated porous medium, the external stresses on the porous media are partly taken by the fluid pressure in the pores and partly by deformation of the solid grains. This can be represented mathematically as:

where σ is the total stress, σ_{e} is Terzaghi’s effective stress, and α is Biot’s coefficient matrix.

The anisotropic Biot’s coefficient matrix reduced to the 2-D form is given by the expression:

The other poroelastic constants such as bulk modulus and compressibility modulus are given as:

where nf is the porosity of the porous media, and K_{s} and K_{f} are the bulk modulus of the solid and fluid, respectively.

### 4.4 Permeability

Permeability enters into the formulation through Darcy’s law which describes the fluid flow in the porous medium as:where κ is the intrinsic permeability tensor, μ is the dynamic viscosity of the pore fluid, q is the flux, and p refers to the pressure in the porous media. Permeability tensor for an anisotropic rock in the principal grain directions is given aswhere permeability in the grain direction κ_{∥} is larger than the permeability perpendicular to the grain direction κ_{⊥}. For permeability in any arbitrary orientation we make use of the transformation matrix:

5 Results

### 5.1 Validation

Since there are no studies which exactly deal with hydraulic fracturing in anisotropic media, we divide the validation into two parts: (1) Mandel’s problem which compares the numerical results with an analytical solution for a transverse isotropic porous medium without fractures and (2) the standard KGD problem which compares the numerical results with the ELP model (Remij et al. 2015) for a hydraulic fracture problem in an isotropic medium.

#### 5.1.1 Mandel’s Problem

**Fig. 6** Geometry and boundary conditions of Mandel’s problem

Abousleiman et al. (1996) provided an analytical solution for Mandel’s problem in a transversely isotropic porous medium. Mandel’s problem (Fig. 6) consists of an infinitely long rectangular block with the left and right ends free from stresses and the fluid is free to flow, whereas an external force is applied on the top and bottom boundaries. The external force is taken as 10.5 MPa. The Young’s moduli along horizontal and vertical directions are 20 and 10 GPa, respectively. In-plane Poisson’s ratio is assumed to be 0.30, whereas the out-of-plane one is taken as 0.20. Similarly, permeability values are taken to be 10^{−19 }and 10^{−17}m^{2} along horizontal and vertical directions. The bulk moduli of the solid grains and the pore fluid are assumed to be 36 and 3 GPa. A time step of 500 s is used.

**Fig. 7** a Pore pressure variation along the centre line at different times, b relative error (%) in pore pressure when compared to the analytical solution and c displacement in the x-direction along the centre line at time t=75,000s

We compare the pore pressure solution from the numerical model with the analytical solution at different time periods in Fig. 7. The numerical pore pressure decay from the centre of the specimen to the free edges is found to be consistent with the analytical solution with relative errors (< 5%). The displacement in the x-direction along the centre line of the specimen is plotted and compared with the analytical solution.

#### 5.1.2 KGD

In this validation case we consider a KGD problem (Fig. 8) which is a standard test case for hydraulic fracture problems. When the rock is assumed to be isotropic, there exists a theoretical solution given by Geertsma and Klerk (1969). The Young’s modulus and Poisson’s ratio are taken as 20 MPa and 0.2. The ultimate tensile strength is assumed to be 2 MPa, while the toughness is 120{N/m}.

**Fig. 8** Geometry and boundary conditions of the standard KGD problem

The permeability and viscosity are given as 10^{−19}m^{2} and 0.1Pas. The initial fracture at the boundary is injected at the rate of 25 mm^{2}/s for a time period of 100 s with a time step of 0.1 s. The KGD problem considered here lies in the viscosity–storage propagation regime.

**Fig. 9** a Comparison of fracture profile variation at different times (t=12.5,25,50,100 s) and b fracture mouth opening pressure variation with time

The fracture propagates on a non-predefined path. In Fig. 9, we compare the fracture profiles at various time steps from the current numerical model and the ELP model (Remij et al. 2015). As observed the current model accurately reduces to the ELP solution for isotropic values of the parameters. The fracture mouth opening pressure variation with time is also plotted and compared. A mesh of 50 × 50 mm^{2 }is used for the purpose.

### 5.2 Vertical Hydraulic Fracture Problem

The test case (Fig. 10) that is used here is very similar to the previous KGD problem, but an initial crack is placed at the bottom of the model, to model the vertical fracture growth representative of the hydraulic fractures. The model is simulated for 10.5 s with a time step of 0.1 s. Also the model is subjected to external stresses of 40 MPa, due to the overburden pressures or in situ stresses existing at depths of around 1.5 km.

**Fig. 10** Geometry and boundary conditions of the hydraulic fracture problem

The fluid is free to flow from the top and the right boundaries, whereas no-flow conditions exist at the left and bottom boundaries. The isotropic values of the parameters used in the model are specified in Table 1. All angles are with respect to the horizontal axis. For cases of anisotropy, these isotropic values are perturbed depending on the degree of anisotropy. The degree of anisotropy is defined as:

**Table 1** Isotropic value of parameters

By using the parameters in Table 1 to obtain the non-dimensional parameter (Mk) described in Bunger et al. (2005), we observe that the hydraulic fracture problem described here lies in the viscosity-dominated regime and closer to the storage edge. In the following Sects. 5.3 and 5.6.1 we try to understand the influence of anisotropy in each individual parameter by keeping all other parameters isotropic. We also look at the possible combination of anisotropy in these parameters in Sects. 5.4 and 5.6.2.

### 5.3 Parametric Anisotropy

In this subsection, we vary one parameter at a time to find out the fracture propagation variation with anisotropy in each individual parameter. In all the considered test cases we assume that the grains are oriented along the horizontal direction (0^{∘}) and the initial fracture is oriented in the vertical direction (90^{∘}).

#### 5.3.1 Anisotropy Due to Young’s Modulus

We consider three possible scenarios for varying Young’s modulus: (1) E-Parallel: anisotropy caused by increasing the Young’s modulus (E_{∥}) parallel to the grain direction alone. (2) E-Perpendicular: anisotropy caused by decreasing the Young’s modulus (E_{⊥}) perpendicular to the grain direction alone. (3) E-Combined: anisotropy caused by varying both parallel and perpendicular values from the isotropic values given by following equation:

Both in-plane and out-of-plane Poisson’s ratios are assumed to have a constant value of 0.2.

From the fracture length plot in Fig. 11, we can observe that E_{∥} has a much greater effect on fracture propagation than E_{⊥}. This is due to the fact that the propagation of the initial vertical fracture is dependent on the stresses which are perpendicular to it. Hence, an increase in E∥ results in higher stresses perpendicular to the initial fracture, which promotes fracture growth significantly, whereas a decrease in E_{⊥} only has a smaller Poisson’s effect on the stress. Since E_{∥}>E_{⊥ }always, the fracture prefers to orient itself perpendicular to the grain orientation which is observed in all the three scenarios in the fracture orientation plot.

Also we plot the pressure at the mouth of the fracture, which is a much more easily measurable quantity in the field. Since all the three scenarios favour fracture propagation in the same initial fracture direction, there is little variation (<5%) in the pressure required to open the fracture.

**Fig. 11** Fracture variation with Young’s modulus anisotropy after 20 s for an equivalent injection volume of 0.012m^{2}

**Fig. 12** Fracture variation with ultimate tensile strength anisotropy after 20 s for an equivalent injection volume of 0.012m^{2}

#### 5.3.2 Anisotropy Due to Ultimate Tensile Strength

Similar to the Young’s modulus variation, we consider the same three scenarios for understanding ultimate tensile strength-induced anisotropy. Fracture propagation is resisted by the ultimate tensile strength perpendicular to the fracture orientation. Hence, the fracture tends to propagate along the direction perpendicular to the minimum ultimate tensile strength. Since τ_{ult}_{⊥ }is always lower than τ_{ult}_{∥}, the fracture tends to propagate parallel to the grain orientation. But since the initial fracture is oriented in an unfavourable direction (perpendicular to the grain direction), the fracture continues to propagate in its initial direction until a threshold level where the effect of anisotropy becomes significant to rotate the fracture as observed from the fracture orientation plot in Fig. 12.