## Abstract

Extracting gas from shale rocks is one of the current engineering challenges but offers the prospect of cheap gas. Part of the development of an effective engineering solution for shale gas extraction in the future will be the availability of reliable and efficient methods of modelling the development of a fracture system, and the use of these models to guide operators in locating, drilling and pressurising wells. Numerous research papers have been dedicated to this problem, but the information is still incomplete, since a number of simplifications have been adopted such as the assumption of shale as an isotropic material. Recent works on shale characterisation have proved this assumption to be wrong.

#### Authors:

Gabriel Hattori^{1}, Jon Trevelyan^{1}, Charles E. Augarde^{1}, William M. Coombs^{1}, Andrew C. Aplin^{2}

^{1}School of Engineering and Computing Sciences, Durham University, South Road, Durham DH1 3LE, UK. ^{2}Science Labs, Department of Earth Sciences, Durham University, Durham DH1 3LE, UK

Received: 2 September 2015, Accepted: 6 January 2016, Published online: 29 January 2016

©The Author(s) 2016.

The anisotropy of shale depends significantly on the scale at which the problem is tackled (nano, micro or macroscale), suggesting that a multiscale model would be appropriate. Moreover, propagation of hydraulic fractures in such a complex medium can be difficult to model with current numerical discretisation methods. The crack propagation may not be unique, and crack branching can occur during the fracture extension.

A number of natural fractures could exist in a shale deposit, so we are dealing with several cracks propagating at once over a considerable range of length scales. For all these reasons, the modelling of the fracking problem deserves considerable attention. The objective of this work is to present an overview of the hydraulic fracture of shale, introducing the most recent investigations concerning the anisotropy of shale rocks, then presenting some of the possible numerical methods that could be used to model the real fracking problem.

## 1 Introduction

Conventional shale reservoirs are formed when gas and/or oil have migrated from the shale source rock to more permeable sandstone and limestone formations. However, not all the gas/oil migrates from the source rock, some remaining trapped in the petroleum source rock. Such a reservoir has been named “unconventional” since it has to be fractured in order to extract the gas from inside. Hydraulic fracture, or “fracking”, has emerged as a alternative method of extracting gas and oil. Experience in the United States shows it has the potential to be economically attractive. Many concerns exist about this type of extracting operation, especially how far the fracture network will extend in shale reservoirs.

King [136] published a review paper about the last 30 years of fracking, and points out four “lessons”:

- No two shale formations are alike. Shale formations vary spatially and vertically within a trend, even along the wellbore;
- Shale “fabric” differences, combined with in-situ stresses and geologic changes are often sufficient to require stimulation changes within a single well to obtain best recovery;
- Understanding and predicting shale well performance requires identification of a critical data set that must be collected to enable optimization of the completion and stimulation design;
- There are no optimum, one-size-fits-all completion or stimulation designs for shale wells.

These points encapsulate well the uncertainties involved. Many models have been proposed over the years but they are either too simplified or they tend to focus on one key aspect of fracking (e.g. crack propagation schemes, influence of natural fractures, material heterogeneities, permeabilities). The scarcity of in-situ data makes the study of fracking even more complicated.

The most usual concerns in fracking are addressed by Soeder et al. [252], where integrated assessment models are used to quantify the engineering risk to the environment from shale gas well development. Davies et al. [55] have investigated the integrity of the gas and oil wells, analysing the number of known failures of well integrity. The modelling of reservoirs is also a difficult task due to the lack of experimental data and oversimplification of the complex fracking problem [177].

Glorioso and Rattia [97] provide an approach more focused on the petrophysical evaluation of shale gas reservoirs. Some techniques are analysed, such as log responses in the presence of kerogen, log interpretation techniques and estimation methods for different volumes of gas in-situ, among others. It is shown that volumetric analysis is imprecise for in-place estimation of shale gas; however, it is one of the few techniques available in the early stages of evaluation and development. The measurement of an accurate density of specimens is an important parameter in reducing the uncertainty inherent in petrophysical interpretations.

