Numerical Simulation of Fracking in Shale Rocks: Current State and Future Approaches

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. 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.


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.

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, sealevel 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 A c is the projected area of contact and S ¼ ðdP=dhÞ h max 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 where C 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: 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.

The Hydraulic Fracturing Process and Its Modelling
The hydraulic fracture or fracking operation involves at least three processes [3]: 1. The mechanical deformation induced by the fluid pressure on the fracture surfaces; 2. The flow of fluid within the fracture; 3. 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.
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].

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 leakoff rate Q 1 is given as [161] where l is the fluid viscosity, k a is the permeability in the a-direction (a ¼ r or a ¼ 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 d is given by where t 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 P d 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 leakoff. 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 natural fracture network 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: leakoff, 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.

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= ffiffi r p , 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 K c , 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 jðvÞ is a scaling factor. When jðvÞ ¼ 1, the crack velocity is zero, whereas jð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: a. The crack extension starts at its tip in radial direction; b. 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 K c , which depends on the materials properties. From [69], the crack propagation angle h p can be obtained from the following relation where K I and K II are the mode I and mode II SIF, respectively, and h 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 G c , 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 halfcrack length. Equation (15) can be expressed in terms of mixed mode SIFs for an isotropic material as where l is the shear modulus. 3. Minimum strain energy: crack propagation occurs at the minimum value of the strain density S defined as [169,247] S ¼ a 11 K 2 I þ 2a 12 where a ij come from the material properties. The direction of propagation goes toward the region where S assumes a minimum value S min . The crack extension r 0 is proportional to the minimum strain energy, such that the ratio S min r 0 is constant along the crack front [169].
One can observe that all these criteria are related to the SIFs. These criteria are well consolidated in the fracture mechanics literature over the years. However they fail in one aspect, since they do not consider the possibility of crack branching, i.e., at some point of the crack propagation process, the crack may bifurcate in two or more new cracks. This issue is especially important when modelling highly heterogeneous materials such as the shale rock.
Yoffe [289] attempted to explain the branching of cracks from an analysis of the problem of a crack of constant length that translates with a constant velocity in an unbounded medium. From this solution she found that the maximum stress acted normal to lines that make an angle of 60 with the direction of crack propagation when the crack velocity exceeded 60 % of the shear wave speed. This fact might cause the crack to branch whenever the crack velocity exceeds that value. However, Yoffe did not consider that the maximum stress would be perpendicular to the crack path, so this assumption is not valid for brittle materials. Moreover, the 60 angles are quite large in comparison with the branching angles observed from experiments [223].
Ravi-Chandlar and Knauss [222,223] have addressed the crack propagation and crack branching problems through several experiments. From [222], the crack branching has the following properties 1. crack branching is the result of many interacting microcracks or microbranches; 2. only a few of the microbranches grow larger while the rest are arrested; 3. the branches evolve from the microcracks which are initially parallel to the main crack, but deviate smoothly from the original crack orientation; 4. the microbranches do not span the thickness of the plate, some occurring on the faces of the plate while others are entirely embedded in the interior of the plate.
Sih [247] made the hypothesis that the instability that occurs in crack bifurcation is associated with the fact that a high speed crack tends to change its direction of propagation when it encounters an obstacle in the material. The excess energy in the vicinity of such a change in direction is sufficient to initiate a new crack. This event occurs so quickly that the crack appears to have been split in two, or bifurcated.
From [223], one can see that the velocity with which the crack propagates is determined by the SIF at initiation. Cracks propagating at low speeds may undergo a change in the crack velocity if stress waves are present. Cracks propagating at high speeds do not change crack velocity, but may exhibit crack branching.
Crack branching formulations can be found in [78,131,247,289], to cite just a few works. In all cases, only the isotropic material case is considered. For anisotropic crack branching, numerical methods have to be employed.

Cohesive Methods
The fracture process is usually considered only at the crack tip. In such cases, the fracture process zone is considered to be small compared to the size of the crack [17,66,67].
Numerical Simulation of Fracking in Shale Rocks: Current State and Future Approaches 287 In Linear Elastic Fracture Mechanics (LEFM), the stress becomes infinite at the crack tip. Since no material can withstand such high stress, there will be a plastification/ fracture zone around the crack tip.
The fracture process zone can be described by two simplified approaches [178]: 1. The fracture process zone is lumped into the crack line and is characterised as a stress-displacement law with softening; 2. Inelastic deformations in the process zone are smeared over a band of a given width, imagined to exist in front of the main crack.
Most of the work done in cohesive cracks makes use of the former approach, otherwise known as the Dugdale-Barrenblatt model, fictitious crack model or stress bridging model [178].
The non-linear behaviour around the crack tip can be considered to be confined to the fracture process zone on the crack surface. Figure 2 illustrates a crack with its corresponding fracture process zone. One can see two tips in this model, the physical tip, where the tractions vanish, and the fictitious tip, where the displacement is zero.
Since there is no singularity at the crack tip, the SIF should vanish. This condition is also called the zero stress intensity factor, and is represented by the superposition of two states where K phys I corresponds to the SIF at the physical crack tip, and K fict I is the SIF at the fracture process zone. Here we consider only the mode I fracture without loss of generality.
The crack propagates when the maximum principal stress reaches the material tensile strength r t , so fracture is initiated at the fracture process zone. The stress on the crack faces depends directly on the relative displacement Du of the crack faces [213]. There are different types of stress-displacement functions which model the behaviour in the fracture process zone. Figure 3 presents two of the most common assumptions Dugdale [66] and Barenblatt [17] models are the basis of many cohesive models. The Dugdale cohesive crack model is very simplistic and is best used for ductile materials. A uniform traction equal to the yield stress is used to describe the softening in the fracture process zone.
Most of the cohesive models are developed for isotropic materials (see [67] for example). However, there are some models for heterogeneous materials [11,216,244] and composite [148,181,261,272] materials. Nevertheless, the material models are quite simple, usually considering different types of isotropic materials instead of a full anisotropic model. To the authors' best knowledge, there is no anisotropic cohesive crack model to this date.
Cohesive models have been also applied in multiscale problems, where cracks are significantly smaller than the RVE. In [210], a microelastic cohesive model is developed for quasi-brittle materials. The stability of crack growth is analysed, and it is concluded that macroscopic strength is not necessarily correlated to crack propagation, and may be caused by unstable growth of cohesive zones ahead of nonpropagating cracks. The initial cohesive zone has a significant influence on the macrostrength of quasi-brittle materials.
A number of different approaches for cohesive models have been proposed over the years. Enriched formulations for delamination problems were analysed by Samimi [236][237][238]. A stochastic approach for delamination in composite materials was proposed in [181], where the imperfections of the material were considered in the cohesive model. The cohesive crack has been extensively studied as can be seen in [59,75,102,152,209,297] to cite just a few works.  Crack propagation in cohesive models was recently discussed in [152,298] for example.
A rudimentary model for hydraulic fracture for isotropic materials using the finite element software ABAQUS was considered in [288]. Papanastasiou [203] has evaluated the fracture toughness in hydraulic fracturing, modelling the rock-fluid coupling through a finite element model with cohesive behaviour. The Mohr-Coulomb criterion is used to take plasticity into account in the rock deformation. The plastic behaviour that develops around the crack tip provides an effective shielding, resulting in an increase in the effective fracture toughness.

Multiscale
The main advantage of multiscale models is to make different hypotheses at different levels within the same problem. For example, material can present distinct degrees of anisotropy depending on the scale of observation (nano, micro or macro). The coupling of different scales can be cumbersome. Some sort of regularisation is commonly used to enforce the coupling between scales. A typical assumption is the use of a Representative Volume Element (RVE), a representative part of the model at the reduced scale so it contains all the distinct properties of the considered scale, and is also defined as the local model. The global model takes the RVE as a homogenised representation of the material's properties at the large scale. An example of an RVE is illustrated in Fig. 4.
Another important part of multiscale modelling is the coupling of stresses and strains from the local and global models. Numerical homogenisation is a popular technique and is an alternative to the traditional analytical homogenisation. It is especially used for monophasic heterogeneous materials, where the balance and constitutive equations are considered at the RVE level. The first work in numerical homogenisation is due to Ghosh et al. [96].
Zeng et al. [292] proposed a multsicale cohesive model for geomaterials. At the macroscopic scale, a sample of polycrystalline material is considered as a continuum made of many material points. The estimation of the material properties at the microscale is performed by statistical homogenisation, since the RVE represents a number of different constituents or phases, as mineral grains and voids, and is therefore composed of randomly distributed constituents.
The Eshelby elastic solution for the spherical inclusion problem [71,72] is used to obtain the local stress and strain fields. Therefore, the strains or stresses in a single crystal are approximated by considering a spherical single crystal embedded in an infinite elastically deformed matrix. The KBW model, named after Kröner [140], Budiansky and Wu [43], extends Eshelby's formulation by taking into account the grain interaction and plastification. By the KBW definition, each crystal is embedded in a Homogeneous Equivalent Medium (HEM) as shown in Fig. 5.
The local stress r and strain e are related to the global stress R and strain E as follows where L is the interaction tensor and is given by where M is the homogenised elastoplastic tangent operator of HEM, S is the Eshelby's tensor and I is a third order identity matrix. Zeng and Li [293] developed a multiscale cohesive zone method, where the local fields are determined through measures of the bond at the atom particle level. The stress relation coupling the local and global fields is given by where X b is the volume of the unit cell, n b is the total number of bonds in a unit cell, /ðr i Þ is the atomistic potential, r i ; i ¼ 1; Á Á Á ; n b is the current bond length for the ith bond in a unit cell and is given by with the deformation gradient F e in element e and the underformed bond vector R. The symbol denotes the outer, or dyadic, product. The strain energy in a given element X e can be written as and therefore the total energy is defined as n a E a ðu a Þ ð 23Þ where n a is a chosen weight and P n a ¼ 1. The energy from each representative atom E a is obtained by the interaction with the neighbouring atoms whose positions are generated using the local deformation.
This formulation is referred to as the local QC method, a simplification of the continuum system when interface and surface energies may be neglected. The general non-local QC potential energy may lead to some non-physical effects in the transition region. The derivatives of the energy functional to obtain forces on atoms and finite element nodes may lead to so-called ghost forces in the transition region between the macro and microscale, and it has several issues that remain to be resolved, such as the computation of approximations in the macroscale far from microscale defects [245] and the correct balance of energy which needs to be ensured between macro and microscales [174].
Since the connections between atoms are modelled through bonds, this multiscale cohesive formulation is able to capture the crack branching behaviour during crack propagation.
In [299], the RVE properties of a hydrogeologic reservoir are averaged through statistical parameters. The main reason is that the heterogeneity of the reservoir can be more easily modelled through the mean and standard deviation of the rock properties. The site scale represents the entire solution domain used for modelling global flow and transport. The layer scale represents geologic layering in the vertical direction. Within a layer, relatively uniform properties are present in both vertical and lateral direction, in comparison with the larger variations between different layers that may vary significantly in thickness. The local scale represents the variation of properties within a hydrogeologic layer.
In [164], a multiscale model for the shale porous network is proposed. Permeability is assumed as an intrinsic porous medium property independent of fluid properties (such as viscosity) or thermodynamic conditions. The porous medium was modelled as networks of pores connected by throats. This simplification neglects the physics of the real porous network. Permeability further depends on the relative size of the void spaces as well as the fraction of pores belonging to each length scale. Unlike absolute permeability in conventional reservoirs, gas permeability depends on absolute pressure values in individual pores (and not only the gradient). Specifically, smaller pressures result in (somewhat counter-intuitively) an increase in permeability.
A number of multiscale models for brittle materials can be mentioned: [2,83,130,209,274] just to cite some of the most recent works.

Discrete Numerical Methods
In this section we will present a brief description of different element-based numerical methods that can be used in the modelling of fracking problems. The boundary element method (BEM) has been used in brittle anisotropic problems including crack propagation. The extended finite element method (X-FEM) has been developed recently and is also a good choice for fracture mechanics problems, and can be easily applied in cohesive models. Meshless methods are becoming popular in fracture mechanics problems. The discrete element method (DEM) is particularly used in problems with rock materials. The phase-field method and the configurational force method are also reviewed in this section.

Boundary Element Method (BEM)
The boundary element method has first appeared in the work of Cruse and Rizzo [52] for elasticity problems, but it was effectively named as BEM in the work of Brebbia and Domínguez [41] and represented a series of advances in comparison to the existent domain discretisation methods as the finite element method (FEM) and the finite differences method (FDM) [109]: • Accurate mathematical representations of the underlying physics are employed, resulting in the ability of the BEM to provide highly accurate solutions; • The problem is defined only at the boundaries, which gives a reduction of dimensionality in the mesh (linear for 2D problems and surface for 3D problems), therefore resulting in a reduced set of linear equations to be solved; • In spite of the boundary-only meshing, results at any internal point in the domain can be calculated once the boundary problem has been solved; • There is a great advantage in certain classes of problem that can be characterised by either (1) infinite (or semiinfinite) domains, or (2) discontinuous solution spaces. These advantages have resulted in the BEM gaining popularity for acoustic scattering, fracture mechanics, re-entry corners and other stress intensity problems, where domain discretisation methods have poorer convergence.
However, there are some drawbacks which may deterred FEM users from migrating to the BEM: • The system of equations is both non-symmetric and fully populated, which may lead to longer computing times (compared to FEM for example), especially in 3D problems. In this case, techniques such as the fast multipole method [228] have been introduced to accelerate the solution in large-scale problems; • A Fundamental Solution (FS) or Green's function, describing the behaviour of a point load in an infinite medium of the material properties is required as part of the kernel of the method. This can make the use of BEM infeasible in problems where a FS is not available; • calculation of the FS must be computationally efficient, which makes explicit FS formulations very desirable in this sense. Dynamic problems usually have implicit formulations, see [60,227,277] for instance, where the FS is expressed in a integral form by means of the Radon transform; • The BEM formulation requires the evaluation of weakly singular, strongly singular and sometimes hypersingular integrals which must be carefully treated. This can be done through a variety of methods, including singularity substraction, e.g. [100], or analytical regularisation, e.g. [91]; • Non-linear problems (e.g. material non-linearities) are difficult to model; The constitutive equations are given as with C ijkl and r ij denoting the elastic stiffness and the mechanical stresses, respectively, and where u i are the elastic displacements. The Einstein summation notation applies in Eqs. (24) and (25). The elastic tractions p ij are given by with n ¼ ðn 1 ; n 2 ; n 3 Þ being the outward unit normal to the boundary.
The time-harmonic equilibrium equations in the absence of body forces can be written as where t is the time and q is the mass density of the material. From Fig. 6, let X be a cracked domain with boundary C, which can be decomposed into two boundaries, an external boundary C c and an internal crack C crack ¼ C þ [ C À represented by two geometrically coincident crack surfaces.
The Dual BEM formulation for time-harmonic loading relies on two boundary integral equations (BIEs), one with respect to the displacements at a point n of the domain X c ij ðnÞu j ðn; tÞ þ Z C p Ã ij ðx; n; tÞu j ðx; tÞ dCðxÞ which follows from the differentiation of the displacement BIE and further substitution into the constitutive laws equation (for details see [90]). N r stands for the outward unit normal to the boundary at the collocation point n; c ij is the free term that comes from the Cauchy Principal Value integration of the strongly singular kernels p Ã ij ; u Ã ij and p Ã ij are the displacement and traction FS and d Ã rij and s Ã rij follow from derivation and substitution into Hooke's law of u Ã ij and p Ã ij , respectively. In most cases, the cracks are considered to be free of mechanical tractions. These boundaries conditions can be summarised as where the '?' and '-' superscripts represents the upper and lower crack surfaces, respectively. Eqs. (28) and (29) can be redefined in terms of the crack tip opening displacement (Du J ¼ u þ J À u À J ) in function of the crack-free boundary C c and one of the crack surfaces, say C þ In this latter equation, the free term has been set to unity due to the additional singularity arising from the coincidence of the two crack surfaces. The inconvenience of this approach is that the BEM formulation will now involve integrals including both strong singularities which require special treatment. Numerous hypersingular approaches have been developed, in particular to anisotropic materials under static [90,91,150,282] and time-harmonic [6,93,94,226,232,283,294] loadings. The use of a hypersingular formulation does not limit at all the crack shape, being valid for curved and branched cracks, for example. However, it is commonplace to make use of discontinuous boundary elements to ensure that all collocation points lie on the smooth surface within the body of an element; this is required to satisfy the Hölder continuity requirement of the hypersingular BIE.
As stated previously, the Stress Intensity Factors (SIF) are the measure of the stress amplification at the crack tip.
They are used extensively when estimating the structural life in a number of applications, from civil engineering structures to aerospace devices. Therefore, a precise calculation of this parameter is essential. The principal difficulty, faced throughout the development of BEM and FEM approaches for modelling LEFM problems, is the use of these discrete techniques to capture the singular stress solution. Traditional finite element piecewise polynomial shape functions are ineffective. We now describe some common approaches to obtain the SIFs: 1. Quarter-point: Developed by Henshell and Shaw [119] and Barsoum [19] for finite elements, it consists in moving the mid-side node of a quadratic boundary element from the centre to 1/4 of the element length from the crack tip. It was shown that the mapping between the element in real space and in the space of the intrinsic coordinates automatically captures the asymptotic displacement behaviour of 1= ffiffi r p present in the vicinity of the crack tip (refer to [231] for further explanations). 2. J-integral: Proposed by Rice [224], a path independent integral (assuming a non-curved crack) is used to evaluate the energy release rate due to the presence of the crack, where n 1 is the component of the outward unit normal vector in the x 1 direction, u i are the displacement and t i are the tractions. The term W ¼ 1 2 r ij e ij is the strain energy density. 3. Interaction integral: the J-integral can be decomposed into 3 parts [110,176] where J ð1Þ is the J-integral of the so-called principal state, which represents the energy release rate of the material; J ð2Þ is the J-integral of the auxiliary state, which depends on the displacements around the crack tip; M ð1;2Þ is the interaction integral containing terms of the principal and auxiliary state, and is defined as where A is the area inside the contour C j surrounding the crack tip, and W ð1;2Þ is given as ij Þ ð 36Þ Let us remark that the indices (1) and (2) correspond to the principal and auxiliary states, respectively.
The quarter-point approach allows a direct extrapolation of the SIF by using the crack opening displacement. The J-integral is more cumbersome numerically since the displacements and tractions at the closed path integral are part of the BEM domain and have to be evaluated first; however it is more accurate than the direct extrapolation.
Chen [47] has analysed mixed mode SIFs of anisotropic cracks in rocks with a definition of the J-integral for anisotropic materials and the relative displacements at the crack tip. In Ke et al. [133], the authors have suggested a methodology to obtain the fracture toughness of anisotropic rocks through experimental measurements of the elastic parameters and further comparison with a BEM code. In another work, Ke et al. [132] have proposed a crack propagation model for transversely isotropic rocks. Let us remark that all the previously mentioned works have used the Lekhnitskii formalism [145] in order to model the anisotropy of the material. The Lekhnitskii formalism is a polynomial analogy form of the matricial Stroh formalism.
Crack propagation problems have also been studied under the BEM framework. Portela et al. [214] used the maximum stress criterion as crack growth criteria in a dual BEM. Quasi-static 3D crack growth is analysed in [169].
Cohesive models have also been developed with the BEM: Oliveira and Leonel [195,196] have proposed a cohesive crack growth model, where the zone ahead of the crack tip is modelled as a fictitious crack model. This formulation gives rise to a volume integral, which must be regularised. The cohesive stresses are dependent on the crack tip opening displacement.
Yang and Ravi-Chandar [286] have proposed a cohesive model where the single-domain dual integral equations are used as an artifice to avoid the mathematical degeneration of the formulation imposed by the crack. In this case, the domain is divided in two sub-domains, where the crack is in the fictional domain division. Moreover, the cohesive zone is modelled as an elastic spring connecting both crack faces. Normal and tangential crack tip opening displacements are considered, and the crack growth is obtained from successive iterations of the non-linear system of equations, where the stiffness of the cohesive zone is taken into account.
Saleh and Aliabadi [233][234][235] and Aliabadi et al. [7] have studied the crack propagation problem in concrete using a fictitious crack tip zone. The cohesive zone is modelled with additional boundary elements at the fictitious crack tip that satisfy a softening cohesive law. A major drawback of this methodology is that the crack growth path has to be known a priori.

Fast Multipole Method (FMM)
The linear system formed in the BEM framework is much smaller than its equivalent with FEM formulation. However, the resulting matrix is full, not sparse like the FEM stiffness matrix, and this considerably increases the computational time required to solve a large problems. In 1985, Rokhlin [228] developed a method to reduce the complexity of solving the system of equations to OðnÞ instead of Oðn 3 Þ, where n is the number of unknowns. This technique was named the Fast Multipole Method (FMM), and generally involves using an iterative solver (such as GMRES [230]) to solve the linear system which comes from the discretisation of Eqs. (31) and (32) where C a is a boundary away from n. This change allows the collocation point n to be independent of the observation point x due to the introduction of a new point x e . Equation (39) has to be evaluated only once for different collocation points. The FMM applied in BEM can be described by the following steps [150]: 1. Discretise the boundary C; 2. Determine a tree structure of the elements. For example, in a 2D domain, define a square containing the entire boundary and call this square the cell of level 0. Then, divide the square into 4 equal cells and call them level 1. Repeat until each cell contains a predetermined number of elements (in Fig. 7, each cell has one element). Cells with no children cell are called leafs. For 3D cases, the same principle applies using cubic cells instead of square cells; 3. Compute the moments on all cells for all levels l ! 2 and trace the tree structure (shown in Fig. 8). The moment is the term from Eq. 6. Obtain the solution of Ax ¼ b. The iterative solver updates the unknown solution of x and goes to step 3 to evaluate the next matrix vector product Ax until the solution converges within a given tolerance.
The FMM has been used in 3D fracture mechanics problems as can be seen in [192,290], and some recent works on GPU can be found in [101,108,278]. The FMM is largely detailed in [149].

Adaptive Cross Approximation (ACA)
The Adaptive Cross Approximation (ACA) approach uses a different technique in order to reduce the complexity of the BEM with respect to the storage and operations. ACA uses the concept of hierarchical matrices introduced by Hackbusch [107], where a geometrically motivated partitioning into sub-blocks takes place, and each sub-block is classified as either admissible or inadmissible according to the separation of the node clusters within them.
The main idea is that admissible blocks are approximated by low-rank approximants formed as a series of outer products of row and column vectors, greatly accelerating the evaluation of the matrix vector product that lies within each iteration of an iterative solver. While the FMM deals with the analytical decomposition of the integral kernels, ACA can evaluate only some original matrix entries, or use a full pivoted form where all terms of matrix are calculated, to get an almost optimal approximation. The approximation of matrix A 2 C tÂS is given by where k is a low-rank compared to t and s. It is important to remark that the low-rank representation can only be found when the generating kernel function in the computational domain of A is asymptotically smooth. It has been shown in [20] that elliptic operators with constant coefficients have this property. In hierarchical matrices, the near and far fields have to be separated. The index sets I for row and J for columns so that elements far away will have indices with a large offset.
By means of a distance based hierarchical subdivision of I and J cluster trees T I and T J are created. In each step of this procedure, a new level of son clusters is inserted into the cluster trees. A son cluster is not further subdivided and is said to be a leaf if its size reaches a prescribed minimal size b min . Usually one of two different approaches is considered. First, a subdivision based on bounding boxes splits the domain into axis-parallel boxes which contain the son clusters. Alternatively, a subdivision based on principal component analysis splits the domain into well-balanced son clusters leading to a minimal cluster tree depth. Now, the hierarchical (H)-matrix structure is defined by the block cluster tree T IJ ¼ T I Â T J using the following admissibility criterion: minðdiamðtÞ; diamðsÞÞ gdistðt; sÞ, with the clusters t & T I ; s & T J , and the admissibility parameter 0 \ g \ 1. The diameter of the clusters t and s, and their distance, are obtained as A block b is said to be admissible if it satisfies this admissibility criterion. Otherwise, the admissibility is recursively verified for each son cluster, until the block becomes admissible or reaches the minimum size.
Finally, the ACA method idea is to split the matrix A 2 C tÂs into A ¼ S k þ R k , where S k is the rank k approximation and R k is the residuum which has to be minimised. We now present the ACA method itself: 1. Define k ¼ 0 where S 0 ¼ 0 and R 0 ¼ A and the first scalar pivot to be found is c 1 ¼ ðR 0 Þ À1 ij , and i, j are the row and column indices of the actual approximation step; 2. For each step t, obtain where the operators ðÞ i and ðÞ j indicate the i-th row and the j-th column vectors, respectively; 3. The next pivot c tþ1 is chosen to be the largest entry in modulus of the row ðR t Þ i or the column ðR t Þ j 4. The approximation stops when the following criterion holds The main advantage comparing to the FMM method is that ACA can be implemented directly into an existing BEM code. Moreover, due to its inherently parallel data structure, parallel programming can be readily implemented, increasing the computational efficiency. However, the original matrix A will not be entirely recovered.
Note that it is not necessary to build the whole matrix beforehand. The respective matrix entries can be computed on demand [20]. Working on the matrix entries has the advantage that the rank of the approximation can be chosen adaptively while kernel approximation requires an a priori choice.
A few recent works on ACA implementation can be found in [81,99]. Use of the method for problems in 3D elasticity can be found in [28,158] and the application of ACA in crack problems was shown for the first time in [137].

eXtended Finite Element Method (X-FEM)
The motivation that lay behind the development of X-FEM was to eliminate some of the deficiencies of standard FEM for crack modelling, most importantly the requirement for highly refined meshing around the crack tips and the mandatory remeshing for crack growth problems. The partition of unity [15] is a general approach that allows the enrichment of finite element approximation spaces so that the FEM has better convergence properties. In X-FEM, the partition of unity method allows element enrichment such that degrees of freedom (dofs) are added to represent discontinuous behaviour. In this framework, the mesh is independent from the discontinuities, so that cracks may now pass through elements rather than being constrained to propagate along elment edges. This gives the FEM much more flexibiility to model crack growth without remeshing.
Two types of enrichment function are applied in the X-FEM: the Heaviside enrichment function, responsible for characterising the displacement discontinuity across the crack surfaces, and a set of crack tip enrichment functions (CTEFs), responsible for capturing the displacements asymptotically around the crack tip. This latter presents complex behaviour, varying for different constitutive laws (see [12,79,193], for some different CTEF). In this sense, it is similar to the FS, necessary in BEM formulations. The displacement approximation u h ðxÞ with the partition of unity can be stated as [176] where N i is the standard finite element shape function associated with node i; u i is the vector of nodal dofs for classical finite elements, and a j and b a k are the added set of degrees of freedom that are associated with enriched basis functions, associated with the Heaviside function HðxÞ and the CTEF F a ðxÞ, respectively. N is the set of all nodes, N H is the set of all nodes lying on crack surfaces, and N CT is the set of all nodes belonging to elements touching a crack tip.
Since the CTEFs describe the displacements at the crack tip zone through the addition of several dofs, the stress concentration around the crack tip can be found more accurately with a significantly coarser mesh compared to the mesh used with standard FEM in a similar problem.
The presence of blending elements, which do not contain the crack but contain enriched nodes, is also important and has to be considered. These elements were analysed by Chessa et al. [48], and some studies have improved the convergence of blending elements (see [84], for instance). The X-FEM convergence rate can also be increased through the use of geometrical enrichment [142], where a number of elements around the crack tip receive the CTEF instead of a single element (this latter named topological enrichment). Figure 9 illustrates an arbitrary elastic body with a cohesive crack. The governing equations for a cohesive crack model are given by [178] where X is the domain, f b is the body force vector, f t is the external traction vector, r is the stress tensor, a is the load factor which controls the load increments, f c is the traction along the cohesive zone, and is a function of the crack opening Du. The discretisation of Eq. (50) yields where B is the finite element strain-displacement matrix, b is the vector containing the body forces and T c ðDuÞ is the cohesive softening law relating the crack surface normal traction f c to the crack opening Du. X-FEM has been widely used with cohesive models in the last few years. Some authors [45,51,175] have used a typical X-FEM formulation to model the cohesive crack, i.e., a Heaviside enrichment function is used to model the jump between the crack surfaces and a crack tip enrichment function is used to model the asymptotic behaviour at the crack tip.
Xiao and Karihaloo [284] have obtained the asymptotic displacement at the cohesive zone for isotropic materials based on the Williams expansion. The authors considered only the case where the crack is traction free and the crack is subject only to mode I. The obtained enrichment functions are where j is the Kolosov constant (for details refer to [284]), h is the crack orientation with respect to the x 1 axis, a 11 is a real constant and comes from the Williams expansion. In this case, Eq. (57) receives a new crack tip enrichment term, as in the X-FEM formulation for linear elastic fracture mechanics (see [110,175]). Zamani et al. [291] uses higherorder terms of the crack tip asymptotic field to obtain an enrichment function based on the Williams expansion. This approach has provided good results for isotropic materials, but it may not be the same for anisotropic materials. An alternative approach is to model the crack with Heaviside elements only [139,180,262,302]. Since there is no discontinuity at the crack tip, there are no SIFs at the Fig. 9 Elastic body with a cohesive crack crack tip, and therefore no crack tip enrichment function is required. The displacement field uðxÞ is given by where N i is the standard finite element shape function associated to node i and a j is the additional set of degrees of freedom associated with the Heaviside enrichment function H, defined as þ1 if it is evaluated above the crack or À1 if below the crack. The sets N and N H denote the standard and enriched nodes, respectively. The crack growth is modelled considering some rules, for example, if the level of stress at the crack tip is above the material tensile strength [178,262].
In [139], a 2D cohesive model for an isotropic material was presented, where both fluid and porous material interact. The pressure inside a crack is also modelled. The Heaviside enrichment function is employed, as well as a pressure enrichment function, which allows the continuity of steep gradients without enforcing this condition. The crack propagation criteria depends on the stress state at the crack tip. The fluid behaviour can retard crack initiation and propagation. A local change of the flow can be seen immediately after crack propagation. The deformation around the crack causes fluid to flow mostly from the crack itself since the crack permeability is much higher than the medium permeability. This flow from the crack to the crack tip causes closing of the crack. However, a delamination test has shown that if the stiffness and permeability are high, the fluid does not influence crack growth.
More methods for crack propagation in X-FEM can be found in [151,167,182,183,225] for brittle fracture and [168,179,291] for cohesive cracks.

Enriched BEM
The extended boundary element method (X-BEM) was first proposed by Simpson and Trevelyan [251] for fracture mechanics problems in isotropic materials. The main idea is to model the asymptotic behaviour of the displacements around the crack tips by introducing new degrees of freedom. The displacements u h ðxÞ are thus redefined as where N and N CT are the sets with non-enriched and enriched nodes, respectively, N i is the standard Lagrangian shape function associated with node i; u i is the vector of nodal degrees of freedom, and a a k represents the enriched basis functions which capture the asymptotic behaviour around the crack tips. In elastic materials, a a k is an 8component vector for two-dimensional problems, since only two nodal variables (u 1 , u 2 ) and four enrichment functions are needed to describe all the possible deformation states in the vicinity of the crack tip [110].
Hattori et al. [110] used the following anisotropic enrichment functions initially developed for the X-FEM where ; r is the distance between the crack tip and an arbitrary position, h is the orientation measured from a coordinate system centred at the crack tip, and A; B and l are obtained from the following eigenvalue problem ð60Þ (no sum on m) with Let us emphasise that the anisotropic enrichment functions can also be used for isotropic materials, since this is a degenerated case from anisotropic materials. For more details please refer to reference [110].
An enriched anisotropic dual BEM formulation using the above enrichment functions [111] for anisotropic materials is similar to the one used by Simpson and Trevelyan [251] for isotropic materials. The extended DBIE and the TBIE can be restated as where C c ¼ C þ [ C À stands for the crack surfaces C þ and C À . Only the element containing the crack tip receives the enrichment function. Strongly singular and hypersingular terms arise from the integration of the p Ã ij , d Ã rij and s Ã rij kernels and they may be regularised in the same way as shown in [92].

Meshless Method
Meshless (or meshfree) methods have been the subject of considerable interest in recent years as alternatives to Finite Element (FE) methods for solid mechanics problems. As the name suggests the main advantage of these methods is their (varying) lack of reliance on a division of the problem domain into a mesh of elements, thus removing issues associated with mesh generation and remeshing (perhaps required following large deformations which would lead to distorted and hence inaccurate elements). However, meshless methods tend to be more computationally expensive largely as a result of the lack of easily consultable connectivity information provided by a mesh, but also because there is greater complexity in the formation of shape functions. The most popular meshless methods for solid mechanics are the Element-free Galerkin method [24] and the Meshless Local Petrov-Galerkin method [14]. The key difference in both of these methods compared to FE methods is the use of shape functions based on a moving least squares (MLS) approximation [77]. Taking the Element-free Galerkin method as an example, the displacement approximation u h at location x is constructed as where N i are shape functions based on the MLS approximation (explained below), u i are nodal values and N is the set of nodes in support at location x, supports being defined using weighting functions centred at nodes. To build the shape functions we choose a polynomial basis, which can be of any order but low orders are usually used, e.g. a quadratic basis in 1D pðxÞ T ¼ 1; x; x 2 È É or in 2D pðxÞ T ¼ 1; x; y; x 2 ; xy; y 2 È É . At any location x we define the matrix P whose rows are the valued basis vectors p T for the nodes in support at x. A least squares minimisation procedure applied to the approximation at node locations and the nodal values then leads to the shape functions as where W is a diagonal matrix of values of node-centred weight functions at location x, which may be splines or exponential functions. Carrying this out in 2D or 2D is simply done via tensor products of the 1D case. Key points of difference as compared to FE methods should be clear, i.e. the shape function formation requires the inversion of a matrix, albeit a small matrix (dimension same size as the number of terms in the basis), the choice of nodal support size is crucial but not easy to define and the use of an MLS approximation contrasts with the interpolation used in FE methods and has the knock-on effect of making the imposition of essential boundary conditions more complicated. Overviews of the various meshless methods for solid mechanics can be found in a number of references [23,85,147,191].

Meshless Methods for Fracture
Ever since their initial development in the 1990s meshless methods have been applied to crack modelling [22,24,197], to dynamic fracture [26] and crack propagation [25].
The key advantage of meshless methods over standard FE methods for fracture is removal of the need to remesh during crack propagation. Another positive feature of meshless methods is that smooth stress results can be obtained for high stress gradients around crack tips [37] thus requiring less effort in postprocessing compared with the X-FEM. As with all numerical methods applied to fracture we have to find ways of dealing with the stress singularities at the crack tips and the discontinuities introduced by the crack surfaces. The former can be dealt in meshless methods by enriching the approximation space just as is done in X-FEM and other enriched methods, e.g. [27], based on the the partition of unity (PU) concept [16,166] where the jump discontinuity is included in the displacement approximation exactly as already laid out for X-FEM above in Eq. (49). ''Extrinsic'' techniques like this have more recently been developed into meshless ''cracked particle'' methods in a number of references [37,219,220,303]. Extrinsic enrichment like this can however lead to an ill-conditioned global stiffness matrix [21] as is the case with many other PU methods, due to the additional unknowns at nodes which do not correspond to the physical degrees of freedom [44]. The cracked particle methods are examples of smeared approaches to modelling cracks, i.e. the exact crack face/surface geometry is approximated, but this clashes with the requirement for an accurate description of the crack geometry since it governs the accuracy of field solution, and hence the crack growth magnitude and direction. Extrinsic approaches which attempt to improve on this have used piecewise triangular facets [37,64] which however suffer from discontinuous crack paths and requires user input to ''repair'' the mesh of facets. Greater promise lies in the use of a level set description of crack geometry combined with a meshless method [65,98,300,301] and an intrinsic rather than extrinsic model of the discontinuity of a crack. Using an intrinsic method in the EFGM there is also no problem of ill-conditioning in the stiffness matrix. Here the displacement jump can be introduced simply by modifying the nodal support via the weight function. A simple way to do this is directly to truncate the nodal support at a crack face. This is the visibility criterion, as shown in Fig. 10. The support of a node is restricted to areas of the domain visible from the node with the crack faces acting as an opaque barrier. If a line between a node and the point of interest intersects a crack, and if the crack tip is inside the support of that node, the node will have no influence on that point, i.e. r I between that point and the node is modified to infinity. (The visibility criterion corresponds to the use of the Heaviside function in the enriched trial functions used in X-FEM). An alternative to the visibility criterion is the diffraction method which works slightly differently as shown in Fig. 11.
The visibility criterion is simpler to implement, especially for 3D problems, but leads to spurious crack extension (thus impairing accuracy) while the diffraction method has no spurious crack extension problem but its implementation leads to high computational complexity especially in 3D or with multiple cracks.
Level sets offer a means accurately to represent crack surfaces and also to track surfaces as crack fronts propagate. The level set method (LSM) is a computational geometry technique for tracking interfaces applicable to many areas in science and engineering [243]. The LSM was first applied to crack modelling using the X-FEM in [98]. Instead of using an explicit representation of a crack, such as line segments in 2D and triangular facets in 3D, the LSM describes the surface implicitly by collecting points at the same distance to the crack into level sets. When the LSM is applied to fracture modelling, two orthogonal level sets, / and w are used: / measures the distance normal to the crack and w measures the distance tangential to the crack (see Fig. 12).
Hence we can fully define the geometry of the crack surface as /ðxÞ ¼ 0; wðxÞ 0 crack surface /ðxÞ ¼ 0; wðxÞ ¼ 0 crack front: As the crack propagates, the level sets are updated to the new crack surface using the procedures in [98] and the corrected update function for / in [65]. Recent work has led to the development of a fracture modelling method for 2D and 3D which uses intrinsic LS representations of cracks using a modified visibility criterion where the crack tips are tied to avoid the spurious propagation problem, and also incorporates enrichment to deal with the stress singularities [300,301].

Phase-Field
The development of the phase-field method provided an alternative formulation when dealing with different interface problems. A phase-field variable is introduced to consider the interface directly into the formulation. The phase-field formulation has been applied to different types of interface problems, including liquid-solid [13], liquidsolid-gas [165], electromagnetic wave propagation [260], analysis in crystal structures [1,58] and more recently in medicine [266], to enumerate some of the applications. The method has been successfully applied to fracture mechanics problems, where the crack is therefore modelled as a different interface in the domain. Figure 13 shows an example of a domain where the damage state is given by the interface parameter. The work of Francfort and Marigo [80] is the first in fracture mechanics to consider a variational formulation where a parameter assumes different values in order to capture the proper interface in the domain. An energy functional Eðu; CÞ, depending on the displacement field u and the crack surface C, is defined as [8,80] Eðu; where E d ðuÞ represents the elastic energy of the body, E s ðCÞ is the energy required to create the crack, considering Griffith's theory, w 0 is the elastic energy density and G c is the material fracture toughness. The work is further extended by [39] which applied a regularised form in order to allow the numerical treatment of the energy functional. The regularised energy functional E ðu; CÞ is given by where s is the phase-field variable, s ¼ 0 representing the undamaged state and s ¼ 1 standing for the fully broken/damaged state, with 0 s 1; [ 0 is a parameter designed to control the width of the transition zone set by the phase-field variable, and k is a small term depending on . The solution of Eq. (69) is found through the minimisation of E ðu; CÞ. To avoid the minimisation problem to be ill-posed, the small term k has been added to the formulation. For more details see [39].
The phase-field formulation has been modified through the years to be more general, consider more cases of interface interaction and different types of loading conditions to the problem. The work of Amor et al. [9] has considered the compression into the formulation, avoiding the interpenetration between crack surfaces. The proposed idea consisted in separating the elastic energy density according to the deviatoric and volumetric contributions.
A different phase-field formulation was proposed by [172,173], defined as a ''thermically consistent'' formulation. The regularised phase-field variable d is defined as 0 for the unbroken state and 1 for the fully broken state.
The stored energy w 0 ðeÞ of an undamaged solid is defined as [173] where w þ 0 ðeÞ is the energy due to tension and w À 0 ðeÞ is the energy due to compression. The positive and negative parts of the energy are given by the following decomposition of the strain tensor e ¼ X 3 a¼1 e a n a n a ð71Þ where a and n a are the principal strain and principal strain direction in the x a -axis, respectively. The standard quadratic energy storage function of an isotropic undamaged material is given by with k [ 0 and l [ 0 are elastic constants.
The phase-field model for fracture in elastic solids is given by where Div represents the divergent, Dd is the Laplacian of the phase-field, l is the width of the transition zone (where 0\d\1), k is a small artificial residual stiffness to prevent the full-degradation of the energy at the fully damaged state d ¼ 1, hxi À ¼ ðjxj À xÞ=2 is a ramp function, _ d is the evolution of the phase-field parameter.
A downside of the phase-field formulation is that it can result in unrealistic solutions. An example analysed by [8] consists of the case when the principal strains are negative, which is not considered in the model of [9] for instance. Nevertheless, a strongly non-linear strain relation is used, which requires higher computational charges as compared to [9].
A history-field variable was introduced in [172] in order to overcome some implementations issues which arose in [173]. Since the w þ 0 term determines the phase-field variable, we have Fig. 12 The level sets description of a crack surface in 3D (from [301]) Fig. 13 Phase-field domain Tðx; tÞ ¼ max Substituting Eq. (75) into (74) and applying a viscous regularisation, the evolution equation can be recast as where g [ 0 is a viscous parameter. The advantage of this new form is that the irreversibility of the crack phase-field evolution is put into a more general form, allowing loading/unloading conditions, besides allowing a better numerical treatment of the phase-field.
Crack branching effects are studied with phase-field in [117] for a 2D fracture problem. The instabilities are seen to appear at the critical crack speed of 0:48c s , where c s is the shear wave speed. It is worth to note that this relation is valid for perfect brittle materials only. Moreover, it was observed that, as the crack speed increases, the curvature of the area around the crack tip increases, splitting into two cracks when a critical value for crack speed is attained. In [118], a 3D study of crack branching stability is performed by means of fractographic patterns. The authors conclude that the instability is either restricted to a portion of the crack front or a quasi-2D branches.
A phase-field model is applied for damage evolution in composite materials in [29]. The evolution equation of the phase-field model was able to include difficult topological changes during damage evolution, such as void nucleation and crack branching and merging. Moreover, no meshing was required by the used phase-field model.
In [141], the formulation used in [39] is complemented by a Ginzburg-Landau type evolution equation, where an additional variable M is responsible for the crack propagation behaviour. If M is too small, the crack propagation may be delayed, while for sufficiently high values, the crack propagation is not affected by M. The FEM was coupled with the phase-field theory. This work was extended by [242] for dynamic brittle fracture.
Numerical aspects of the phase-field models used with finite differences, FEM and multipole expansion methods are discussed in [211].
More information about phase-field methods can be found in [38,50,217,242,253,275].

Configurational Force Method
Numerical implementations of brittle fracture propagation are relatively rare in the computational mechanics literature. One of the most promising numerical techniques developed within a conventional finite-element framework over the last decade is based on configurational forces. Within this setting, the most recent application of the configurational force methodology to the modelling of fracture is the work of Kaczmarczyk et al. [129], which focuses on large, hyperelastic, isotropic three-dimensional problems.
Kaczmarczyk et al.'s paper [129] is largely based on the work of Miehe and co-workers [103,170,171]. Miehe and Gürses [170] presented a two-dimensional large strain local variational formulation for brittle fracture with adaptive R-refinement, the simplification of this framework to small strain problems was presented by Miehe et al. [171]. The approach was extended to three-dimensions for the first time by Gürses and Miehe [103].
All of the works in this area are based on Eshelby [70,73] and Rice's [224] concept of material configurational forces acting on a crack tip singularity. A more general overview can be obtained from several sources [104,105,135,162,256]. Within this setting several local variational formulations have been proposed, for example see the works of [163,258], and fracture initiation defects of the classical Griffith-type brittle fracture overcome by global variational formulations [54,80]. Several researchers have numerically determined the material configurational forces at static fracture fronts [61,116,185,255]. Before the works of Miehe and co-workers [103,170,171], there were several other attempts towards the implementation of fracture propagation in the configurational mechanics context, including: Mueller and Maugin [186] within the conventional finite-element context, Larsson and Fagerström [74,143] in X-FEM and Heintz [115] within a discontinuous Galerkin (DG) setting. The framework has also recently been applied to materials with non-linear behaviour, see for example the works of Runesson et al. [229] and Tillberg and Larsson [265] on elasto-plasticity and Näser et al. [189,190] on time-dependent materials and the review by Ö zenç et al. [202]. In the following a configurational force approach to modelling fracture propagation is outlined based on the notation of Kaczmarczyk et al. [129].
The method can be cast within an Arbitrary Lagrangian-Eulerian (ALE) description of motion, where the deformation of the body is decoupled from the development of an advancing crack front (see Fig. 14). This approach requires the specification of three configurations: a reference state, B 0 ; and two current states: a material configuration, B t , containing the evolution of the crack surface; and a spatial configuration, X t , containing the physical deformation of the body. A conventional finite-deformation mapping, uðX; tÞ, connects the spatial, x, and material, X, configurations. Similarly the material and reference, N, frames are linked by a deformation mapping, Nðv; tÞ, that contains the structural change of the material. The crack surface is denoted as C 2 B t and the crack front, oC, as shown in Fig. 15.
Numerical Simulation of Fracking in Shale Rocks: Current State and Future Approaches 301 From the first law of thermodynamics, equilibrium of the crack front is governed by where _ W is the crack front velocity, A oC is a kinematic state variable that defines the current crack front direction and c is the surface energy. The configurational force at the crack front, oC, is given by where N is the normal to the surface encircling oC; R is the Eshelby stress tensor, L is the length oC; L n is the curve orthogonal to oC that defines the crack front encircling surface (as shown in Fig. 15). The Eshelby stress tensor, R, is defined as where WðFÞ is the free-energy function, F the deformation gradient, 1 is the second order identity tensor and P ¼ oWðFÞ=oF is the first Piola-Kirchhoff stress. As noted by Kaczmarczyk et al. [129], three possible solutions to Eq. (77): zero crack growth with _ W ¼ 0; force balance ðcA oC À G oC Þ ¼ 0; or that the crack front velocity is orthogonal to ðcA oC À G oC Þ. However, there is insufficient information in Eq. (77) to dictate the evolution of the crack front. Such an evolution law can be obtained by considering the second law of thermodynamics, supplemented by a material constitutive law and the principal of maximum energy dissipation.
Starting from a Griffith-type criterion for crack growth where g c is a material parameter controlling the critical threshold of energy release per unit area. Combining this with the principal of maximum dissipation, and through the application of Lagrange multipliers, it is possible to arrive at the condition that Therefore, the direction of crack propagation is constrained to be coincident with the configurational force direction. In addition, the configurational force approaches based on the work of Miehe and co-workers [103,129,170,171] utilise R-adaptive mesh alignment. This method aligns the propagating crack front with the direction of the configurational force by modifying the position of the node(s) attached to the element faces to be split.
In the work of Kaczmarczyk et al. [129], this fracture methodology was combined with a mesh quality control algorithm based on the work of Scherer et al. [240]. Within this, the nodal positions of the elements are modified based on a shape-based (volume to length) measure of element quality through the determination of a pseudo force vector. This pseudo force features in the discretised material nodal force equilibrium equation and is solved using a Newton-Raphson process. Note, that this modification to the discrete equilibrium equation only influences the stability of the solution and not the crack propagation criterion [129]. This mesh quality control procedure reduces the progressive degradation of the solution with fracture propagation.
Kaczmarczyk et al. [129] note that their approach could easily be extended to anisotropic materials. However, one limitation of the approach is that it is currently unable to capture non-smooth crack kinking [171]. Also, crack branching and multiple crack coalescence has yet to have been demonstrated, or even formulated. The discrete element method (DEM) has been initially developed for materials which have particle-like behaviour, such as soil and rocks [146]. The method was formally proposed by Cundall and Hart [53] and consisted of modelling of the interaction between elements using contact. This was later called the bonded-particle approach and is illustrated in Fig. 16 for two arbitrary bodies X 1 and X 2 having a normal contact stiffness K n . However, one of the main restrictions of this bonded-particle approach is that it did not allow rotations, and therefore does not consider momentum. To overcome this restriction, shear contact stiffness K s has been introduced to the formulation and can be seen in Fig. 17.
The DEM is characterised by the following properties [30,146,215]: • Finite displacements and rotations of the bodies is permitted, which includes complete detachment; • New contacts (or the absence thereof) are recognised automatically as the calculation progresses.
In practice, DEM is used in problems with a large number of elements, each element representing a body in contact. The formulation itself can be quite simplified compared to other discretisation methods, but it allows the simulation of complex behaviour, including material heterogeneities.
The DEM can be decomposed into several subclasses, which differ in some aspects such as the contact treatment, material models, number of interacting bodies, fracturing, and integration schemes [30].
In this framework, each element is a particular body which can be in contact with a number of surrounding elements. This implies that contact detection is one of the main problems that can arise, since missing a contact between elements can result in non representative behaviour of the model. Moreover, inspecting the elements for possible contact can require large amounts of computational processing time. The most common contact search algorithms are based on so-called body based search, where the vicinity of a given discrete element is searched for possible contact, and repeated after a number of iterations to check if the elements are still in contact. The Region Search algorithm [263] is an example of this kind of contact detection. Other contact detection algorithms use space search rather than a body search, and some examples are based on binary trees [30,36,208].
The next step is to obtain the contact forces. The calculation is usually performed with penalty based methods or Lagrange multiplier based methods. A review of contact algorithms evaluation can be found in [112].
The modelling of fracture using DEM has been mostly confined to element interfaces, where the breakage of the link between elements determines the appearance or propagation of the damage [30]. Particles can be bonded into clusters, where the bond stiffnesses are the equivalent to the continuum strain energy. Bond failure is assumed when the strength has exceed the maximum tension the bond can handle. Consistent breakage of the particle bonds define the fracture shape in the material. In [18,187], a combination of the FEM with DEM has been used to model fracture starting from a continuum representation of the finite elements, and as the damage appears it is restated in the discrete element framework. A multifracture FEM/ DEM scheme has been proposed by [212], where sliver elements arising from poor intra-element fracturing were avoided using local adaptive mesh refinement.
Discontinuous deformation analysis (DDA) is a variation of the DEM proposed originally by Shi and Goodman [246] to simulate the dynamics, kinematics, and elastic deformability of a system contacting rock blocks. While each block is treated separately in DEM, in DDA the total energy of the system is minimised in order to obtain a solution; a linear system of equations is obtained, resembling the finite element formulation. In fact, displacements and strains are taken as variables and the stiffness matrix of the model is assembled by differentiating several energy contributions including block strain energies, contacts between blocks, displacement constraints and external loads [146]. In the basic DDA implementation, each block is simply deformable as the strain and stress fields are constant over the entire block area, while the contacts are solved using regular contact algorithms that allow interpenetration between bodies [112]. To conclude, DDA is an  implicit formulation while DEM uses an explicit procedure to solve the equilibrium equations. DDA has been used extensively in rock mechanics applications, as can be seen in [113,114,155,267] for example. The influence of the bond parameters defined at the microscale and how they affect the response on the macroscale are analysed in detail in [49] for rock model analysis. It is shown that using a clumped-particle model, i.e. the particles rotate in a cluster instead of each particle being allowed to rotate, can reduce the limitations of the model, such as the overestimated ratio between tensile and compressive strengths, and the friction angles of the failure envelope.
A combined Lattice Boltzmann method (LBM) and DEM have been used to simulate fluid-particle interactions by [76]. The fluid field is solved by an extended 3D LBM with a turbulence model, while particle interactions are modelled using the DEM. Simulation results have matched experimental measurements.
There are available codes for the DEM, as the universal distinct element code (UDEC) [124], the ELFEN [218], the Yade [138] and Y-Geo [159]. More information on the discrete element framework can be found on [30,146] and some applications in [33,127,128].

Peridynamics
We will now introduce a new numerical method called peridynamics, which appears to be very promising for fracking problems. The main difference between the peridynamic theory and classical continuum mechanics is that the former is formulated using integral equations as opposed to derivatives of the displacement components. This feature allows damage initiation and propagation at multiple sites, with arbitrary paths inside the material, without resorting to special crack growth criteria. In the peridynamic theory, internal forces are expressed through non-local interactions between pairs of material points within a continuous body, and damage is a part of the constitutive model. Interfaces between dissimilar materials have their own properties, and damage can propagate when and where it is energetically favourable for it to do so.

Definitions
The peridynamics formulation was first developed by Silling [248], where he tried to overcome the limitation of current theories dealing with discontinuity, such as in fracture mechanics problems. The main argument was that the difficulty of existing theories was due to the presence of partial derivatives in the formulation to represent the displacement and forces, making necessary specific approaches to eliminate the singularities which would arise. Silling proposed a new formulation based on particular interactions as in molecular dynamics, but applied to continuum mechanics. The term peridynamics was adopted to describe this formulation, and it comes from the Greek roots for near and force. The pairwise interaction between two particles can be defined as [249] q € uðx; tÞ ¼ Z H fðuðx 0 ; tÞ À uðx; tÞ; x 0 À xÞ dV x 0 þ bðx; tÞ where q is the mass density, f is the pairwise force function that the particle x 0 exerts on the particle x; H is the neighbourhood of x, u is the displacement vector field, b is a prescribed force vector field (per unit volume). It is usual to adopt the relative position n of the two particles in the reference configuration as Analogously, the relative displacement g is stated as The current relative position can be easily given as g þ n. The function f must satisfy two conditions fðÀg; ÀnÞ ¼ Àf ðg; nÞ ð 85Þ which represents Newton's third law and enforces conservation of linear momentum, and ðn þ gÞ Â fðg; nÞ ¼ 0; 8g; n ð86Þ which assures conservation of angular momentum. The interaction between particles is defined as a bond, which in continuum mechanics could also be considered as a spring connecting two particles. This definition is fundamentally the difference between the classical theory and peridynamics, where the main idea is the direct contact between two particles. The area of influence of a particle is defined as the horizon d and is stated as Figure 18 illustrates the horizon d in an arbitrary body. Outside the horizon d, a particle has no influence on the other particles. For this reason, the peridynamics formulation is considered as a non-local model.
A material is microelastic if the pairwise function can be obtained through derivation of a scalar micropotential w such as The micropotential w is the energy present in a single bond (in terms of energy per unit volume squared). Thus, the local strain energy density is defined as where the factor 1/2 is present since each particle possesses half of the energy of the bond between them. If a material is microelastic, then every pair of particles x and x 0 is connected by a spring. The force in the spring depends only on the distance between the particles in the deformed configuration. Hence, there is a scalar functionŵ such that wðy; nÞ ¼ wðg; nÞ 8g; n; y ¼ jg þ nj ð 90Þ From Eqs. (88) and (90), the pairwise function f is restated as fðg; nÞ ¼ n þ g jn þ gj f ðjn þ gj; nÞ 8g; n ð91Þ with f ðy; nÞ ¼ oŵ oy ðy; nÞ 8y; g ð92Þ From Eqs. (82) and (91), the peridynamics model is fully defined for a non-linear microelastic material. However, a linearised theory of the peridynamics microelasticity can be defined as where C is the material's micromodulus function. It will be seen that the micromodulus has similar function to the material constitutive law. For more details, see reference [248]. Boundary conditions in peridynamics are not completely alike to the classical theory. Although the essential boundary condition is still present (displacements), there are no natural boundary conditions (tractions) in the peridynamics framework. Forces at the surface of a body must be applied as body forces b acting through the thickness of some layer under the surface. Usually, the thickness is taken to be the horizon d. The displacement boundary conditions also have to be imposed as a volume rather than a surface. For more details see [248].

Constitutive Modelling
We assume that the bond force f depends only on the bond stretch s, defined as As expected, s is positive only when the bond is under tension. Failure is introduced into the peridynamics model through breakage of the bonds connecting two particles over some stretching limit. Once a bond fails, it never becomes reconnected (i.e. no healing is considered). An example of a history dependent model is given by the prototype microelastic brittle (PMB) material, and is given by f ðyðtÞ; nÞ ¼ gðsðt; nÞÞlðt; nÞ ð 95Þ where gðsÞ ¼ cs; c is a constant and l is a history-dependent scalar-valued function, assuming either the values 0 or 1 according to In this case, s 0 is the critical stretch for bond failure. The local damage at a point can be defined as where x has been included as a reminder that the history model also depends on the position in the body. One can see that 0 u 1, 0 representing the undamaged state and 1 representing full break of all the bonds of a given particle to all other particles inside the horizon d. The broken bonds will eventually lead to some softening material response, since failed bonds cannot sustain any load. There are only two parameters that define the PMB material, the spring constant c and the critical stretch s 0 . Assuming g ¼ sn and substituting in Eq. (89), the local strain energy can be expressed as This relation must be identical to its equivalent in the classical theory, W ¼ 9ks 2 =2, where k represents the material bulk modulus [249]. The spring constant of the PMB material model is obtained as Now we describe the bond breakage formulation. Let the work G 0 necessary to break all the bonds per unit fracture area be given as The elements of Eq. (100) are depicted in Fig. 19. Equation (100) is the energy to break all points A, where 0 z d from the points B. After evaluation of the integrals we obtain

Anisotropic Materials in Peridynamics
The peridynamics formulation was initially presented for isotropic materials, in order to make some simplifying assumptions. It is expected then that the spring stiffness of the bonds does not vary over the direction of n. It was demonstrated in detail in [248] that for isotropic materials, the Poisson's ratio in the peridynamics formulations is constrained to take the constant value of 1/4. The constant Poisson's ratio is a consequence of the Cauchy relation for a solid composed of a lattice of points that interact only through a central force potential [153]. Refinements of the peridynamics theory can allow the dependence of strain energy density on local volume change in addition to two-particle interactions [154].
A composite material is formed by different materials, commonly a brittle and stiff material (fibre) embedded into a ductile one (matrix). In [201], the micromodulus C is redefined in order to accommodate the new variables arising from the material's anisotropy, including the fibre and matrix bonds for a laminate, and the shear and interlayer bonds present between two different laminates. However, in real composite materials, the fibre and matrix present properties vary significantly with the direction, which was not the case in this work. Instead, different isotropic materials were employed to form the composite fibre and matrix. In [122], the fracture in fibre-reinforced composites is tackled with more attention to the material modelling, where the differences between the fibre and matrix bonds are specifically defined. Moreover, the effect of arbitrary fibre orientation in the peridynamic model is taken into account, and it is shown that for a given particle x, the number of fibre bond particles within the horizon d can vary considerably, which leads to large variation of the strain energy density, the parameter which describes the bond stiffness. To consider this modelling issue, a semianalytical model was deduced for fibre orientation of 45 , and also for random fibre orientation.
A recent work [95] has deduced a peridynamic formulation for orthotropic media. The micromodulus C is defined in terms of the orientation of the angles u, as illustrated in Fig. 20. The dependency on the h orientation can be suppressed since the material properties do not change over h for a transversely isotropic material. After some mathematical manipulation, the new definition of the micromodulus is given as where A n0 represents constant coefficients and P m n are the associated Legendre functions of degree n and order m P m n ðcos uÞ ¼ ðÀ1Þ m 2 n n! ð1 À cos 2 uÞ m=2 d nþm dðcos uÞ nþm ðcos 2 u À 1Þ n ð103Þ Fig. 19 Fracture energy evaluation  (103) can be further simplified into Assuming cð0Þ ¼ c 1 and cðp=2Þ ¼ c 2 , it can be shown that Eq. (104) is also equivalent to where c 1 and c 2 are constants of the material model and are given by where C 11 ; C 22 ; C 16 and C 66 are elements of the constitutive matrix given in the Voigt notation. Note that an orthotropic material has only 2 independent material constants in the peridynamic model instead of the normal 4 independent constants. This restriction is used linked to the fact that a point is only able to interact to another one individually, while in the classical theory this condition does not apply (a disturbance in a continuous point will automatically induce some disturbance on the points around the body). This restriction on peridynamics theory has been addressed by Silling et al. [250] and will be detailed in the next section.
The critical bond stretch also depends on the direction of n and is given by where B n0 are constants and are detailed in [95]. The critical strain energy release rates for mode I crack propagation in the planes normal to the principal axes 1 (G Ic1 ) and 2 (G Ic2 ) can be obtained from the following relations After integration of Eqs. (110) and (111), the critical stretches s 01 and s 02 are given by The fracture behaviour of the material is fully defined by using the mode I energy release rates. Hence, mode II energies are not independent from mode I, which is another consequence of the bond-based peridynamic theory.
An important issue has been highlighted in [95,201], concerning the use of ''unbreakable'' bonds near to the regions where a traction boundary condition is applied. The possible reason for this would be crack initiation and propagation close to these regions, due to the high stresses that could be present. It is important to understand the physics of the analysed problem properly in order to use this type of assumption during a peridynamic simulation.

State-based Formulation
The peridynamics formulation assumes that any pair of particles interacts only through a central potential which is independent of all the other particles surrounding it. This oversimplification has led to some restrictions of the material's properties, such as the aforementioned fixed Poisson's ratio of 1/4 for isotropic materials. Also, the pairwise force is responsible for modelling the constitutive behaviour of the material, which is originally dependent on the stress tensor. To overcome this limitation, Silling et al. [250] have extended the peridynamics formulation to include vector states. The vector states allow us to consider not only a particle, but a group of particles in the peridynamics framework. Moreover, the direction of the vector states would not be conditioned to be in the same direction of the bond, as in the bond-based theory. This property is fundamental to consider truly anisotropic materials.
Let A be a vector state. Then, for any n 2 H, the value of Ahni is a vector in R 3 , where brackets indicate the vector on which a state operates. The set of all vector states is denoted V. The dot product of two vector states A and B is defined by The concept of a vector state is similar to a second order tensor in the classical theory, since both map vectors into vectors. Vector states may be neither linear nor continuous functions of n. The characteristics of the vector states are Numerical Simulation of Fracking in Shale Rocks: Current State and Future Approaches 307 listed in [250], and they imply the vector states mapping of H may not be smooth as in the usual peridynamic model, including the possibility of having a discontinuous surface.
In the state theory, the equation of motion (82) is redefined as qðxÞ € uðx; tÞ ¼ Z H fT½x; thx 0 À xiT½x 0 ; thx À x 0 ig dV x 0 þ bðx; tÞ ð 115Þ with T as the force vector state field, and square brackets denote that the variables are taken in the state vector framework.
To ensure balance of linear momentum, T must satisfy the following relation for any bounded body B Z B q € uðx; tÞ dV x ¼ The balance of angular momentum for a bounded body B is also required Z B yðx; tÞ Â ðq € uðx; tÞ À bðx; tÞÞ dV where yðx; tÞ ¼ x þ uðx; tÞ ð 118Þ The deformation vector state field is stated as Y½x; thni ¼ yðx þ n; tÞ À yðx; tÞ; 8x 2 B; n 2 H; t ! 0 ð119Þ The non-local deformation gradient for each individual node is given by where denotes the dyadic product of two vectors, and xðjnjÞ is a dimensionless weight function, used to increase the influence of the nodes closes to x. The use of this factor is still under study [280], but the assumption of xðjnjÞ ¼ 1 has been seen to provide good results. The discretisation of Eqs. (120) and (121) can be expressed as a Riemann sum as [280] Bðx j Þ ¼ X m n¼1 xðjx n À x j jÞððx n À x j Þ ðx n À x j ÞÞV n " # À1 Fðx j Þ ¼ X m n¼1 xðjx n À x j jÞðYhx n À x j i ðx n À x j ÞÞV n " # where m is the number of nodes with the horizon of node j.
x j must be connected to at least three other nodes in the system to ensure that Bðx j Þ will not be singular. In state vector peridynamics, there are two ways to determine how the force state depends on the deformation near a given point. The first consists of formulating a constitutive model in terms of the force vector T and the deformation state Y½x; t. In this case, the force state is defined as where W is the strain energy density and r indicates the Fréchet derivative, which is defined as any infinitesimal change in the deformation state dY resulting in a change of the strain energy density dW such as with H x being a sphere centred at the point x with radius equal to the horizon d. Note that the Fréchet derivative can be seen as an equivalent of the tensor gradient in classical theory.
The second approach to relating the force and deformation in a state vector framework is to adopt a stressstrain model as an intermediate step [42,280]. For a strain energy density WðFÞ, the stress tensor can be expressed as The force vector is redefined as [250] T ¼ rW ¼ oW oF rF ð127Þ After evaluation of the Fréchet derivative, the force vector can be defined explicitly as Thx 0 À xi ¼ xðjx 0 À xjÞ½rðFÞ t :B:ðx 0 À xÞ ð 128Þ The processing of mapping a stress tensor as a peridynamic force state is the inverse of the process of approximating the deformation state by a deformation gradient tensor. A peridynamic constitutive model that uses stress as an intermediate quantity results in general in bond forces which are not parallel to the deformed bonds. This type of modelling was called ''non-ordinary'' by Silling [250].

Numerical Discretisation
The discretisation of the peridynamics model is quite straightforward. Equation (82) can be rewritten as a finite sum where n is the time step and subscripts denote the node number, i.e., u n i ¼ uðx i ; t n Þ; V p is the volume of node p. Equation (129) is taken over all p nodes which satisfy jx p À u i j d. The grid spacing Dx is also an important parameter in the peridynamics discretisation. The discretised form of the linearised peridynamics model is given by The displacements u n i are obtained using an explicit central difference formulation, with Dt as the time step. Some studies of the stability of the numerical discretisation were described in [154,249]. It has been established that the time step must not exceed a certain value in order for the numerical discretisation to be stable. Moreover, the error associated with the discretisation depends on the time step with (OðDtÞ) and the grid spacing with (OðDx 2 Þ), Convergence in peridynamics is affected by two parameters: the grid spacing Dx and the horizon d. Reducing the horizon size for a fixed grid spacing will lead to the peridynamics solution approximating the solution using classical theory. However, fixing the horizon size while increasing the grid spacing will lead to the exact nonlocal solution for that particular horizon size [122]. As for domain discretisation methods, it is important to balance the size of the horizon so the damage features in the analysed body are properly considered, and the grid spacing should be sufficiently small for the results to converge to the non-local solution. Usually, it ranges from 1/3 to 1/5 of the size of the horizon.
In recent works, the peridynamics formulation is used conjointly with other discretisation methods, such as meshless formulation [249] and finite element formulation [154].
In [201], peridynamics is used only to obtain the prediction of failure of the composite material, where an FEM code is employed to solve the global problem. This type of combined approach is often necessary since the peridynamic formulation can demand significant computational power, a common problem in molecular dynamics simulations as well.

Conclusions and Prospective Work
We have seen that the hydraulic fracture problem presents several characteristics which makes its study complicated: the shale is not a homogeneous material, it is not isotropic, the nanoporosity may retard crack propagation as the fluid penetrates the rock, and a large fracture network has to be considered in which cracks develop at multiple length scales, all of which can greatly increase the computational solving time. Moreover, most current analytical and numerical methods do not take into account crack branching, a key factor in order to obtain a correct estimation of the extended fracture network.
The current fracture models for brittle rocks and fracking have been useful as a first step in offering a more realistic fracking model. There are of course other limitations attached to each of the numerical models discussed earlier: for instance, in cohesive models, the cohesive zone model is not a parameter to be found, so the crack propagation path is already known a priori. Most works on X-FEM and BEM models consider that the crack propagation path is unique; only recently have some works appeared considering crack branching [184,244,285].
Fracking models developed so far have not considered the full complexity of shale rocks. Ulm and co-workers [268][269][270] have established that shales are likely to be transversely isotropic materials, with the direction perpendicular to the bedding planes taken as the symmetry axis. This is mainly due to the deposition process. It was also stated that the shale anisotropy is due more to the interaction between the particles than the elastic behaviour of the shale components.
It was seen in [139] that the fluid penetrating the crack may retard crack propagation, so the material's porosity has to be taken into account in the numerical model.

Future Works
The main challenges researchers are facing with respect to the development of a new numerical formulation for modelling hydraulic fracture are: (1) the multiscale characteristic of the fracking in shale rocks, and (2) the requirement for the numerical method to deal with a large number of cracks simultaneously propagating and possibly branching.
For crack propagation and crack branching, the peridynamics formulation has been shown to have excellent results. A few issues have been raised about the method, such as how to choose the grid spacing (interval between particles) and the horizon (area of influence of a given particle). Even though an orthotropic formulation for 2D materials was developed by [95], there are some limitations over these formulations, since a direct bond force formulation is used. To overcome this limitation, a state-based formulation for anisotropic materials should be developed.
A multiscale model must be able to consider how a crack entering the RVE interacts with the voids that are present. Moreover, there must be a coupling between the Numerical Simulation of Fracking in Shale Rocks: Current State and Future Approaches 309 microscale (anisotropic) and the macroscale (transversely isotropic). The peridynamics formulation could be used to model the microscale, so the crack branching inside the RVE can be properly considered. Once the crack propagation path is obtained, another numerical method (X-FEM/X-BEM) can be employed to model the crack in the macroscale. Crack branching has already been considered in peridynamics in [106]. A comparison against experimental results of X-FEM, cohesive models and peridynamics in dynamic fracture is done in [5], where it is observed that the peridynamics model is able to capture the physical behaviour seen in experiments. A stochastic approach is likely to be the most useful way to model the extended fracture network, since the natural variability in geological conditions makes us unlikely to be able to obtain a deterministic model of the fracture system induced around any particular well. Moreover, the crack propagation obtained with the peridynamics formulation may change significantly if changes to the grid spacing or horizon size are made.