This paper provides an overview of the current state of fracking research. A state-of-the-art review of fracking is performed, and several points are analysed such as the models employed so far, as well as the underlying numerical methods. Special attention is given to problems involving brittle materials and the dynamic crack propagation that must be taken into account in the fracking model. The hydraulic fracture modelling problem has been tackled in several different ways, and the shale rock has mostly been assumed to be isotropic. This simplification can have serious consequences during the modelling of the fracking process, since shale rocks can present high degrees of anisotropy.

This paper is organised as follows: a description of the shale rock including the most common simplifications is presented in Sect. 2, followed by the description of the fracking operation in Sect. 3. Section 4 presents a review of the analytical formulations for crack propagation and crack branching. Different types of models such as cohesive methods and multiscale approaches are tackled in Sects. 5 and 6. Numerical aspects are discussed in Sect. 7, including the boundary element method, the extended finite element method, the meshless method, the phase-field method, the configurational force method and the discrete element method. A recently proposed discretisation method is discussed in Sect. 8. The paper ends with conclusions and a discussion of possible future research directions in Sect. 9.

## 2 Description of the Shale Rock

Shale, or mudstone, is the most common sedimentary rock. It can be viewed as a heterogeneous, multi-mineralic natural composite consisting of sedimented clay mineral aggregates, organic matter and variable quantities of minerals such as quartz, calcite and feldspar. By definition, the majority of particles are less than 63 microns in diameter, i.e. they comprise silt- and clay-grade material. In the context of shale gas and oil, organic matter (kerogen) is of particular importance as it is responsible for the generation and, in part, the subsequent storage of oil and gas.

Mud is derived from continental weathering and is deposited as a chemically unstable mineral mixture with 70–80 % porosity at the sediment-water interface. During burial to say 200 °C and 100 MPa vertical stress, it is transformed through a series of physical and chemical processes into shale. Porosity is lost as a result of both mechanical and chemical compaction to values of round 5 % [31, 32, 287]. At temperatures above 70 °C, clay mineral transformations, dominated by the conversion of smectite to illite (e.g. [121, 254]), lead to a fundamental reorganisation of the clay fabric, converting it from a relatively isotropic fabric to one in which the clay minerals are preferentially aligned normal to the principal (generally vertical) stress [56, 57, 120].

Although quantitative mechanical data are scarce for mineralogically well-characterised samples, it is likely that the clay mineral transformations strengthen shales [206, 264]. In muds which contain appreciable quantities of biogenic silica (opal-A) and calcite, the conversion of opal-A to quartz [134, 281], and dissolution-reprecipitation reactions involving calcite [259], will also strengthen the shale. Indeed, it is generally considered that fine-grained sediments which are rich in quartz and calcite are more attractive unconventional oil and gas targets compared to clay-rich media, as a result of their differing mechanical properties (e.g. [204]).

Shales with more than ca. 2 % organic matter act as sources and reservoirs for hydrocarbons. Between 100 and 200 °C kerogen is converted to hydrogen-rich liquid and gaseous petroleum, leaving behind a carbon-rich residue (e.g. [126, 144, 207]). The kerogen structure changes from more aliphatic to more aromatic, and its density increases [194]. Changes in the mechanical properties of kerogen with increasing maturity are not well documented. However, they may be quite variable, depending on the nature of the organic matter. For example, Eliyahu et al. [68] performed PeakForce QNM® tests with an atomic force microscope to make nanoscale measurements of the Young’s modulus of organic matter in a single shale thin section. Results ranged from 0-25 GPa with a modal value of 15 GPa.

Shales are heterogeneous on multiple scales ranging from sub-millimetre to tens of metres (e.g. [10, 204]). Hydrodynamic processes associated with deposition often result in a characteristic, ca. millimetre-scale lamination [35, 157, 241], which can be disturbed close to the sediment-water interface by bioturbation [63]. On a larger, metre-scale, parasequences form within mud-rich sediments, driven by orbitally-forced changes in climate, sea-level and sediment supply [35, 156, 157, 204].

Parasequence boundaries are typically defined by rapid changes in the mineralogy and grain size of mudstones, with more subtle variations within the parasequence. Stacked parasequences add further complexity to the shale succession and result in a potentially complex mechanical stratigraphy which depends on the initial mineralogy of the chosen unit and the way that burial diagenesis has altered physical properties on a local scale.

During the shale formation process bedding planes are formed, which may present sharp or gradational boundaries. This is the most regular type of deposition that occurs in shales. Deposition may not be uniform during the whole process, presenting discontinuities at some points or other type of deposition patterns. This makes the mechanical characterisation of shale a complex issue. Moreover, not all shale rocks are the same, so a prediction made for an specific shale rock probably is not valid elsewhere.

The works of Ulm and co-workers about nanoindentation in shale rocks [34, 198, 199, 200, 268, 269, 270] have been important developments in our ability to characterise the mechanical properties of shale rocks. From [268], it is seen that shales behave mechanically as a nanogranular material, whose behaviour is governed by contact forces from particle-to-particle contact points, rather than by the material elasticity in the crystalline structure of the clay minerals. This assumption is valid for scales around 100 nm.

The indentation technique consists of bringing an indenter of known geometry and mechanical properties (typically diamond) into contact with the material for which the mechanical properties are to be known. Through measurement of the penetration distance h as a function of an increasing indentation load P, the indentation hardness H and indentation modulus M are given by

where Ac is the projected area of contact and S = (dP/dh)_{hmax} is the unloading indentation stiffness. For the case of a transversely isotropic material, where x^{3} is the axis of symmetry, the indentation modulus is given by [268, 269]

and for the x_{1} ; x_{2} axis by

_{ij}come from the constitutive matrix and are given in the Voigt notation [276].

From [269], it was seen that the level of shale anisotropy increases from the nanoscale to the macroscale. Macroscopic anisotropy in shale materials results from texture rather than from the mineral anisotropy. The multiscale shale structure can be divided into 3 levels:

- Shale building block (level I – nanoscale): composed of a solid phase and a saturated pore space, which form the porous clay composite. A homogeneous building block, which consists in the smallest representative unit of the shale material, is assumed at this scale. The material properties are composed of two constants for the isotropic clay solid phase, the porosity and the pore aspect ratio of the building block.

- Porous laminate (level II – microscale): the anisotropy increases due to the particular spatial distribution of shale building blocks (considering different types of shale rocks). The morphology is uniform allowing the definition a Representative Volume Element (RVE).
- Porous matrix-inclusion composite (level III – macroscale): shale is composed of a textured porous matrix and (mainly) quartz inclusions of approximately spherical shape that are randomly distributed throughout the anisotropic porous matrix. The material properties are separated into six indentation moduli plus the porosity.

One can observe that the heterogeneities are manifested from the nanoscale to the macroscopic scale, and combine to cause a pronounced anisotropy and large variety in shale macroscopic behaviour.

Nanoindentation results provide strong evidence that the nano-mechanical elementary building block of shales is transversely isotropic in stiffness, and isotropic and frictionless in strength [34]. The contact forces between the sphere-like particles activate the intrinsically anisotropic elastic properties within the clay particles and the cohesive bonds between the clay particles.

The determination of the mechanical microstructure and invariant material properties are of great importance for the development of predictive microporomechanical models of the stiffness and strength properties of shale.

## 3 The Hydraulic Fracturing Process and Its Modelling

The hydraulic fracture or fracking operation involves at least three processes [3]:

- The mechanical deformation induced by the fluid pressure on the fracture surfaces;

- The flow of fluid within the fracture;

- The fracture propagation

The shale measures in question are usually found at a distance of 1 to 3 km from the surface. A major concern relating to fracking is that the fracture network may extend vertically, allowing hydrocarbons and/or proppant fluid to penetrate into other rock formations, eventually reaching water reservoirs and aquifers that are found typically approximately 300 m below the surface.

Fracking can occur naturally, such as in magma-driven dykes for example. In the 1940s, when fracking started commercially in US, the hydraulic fracture was applied through a vertical drilling. In that case, the pressurised liquid was applied perpendicular to the bedding planes. It was known that the shale was a layered material due to its formation process, but technology of that time was very limited.

**Fig. 1** Fracking example.

In the last 15 years, recent engineering advances have allowed engineers to change the direction of the drilling, making it possible to drill a horizontal well and consequently, to pressurise the shale rocks in the same horizontal plane of the bedding plane, making the fracking process much more effective. Figure 1 illustrates the structure of the well’s drilling, and the natural fracture network that can be found. In detail it is a sketch of the pressurised liquid entering a crack, resulting in the application of a pressure P over the crack surfaces and the crack opening w.

The horizontal drilling was not new to the industry, but it was fundamental for the success of shale gas developments. From 1981 to 1996, only 300 vertical wells were drilled in the Barnett shale of the Fort Worth basin, north central Texas. In 2002, horizontal drilling has been implemented, and by 2005 over 2000 horizontal wells had been drilled [40]. The Barnett shale formation found in Texas produces over 6 % of all gas in continental United States [273]. The application of this new drilling technique has turned the United States from a nation of waning gas production to a growing one [221].

To optimise the fracking process of shale, it is important to detect accurately the location of natural fractures. The anisotropic behaviour of the shale generates preferential paths through the shale fabric [136, 279]. Moreover, the alignment of the natural fractures can also induce anisotropic patterns of the fluid flow [86, 87].

### 3.1 Modelling of the Shale Fracture

Much of the work done so far in attempting to model shale fracture is very simple, taking into account only the influence of the crack and not the fluid. Only recently have a few researchers [3, 4, 62, 160, 161, 188, 205, 296] successfully developed more sophisticated methods including the fluid-crack interaction.

The usual assumption in hydraulic fracture is that the fracture is embedded within an infinite homogeneous porous medium, where flow occurs only perpendicular to the fracture plane, which was first defined by Carter [46]. Moreover, the injection pressure does not propagate beyond the current extent of the fracture. Carter’s model can lead to an overestimate of the fracture propagation rate by a factor of 2 as compared to a 3D model [161]. The reason is that the pressure increases beyond the length of the hydraulic fracture, causing an increasing of the leak-off and a corresponding reduction in fracture growth. The leak-off rate Q_{1} is given as [161]

where µ is the fluid viscosity, k_{α} is the permeability in the α-direction α= r or α= z, P is the hydraulic pressure, a(t) is the hydraulic fracture radius, dependent of time t, r and z are the distances parallel and normal to the fracture plane, respectively. The hydraulic pressure is defined by the boundary value problem

where S is a storage coefficient of the porous medium. The solution of Eq. (6) can be obtained using a standard finite volume method, as used in [161].

Assuming that the faces of the fracture are loaded by a uniform pressure P_{d}, the displacement of the fracture face normal to the fracture plane δ is given by

where υ is the Poisson’s ratio and E is the Young’s modulus. The fracture volume V is found from

The energy release rate G of the rock is obtained from the following expression

and is related to the mode I stress intensity factor K_{I} through the expression

From Eqs. (8) and (9), it is possible to write Pd and a in terms of V as

In the early stages following initial pressurisation, the volume of injected fluid is sufficiently small such that the the porous formation do not absorb the incoming fluid. As injection process continues, the fluid is accommodated locally in the pore space and consequently predicts leak-off. Once the system reaches steady-state, it again becomes independent of porosity system. This analytical formulation have issues when predicting the behaviour during the transient state [161].

Even though these models can represent complex processes occurring during fracking, they are still far from being accurate, mainly because shale is considered to be isotropic, which has been seen not to be true [268], and since the material presents nanoporosity, it is difficult to accurately model the mechanical properties of shale.

Some other open questions are [3]:

- how to appropriately adjust current (linear elastic) simulators to enable modelling of the propagation of hydraulic fractures in weakly consolidated and unconsolidated “soft” sandstones;
- laboratory and field observations demonstrate that mode III fracture growth does occur, and this needs to be further researched.

Some works have analysed the crack propagation path in shales, including refracturing sealed wells. For example, Gale and co-workers [87] found that propagation of the hydraulic fracture over a natural fracture will cause delamination of the cement wall and the shale. The fluid enters the fracture and causes further opening of the fracture in a direction normal to the propagating hydraulic fracture while the pressure inside the fracture increases. After the fracture propagation at the natural fracture reaches a sealed fracture tip, the hydraulic fracture resumes growth parallel to the direction of maximum shear stress.

In an analytical work, Vallejo [271] has investigated the hydraulic fracture on earth dam soils, where shear stresses were seen to promote crack propagation on traction free cracks. Other analytical study about re-fracturing was carried out by [295], where the dynamic fracture propagation is characterised in low-permeability reservoirs. The results are comparable to an experimental test with the same material parameters.

In summary, research works in hydraulic fracture formulation have considered a large number of variables and processes which occur during the actual operation: leak-off, shale permeabilities, crack opening and fluid interaction over a crack surface. However, the current analytical theories for hydraulic fracture do not include crack propagation conditions, especially dynamic crack propagation, neither crack branching, since material instabilities at the crack tip during crack propagation may cause the propagation path not to be unique. These concerns are summarised in the next section.

## 4 Crack Propagation and Crack Branching

Consider a homogeneous isotropic body under a known applied loading. The resulting elastic stress distribution over the body due to the applied force is generally smooth. However, introducing a discontinuity such as a crack imposes a singular behaviour to the stress distribution. It can be shown that the stress increases as it is measured closer to the crack tip, varying with 1/r√, where r is the distance from the crack tip. Irwin [125] proposed that the asymptotic stress field at the crack tip is governed by parameters depending on the geometry of the crack and the applied load. These parameters are known as Stress Intensity Factors (SIFs) and have been widely used as criteria for crack stability and propagation. The three SIFs, K_{I},K_{II},K_{III}, each correspond to one of three modes of crack behaviour: mode I (opening), mode II (sliding) and mode III (tearing). In this paper we will confine ourselves mostly to mode I.

It can be postulated that crack growth will begin if the value of the SIFs increase to a certain value. If the SIF is higher than a critical fracture toughness parameter Kc, which depends on the material properties, then the crack will propagate through the body. The situation becomes more complicated when the load is applied rapidly so that dynamic effects become important. This does not imply that the value of the dynamic fracture toughness will be independent of the rate of loading or that dynamic effects do not influence the fracture resistance in other ways [82].

In some cases, the toughness appears to increase with the rate of loading whereas in other cases the opposite dependence is found. The explanation for the shift must be sought in the mechanisms of inelastic deformation and material separation in the highly stressed region of the edge of the crack in the loaded body [82]. The dynamic crack propagation formula can be defined as

where K_{d}^{I} is the mode I dynamic SIF, K_{I} is the static SIF, v is the crack velocity and κ(v) is a scaling factor. When κ(v)=1, the crack velocity is zero, whereas κ(v)=0 indicates that the crack velocity is equal to the Rayleigh wave speed.

The theoretical limiting speed of a tensile crack must be the Rayleigh wave speed. This was anticipated by Stroh [257] on the basis of a very intuitive argument [82].

Gao et al. [89] studied crack propagation in an anisotropic material, and presented expressions for the dynamic stresses and displacements around the crack tip. These predict that larger crack propagation velocities induce higher stress and displacement fields at the crack tip. The limiting speed in crack propagation is analysed in [88], where a local wave speed resulting from the elastic response near the crack tip also changes with the crack propagation velocity. A molecular dynamic model is used in this work, so crack propagation is modelled as bond breakage between the particles. The crack velocity is expressed using the Stroh formalism.

There are three types of criteria for brittle crack propagation:

1-Maximum tangential stress: This criteria was defined by Erdogan and Sih [69] and is based in two hypothesis:

- The crack extension starts at its tip in radial direction;
- The crack extension starts in the plane perpendicular to the direction of greatest tension.

The crack propagates when the SIF is higher than a critical SIF Kc, which depends on the materials properties. From [69], the crack propagation angle θp can be obtained from the following relation

where K_{I} and K_{II} are the mode I and mode II SIF, respectively, and θp is taken with respect to the horizontal axis. This crack propagation criteria was extended to anisotropic materials in [239].

2-Strain energy release rate: In this criteria, the crack propagates when energy release rate G reaches some critical value Gc, taking the direction where G is maximum [123]. The energy release rate is defined as

where W represents the strain energy and a is the half-crack length. Equation (15) can be expressed in terms of mixed mode SIFs for an isotropic material as

where μ is the shear modulus.

3-Minimum strain energy: crack propagation occurs at the minimum value of the strain density S defined as [169, 247]

where a_{ij} come from the material properties. The direction of propagation goes toward the region where S assumes a minimum value Smin. The crack extension r_{0} is proportional to the minimum strain energy, such that the ratio S_{minr0 }is constant along the crack front [169].