## Abstract

Hydraulic fractures are tensile cracks that propagate in pre-stressed solid media due to the injection of a viscous fluid. Developing numerical schemes to model the propagation of these fractures is particularly challenging due to the degenerate, hypersingular nature of the coupled integro-partial differential equations. These equations typically involve a singular free boundary whose velocity can only be determined by evaluating a distinguished limit. This review paper describes a class of numerical schemes that have been developed to use the multiscale asymptotic behaviour typically encountered near the fracture boundary as multiple physical processes compete to determine the evolution of the fracture. The fundamental concepts of locating the free boundary using the tip asymptotics and imposing the tip asymptotic behaviour in a weak form are illustrated in two quite different formulations of the governing equations. These formulations are the displacement discontinuity boundary integral method and the extended finite-element method. Practical issues are also discussed, including new models for proppant transport able to capture ‘tip screen-out’; efficient numerical schemes to solve the coupled nonlinear equations; and fast methods to solve resulting linear systems. Numerical examples are provided to illustrate the performance of the numerical schemes. We conclude the paper with open questions for further research.

This article is part of the themed issue ‘Energy and the subsurface’.

## 1. Introduction

Hydraulic fractures (HFs) are a class of brittle fractures that propagate in pre-stressed solid media due to the injection of a viscous fluid. On a geological scale, HF occur as kilometres-long vertical dikes that bring magma from deep underground chambers to the Earth’s surface [1–5]. Industrial applications of HF include: accelerating the waste remediation process [6], waste disposal [7] or preconditioning in rock mining [8,9]. One of the most common industrial applications is for oil and gas reservoir stimulation to enhance the recovery of hydrocarbons by the creation of permeable pathways [10].

In less than a decade, HF have become the focus of considerable public attention as new technology to develop multiple HF from deep horizontal wells has enabled the extraction of hydrocarbons from impermeable shale deposits. During this time the development of numerical models to capture HF propagation has become an area of intense research. These models can be classified into two broad categories: continuum models, which are the more classical approach to HF modelling that treat the solid and the fluid as continuous media; and discrete models developed more recently that treat the solid as a collection of particles in a lattice connected by springs, and the fluid is modelled either as a continuum or explicitly by discrete particles [11]. Hybrid discrete–continuum models [12] have also been developed. Recently, a class of phase field or smeared crack models has also been developed, in which fractures are not modelled explicitly but rather by distributed damage that is represented by a field variable [13–15]. Because of space constraints, in this paper, we will focus on continuum models that use explicit representations of fractures that propagate according to linear elastic fracture mechanics (LEFM).

Inspired by the seminal 1987 paper by Spense & Sharp [16] and the 1994 discovery of the so-called ‘Schlumberger Cambridge Research (SCR) viscous asymptote’ [17], the beginning of this millennium heralded intense research into the propagation of HF in the special case of plane strain. This research established the pivotal role of the asymptotic behaviour of the solution in the vicinity of the fracture tip [18–27] and the complex multiscale structure of the solution in which two or more physical processes compete at multiple length scales and on multiple time scales (see the 2016 review of this corpus of work in [28]).

However, progress on incorporating this multiscale structure into numerical models has lagged behind these developments. Despite significant progress made since the first algorithms were developed in the 1970s [29,30], many three-dimensional HF algorithms still use propagation criteria based on the LEFM of dry cracks (see the 2007 review of numerical models in [31]). This lag in development is due to the intrinsic complexity of the coupled HF equations, which involve a degenerate hypersingular integro-delay-PDE along with a singular free boundary problem. It was not clear how to develop numerical models of HF in three-dimensional elastic media with arbitrarily shaped boundaries able to capture the multiscale structure without prohibitively expensive and complex re-meshing. However, in 2008 a methodology was developed [32] that is able to capture a given asymptotic behaviour and use it to accurately locate the singular free boundary while using a relatively coarse structured mesh. The scope of the present paper aims to review the subsequent developments in this multiscale modelling effort.

In §2, we state the governing equations for the elastic medium, fluid flow, leak-off and proppant transport, boundary conditions, and tip asymptotics; in §3, we discuss discretization for the displacement discontinuity (DD) method and the extended finite-element method (XFEM), solving the coupled equations and locating the free boundary using the implicit level set method; in §4, we provide numerical examples and in §5, we provide interesting questions for further research.

## 2. Mathematical model

### (a) Assumptions

The equations governing the propagation of an HF in a reservoir have to account for the dominant physical processes taking place during the treatment, namely the deformation of the rock, the creation of new fracture surfaces, the flow of the fracturing fluid in the crack, the leak-off of the fracturing fluid into the reservoir and the flow of proppant particles within the fluid. In this section, we consider the governing equations for HF propagation that involve continuum descriptions of the fluid flow within the fracture and the deformation of the solid medium in which the fracture is advancing. For the solid medium, we choose to discuss two distinct formulations of the continuum equations namely: a DD boundary integral formulation and a formulation based on the XFEM. To that end we start with the general equilibrium equations and elastic constitutive equations, which form the basis for the XFEM model and describe how to arrive at an equivalent DD integral equation formulation. The integral equation formulation provides a succinct way to explore some of the unique properties of the coupled fluid flow–solid deformation equations via asymptotic analysis. While we start with a fairly general description, applicable to a heterogeneous elastic medium, we will ultimately, for the purposes of illustration, assume: (i) that the rock is homogeneous, linear elastic, brittle, permeable and of infinite extent whose stiffness is characterized by the Young’s modulus *E* and Poisson’s ratio *ν*, while the amount of energy required to break the rock is represented by the fracture toughness *K*_{Ic}—all assumed to be spatially constant; (ii) the fracturing fluid is assumed to be incompressible and Newtonian having a dynamic viscosity *μ*_{f}; (iii) Carter’s leak-off [33] is used to model loss of fluid to the reservoir and characterized by the coefficient *C*_{L}, which is also assumed to be constant; (iv) the minimum *in situ* confining stress field *σ*_{0}(*x*,*y*) is independent of *z* so that the fracture evolves in a plane orthogonal to the *z*-direction (figure 1*a*); (v) most of our discussion will focus on the situation in which the fluid front has coalesced with the crack front, which occurs rapidly under the high confinement conditions typically encountered in reservoir stimulation [34–36]; and (vi) gravity is neglected in the lubrication equation (we briefly relax this to discuss proppant transport and point to a recently developed proppant model that can include the effect of gravitational settling and autonomous proppant packing).

The solution of the HF problem involves determining the fracture aperture *w*(*x*,*y*,*t*), the fluid pressure *p*_{f}(*x*,*y*,*t*), the fluid flux and the position of the front , where *t* denotes the time and *x*,*y* are coordinates referenced to the injection point (figure 1*a*). The solution depends on the injection rate *Q*(*t*), the far-field compressive stress perpendicular to the fracture plane *σ*_{0}(*x*,*y*) (a known function of position) and the four material parameters *μ*′, *C*′, *E*′ and *K*′ defined as
2.1where *E*′ is the plane strain modulus and the alternate viscosity *μ*′, leak-off *C*′ and toughness *K*′ are introduced to keep the equations uncluttered by numerical factors.

### (b) Governing equations

#### (i) Elasticity

*Equilibrium equations and Hooke’s law.* Assuming that the fracture is propagating quasi-statically in an elastic medium, the equilibrium equations, Hooke’s law and the strain-displacement gradient relations can be expressed in the form
2.2where *u*_{i,j}=∂*u*_{i}/∂*x*_{j}, we have assumed summation over repeated indices, *σ*_{ij} is the stress tensor, *f*_{i} is the body force vector, *C*_{ijkl} is the tensor of elastic constants and *u*_{i} are the displacement components. It is convenient to represent a crack in terms of the components of the displacement jump (or DD) across the fracture surface in which are the displacements of the elastic body on the positive/negative sides of the normal *n*_{i} to the surface.

*Displacement discontinuity boundary integral formulation.* Assuming the elastic medium is homogeneous or piecewise homogeneous with a layered structure, it is possible to reduce the PDE relating the displacement and stress field in the solid (2.2) to an equivalent integral equation that relates the crack opening displacement *w*=[*u*_{i}]*n*_{i} to the fluid pressure *p*_{f}
2.3where denotes the fracture footprint enclosed by the crack front , *σ*_{0}(*x*,*y*) is the prescribed *¡textit¿in situ*¡/textit¿ stress field and *p* is the net pressure. The Green’s function *G*(*x*,*y*;*x*′,*y*′) is obtained by solving (2.2) subject to an appropriate singular forcing [37] or equivalent body force [38]. For example, the point DD vector [*u*_{i}]=Δ*uν*_{i} located at (*x*_{s},*y*_{s},*z*_{s}) across a surface with a unit normal *n*_{j} can be obtained by solving (2.2) subject to the equivalent body force [39]:
2.4Using this body force representation and a double Fourier transform the Green’s function *G*(*x*,*y*;*x*′,*y*′) appropriate for modelling planar fractures in a layered elastic medium on an elastic half-space has been constructed [40–42]. If the elastic medium is assumed to be homogeneous and infinite then the Green’s function reduces to [43,44]

*XFEM variational formulation.* The primary utility of the XFEM for HF problems is to solve the equilibrium equations for a fractured elastic medium. Two properties of the HF models substantially affect the XFEM formulation that is required. In §2b(iii), we will see that if the fluid and fracture fronts are distinct (i.e. there is a finite fluid lag), the fluid pressure is finite while if they coalesce the pressure at the tip tends to . In addition, a local analysis of (2.3) shows that the action of the integral operator on *w* is similar to taking a derivative so that the net pressure field has the same regularity as the derivative of *w*. For this reason, the XFEM will yield accurate results (*O*(*h*^{2}) with the correct enrichment [45]) if the pressures are specified and the XFEM is used to solve for the corresponding fracture widths—i.e. determining the so-called Neumann-to-Dirichlet (ND) map, which is equivalent to numerical integration. Conversely, the XFEM will yield poor results (at best *O*(*h*) accuracy) if the fracture widths are specified and the XFEM is used to determine the corresponding net pressure field—i.e. determining the so-called Dirichlet-to-Neumann (DN) map, which is equivalent to numerical differentiation. Thus, if there is a finite fluid lag then, as the fluid pressure field is finite, it is possible to use an XFEM formulation that yields the ND map. However, if there is no lag, so the fluid and fracture fronts coalesce, then it is necessary to use a mixed variational formulation for the XFEM in which it is possible to specify the fracture widths near the fracture front at which the pressure field becomes infinite while away from we specify the pressure field and solve for the fracture width field. This type of geometric decomposition of the fractured region into a so-called ‘channel region’ away from the tip and ‘tip region’ between the channel region and the fracture boundary has been used to great advantage in a number of contexts that yield highly accurate numerical schemes. We will report on a number of such schemes throughout this paper. We note that the formulation we give here is for an HF in a three-dimensional elastic medium. The formulation is virtually the same for a two-dimensional plane strain HF in which is just a curve and the boundary curve reduces to the two points at the tips of the fracture.

We now describe the variational formulation [46] of (2.2) for a fractured elastic medium that is suitable for devising an XFEM scheme in which the fluid front and the fracture front coalesce. Within the tip region the fracture width *w*_{t} is prescribed while in the channel region the normal stress *σ*_{n} is prescribed, while the shear stresses *σ*_{sk} are prescribed for the complete fracture. Assume that the fracture surface is entirely contained within the domain having an outer boundary *Γ*. Assume that is discretized into a mesh of non-overlapping finite elements *e* each of which occupies the region such that . Let be the subset of elements that intersect , i.e. . Thus is partitioned into the two volumetric regions and and .

Let *U* be a finite dimensional subspace of the Sobolev space in two dimensions or *H*^{1}×*H*^{1}×*H*^{1} in three dimensions, then assuming Dirichlet boundary conditions on the remote boundary
the displacement in is approximated by elements of the trial space while variations are taken from elements of the test space . For each element , the stress is introduced as an auxiliary tensor for which Hooke’s law from (2.2) is weakly imposed. We approximate ** σ** by using trial and test functions from the symmetric tensor function space where

*H*

^{−1,h}is a finite-dimensional subspace of the space of functions square-integrable in each element in and discontinuous at element edges and normal to the crack .

The weak formulation [46–48] of the elasticity problem (2.2) seeks to find such that for all : 2.5where and {⋅} defines the averaging operator . Finally, we observe that the XFEM formulation, appropriate for the ND map, involves only the first three terms shown on the first line of (2.5) in which in the first integral is replaced by and in the second integral is replaced by .

#### (ii) Fluid transport

*Poiseuille flow and fluid conservation.* The lubrication equations are obtained by combining the average flux for Poiseuille flow
2.6and the continuity equation
2.7to yield the Reynolds equation
2.8where *w* is the fracture width, *p*_{f} is the fluid pressure, *μ*′ and *C*′ are the scaled dynamic viscosity and Carter leak-off coefficients, respectively, which are defined in (2.1), *H* is the Heaviside function, *δ*(*x*,*y*) is the Dirac delta function, and *t*_{0}(*x*,*y*) represents the trigger-time at which the unknown fracture front passes the point (*x*,*y*). We observe that including the history term *t*_{0}(*x*,*y*), which is not known *a priori*, turns (2.8) into a delay partial differential equation.

*Leak-off, poroelasticity and Carter’s model.* When fluid leaks from an HF into the porous host rock an impermeable *filter cake* layer forms as polymer molecules are deposited by the leaking fluid. Thereafter, an *invaded zone* forms that penetrates the filter cake, which releases filtrate that is transported along with pore fluid through the porous medium. This fluid leak-off causes an increase of pore pressure resulting in dilation of the rock that can, in turn, potentially inhibit the propagation of the fracture. Biot’s theory of linear poroelasticity [49] has been used [50,51] to investigate the feedback coupling between the increased pore pressure due to the infiltrating fluid and concomitant rock dilation. The linear poroelasticity equations [51] involve a Navier type equation for the displacement field *u*_{i} and a diffusion equation for the pore pressure . These equations are coupled by a pore pressure gradient forcing term in the Navier equation and a volumetric strain rate term in the diffusion equation each multiplied by the Biot constant. Detournay & Cheng [50] have demonstrated that, under hydraulic boundary conditions involving the pore pressure, the fluid exchanged between an HF and the host rock that is calculated using poroelasticity theory is nearly identical to that computed using uncoupled diffusion.

Carter’s leak-off model, which is represented by the *C*′ term in (2.8), makes the further assumption that the pressure diffusion is one dimensional. In addition, Carter’s model assumes that the hydraulic load that drives the leak-off process is approximately constant and lumped into the coefficient *C*′, where *p*_{0} is the reservoir pore pressure. This assumption can be justified for high confinement reservoirs in which . Thus, the hydraulic load is much larger than the net pressure and is approximately constant, i.e. Δ*p*≈*σ*_{0}−*p*_{0}. In addition, Carter’s model can also be used [18] to capture under a single lumped coefficient *C*′ the following three processes taken in series: filter cake formation, invaded zone growth and one-dimensional pressure diffusion into the reservoir.

While the characteristic fracture length *L* is much larger than the diffusion length scale , where *c* is the diffusion coefficient and *T* is the treatment time, the pore pressure diffusion around the fracture is one dimensional. Indeed, Carter’s model, which is considered in this paper, has been shown to be appropriate for HF treatments of hydrocarbon reservoirs [10], which take place on the order of hours. Moreover, in [52] it was shown that, provided *σ*_{0}≫*p*_{0} then a Carter-like leak-off model is valid even if the diffusion is not one dimensional. However, if , then the diffusion pattern is three dimensional. For example, models of water-flooding treatments that take place over months or years must account for large-scale three-dimensional diffusion. Kovalyshen [52] exploited the decoupling of the pure elastic deformation due to the mechanical loading of the crack from pore pressure diffusion to investigate the evolution of radially symmetric fluid-driven fractures that involve fully developed three-dimensional diffusion. In this paper, to develop the asymptotic analysis required to locate the free boundary, we exploit the fact that the tip region of a finite HF can be represented by a semi-infinite fracture, which was shown [52] only to be true for a permeable medium if . However, when three-dimensional diffusion is fully developed around an HF, a new approach will need to be developed to analyse the tip region.

*Proppant transport.* Proppant transport has typically been modelled [31] by representing the proppant-fluid mixture as a slurry. To model the proppant transport the evolution of the volumetric concentration *ϕ* is captured by the advective mass conservation equation
2.9The proppant velocity in (2.9) is determined from the slurry velocity and the slip velocity via
2.10in which the slip velocity is typically assumed to be due to the gravitational settling via the Stokes settling velocity [31]; and finally, the slurry viscosity used in (2.6) is assumed [31] to increase from the nominal clear fluid Newtonian viscosity *μ*_{f} as the proppant concentration *ϕ* increases until it approaches its maximum value *ϕ*_{m}=0.585 according to the relation
2.11While the model (2.9)–(2.11) is able to capture proppant transport for moderate concentrations and Stokesian settling, as the adjusted viscosity becomes infinite so that some form of cut-off viscosity is typically imposed. Physically this situation occurs when the proppant starts to pack so the advective-Poiseuille model is incorrect as the appropriate description of the fluid motion in this case is that of a clear fluid infiltrating a porous medium as described by Darcy’s law.

Based on recent experiments [53] and using appropriate averaging over the fracture width there has been considerable progress [54,55] in the formulation of an effective lubrication theory that provides a more realistic proppant transport model. Without making any assumptions about the slip velocity the effective model [55] was derived using a matched asymptotic analysis featuring two-dimensional flow and a boundary layer, and is able to capture gravitational settling and can autonomously capture the transition from Poiseuille flow for moderate *ϕ*/*ϕ*_{m}<1 to Darcy flow as when the particles start to pack. In this model the flux vector in (2.7) is replaced by the slurry flux while the advection equation in this case involves the proppant flux
2.12where is the normalized, width-averaged proppant concentration, and the slurry and proppant fluxes and are given by
and
where *μ*′ is the scaled clear fluid viscosity, is the fluid pressure corrected by the hydrostatic pressure, *ρ*^{p}−*ρ*^{f} is the difference between the proppant and slurry mass densities, *g* is the gravitational acceleration, *a* is the proppant particle radius, *B* is a ‘blocking function’ introduced to capture bridging when *w*=*O*(*a*) and , and are given functions. We observe that can be decomposed into two terms, the first *Q*^{s}, when multiplied by 1/*μ*′, acts like the reciprocal of the effective fluid viscosity analogous to (2.11) and tends to zero as , while the second term tends to , which yields the limiting Darcy flux vector , where *D* represents the permeability of packed particles. Thus, the variation of with enables the slurry flux to transition autonomously from Poiseuille flow to Darcy flow; represents the proppant motion caused by slurry flow and captures the gravitational settling of particles.

This proppant transport formulation has demonstrated the capacity to capture proppant packing that occurs in ‘tip-screenout’ and gravitational settling in implementations involving fully coupled plane strain two-dimensional and pseudo-three-dimensional [56–58 HF models [59] and, more recently, a three-dimensional HF simulator [60]. The numerical solution of the advection equation (2.12) presents considerable challenges particularly when in which limit the proppant starts to pack leading to shock waves that emanate from the packed regions in the upstream direction. This represents a particularly challenging class of shock capture problems in which the flux function is both spatially varying and a nonlinear function of the solution of a non-local operator equation through the coupling to the elasticity-fluid-flow equations (2.3)–(2.7). Because of the complex shock behaviour exhibited by (2.12) delicate shock-capturing disretizations are required that use explicit time stepping, while the coupled elasticity-fluid-flow equations (2.3)–(2.7) are extremely stiff so that explicit schemes require [31,61] Courant-Friedrichs–Lewy conditions of the form Δ*t*=*O*(Δ*x*^{3}). Thus, the coupled elasticity-fluid-flow equations (2.3)–(2.7) are typically solved implicitly allowing for relatively large time steps while explicit shock capturing schemes are used to solve (2.12) on sub-time steps of the implicit scheme. In order to avoid having to use these delicate shock capturing schemes, a novel Lagrangian formulation of (2.12) has been developed [62], which involves following packets of particles and shows superior performance to the grid-based shock-capturing schemes.

#### (iii) Boundary and propagation conditions

For this complex nonlinear problem, for which there exist multiple equilibrating and volume-conserving fluid pressure and width fields, imposing the appropriate boundary and propagation conditions is key to selecting the desired global solution. Since the development of the first HF models in the mid-1950s many analytic models for plane strain and radial fractures were based on *ad hoc* assumptions while many numerical models inherited propagation conditions based on the LEFM of dry cracks. It is only with the more recent tip asymptotic analysis that the importance of the appropriate tip boundary conditions has been fully appreciated [36].

*Finite lag between fluid and fracture fronts.* When the fluid front is entirely contained within but distinct from the fracture front so that there is a finite lag between the two fronts, the following boundary and propagation conditions apply:
2.13and
2.14In each of (2.13) and (2.14), the first conditions prescribe boundary conditions on the fluid pressure and width fields, while the second conditions in (2.13) and (2.14) provide propagation conditions for the moving fluid front and fracture front , respectively. In particular, the second condition in (2.13) is the classic Stefan condition, typical in free-boundary problems, in which the front velocity is equated to the average fluid velocity given by . The second condition in (2.14) imposes the condition that at all points of that are propagating, i.e. for which , the mode I stress intensity factor *K*_{I} matches the fracture toughness *K*_{Ic}, while those points of for which the stress intensity factor falls below the fracture toughness the boundary does not move, i.e. . Here is the velocity of the crack front and, for simplicity, we exclude the situation of fracture recession. From a computational point of view, even though a finite lag involves tracking the two distinct free boundaries and , a finite lag problem is simpler to model than the vanishing lag situation in which the two fronts coalesce to yield a singular free boundary problem.

*Coalescent fluid and fracture fronts.* Most HF treatments take place in environments in which the confining stress *σ*_{0} is large compared with a characteristic driving stress in the problem in which case the lag between the fluid and fracture fronts becomes negligible. Rather than trying to capture this vanishingly small lag with an extremely fine mesh it is preferable and more efficient to derive the appropriate boundary conditions that hold in this limiting situation. We consider the limit in which the two fronts coalesce in the simplest case of an impermeable medium *C*′=0 (for further details, see [36]). Assume that the confinement *σ*_{0} is large compared with a reference stress so that *σ*_{0}/(*μ*′*V* *E*^{′2}/*K*^{′2})≫1 then the fluid lag λ decays exponentially as this ratio increases [34] and vanishes as (here ). In this limit for , the Stefan condition , combined with , and the fact that the crack front velocity is finite, implies that and . Thus, the zero pressure boundary condition at the fluid front (2.13) is replaced by a zero flux boundary condition at the crack front and from (2.3), it follows that the net pressure . Thus the boundary conditions that apply for coalescent fluid and fracture fronts are
2.15We observe that as the first condition in (2.15) implies that is a level set for *w* and the limiting process implies that it is also a level set for , thus (2.6) implies that the flux vector is orthogonal to the fracture front . LEFM [63] implies that the width field in the vicinity of the crack front has the asymptotic behaviour
2.16and
2.17where *s* denotes the distance from the crack front (with the *s*-axis directed inwards as shown in figure 1*a*), and *K*_{I} is the unknown stress intensity factor. Thus, the stress intensity factor is assumed to be in limit equilibrium with the fracture toughness (i.e. *K*_{I}=*K*_{Ic}) if it propagates so that . However, if the fracture does not propagate, i.e. , the stress intensity factor is below its critical value *K*_{I}<*K*_{Ic}, which implies partial closure of the crack in the region near the tip. Complete fracture closure *K*_{I}=0 implies that the fracture starts to recede.

We observe that in the limit , the asymptotic relation yields the first two conditions in (2.15). As the front velocity for an impermeable medium is equal to the average fluid velocity as , it follows that 2.18which shows that is the limit of an indeterminate form, which is difficult to evaluate numerically. For a permeable medium, in the region of dominant leak-off the average fluid velocity tends to infinity as , and the crack front velocity needs to be determined by detailed asymptotic analysis [36]. Thus established front location algorithms [64–66], which are based on knowing the front velocity explicitly, cannot be applied directly to this singular free-boundary problem.

#### (iv) Tip asymptotics, vertex solutions and generalized asymptotes

A propagating HF with zero lag is governed by two competing dissipative processes: one associated with viscous energy losses and the other with the energy expended in breaking the rock, which is characterized by the fracture toughness; and two competing fluid-balance components: fluid storage in the fracture and storage in the surrounding rock—associated with leak-off. Limiting regimes can be identified in which just one of the dissipative processes and one of the storage mechanisms is dominant. Associated with each of these limiting regimes are so-called vertex solutions with reference to vertices of a multi-process phase space within which the dissipative and storage mechanisms compete. For example, for the storage-viscosity regime known as the *m*-vertex, leak-off is negligible compared with the fluid stored in the fracture and the energy expended breaking the rock is negligible compared with viscous dissipation, so that *K*′≈0≈*C*′. An HF having radial symmetry and for which *K*′≠0≠*C*′ will change its regime of propagation over time during which multiple length and time scales are important (see [28,32,67] for further details). The key tool to characterize this multiprocess, multiscale environment is an analysis of the tip asymptotics and the so-called semi-infinite fracture problem.

*Tip asymptotics and semi-infinite crack problem.* Provided the fracture front is sufficiently smooth that a finite radius of curvature can be defined at each point along its perimeter and assuming a homogeneous elastic medium it can be shown [32] that the governing equations (2.2), (2.8) and (2.16) for the aperture *w* and net pressure *p* in the vicinity of the fracture front reduce to those for a semi-infinite fluid-driven fracture steadily propagating in a state of plane strain at a constant velocity *V* and characterized by zero lag (see [22,34]), namely:
2.19where *s* represents a coordinate located on the moving tip and pointing toward the interior of the fracture. The consequence of this reduction is profound in that the tip asymptotic solution for a *finite fracture* at any time is given by the solution of the stationary semi-infinite crack in a state of plane strain whose constant tip velocity corresponds to the *current* propagation speed of the finite fracture. The tip solution is thus autonomous. Note that the spatial variation of the far-field stress can be ignored when viewed at the tip scale, unless the stress field is discontinuous (in which case, the tip solution outlined here is not relevant).

*Vertex solutions.* Assuming the fracture width is given by a power law and recognizing that in this case the integral operator (2.19) reduces to a Mellin transform, it follows from Oberhettinger [68] that
2.20Thus for this range of values of *α* the power law is an eigenfunction of the integral operator in (2.20).

This power-law solution pair can be used to identify the three limiting regimes of propagation also known as vertex solutions, namely: the viscosity-dominated solution, the leak-off-dominated solution and the toughness-dominated solution. Choosing *α*≠1/2 and substituting the power-law behaviour and the corresponding power-law behaviour for given by (2.20) into the lubrication equation in (2.19*a*), we can obtain two different vertex solutions by identifying two distinct dominant balances. Choosing the dominant balance between the term on left side and the first term on the right, we find that *α*=2/3 and matching coefficients yields the so-called *m*-vertex solution [17]
2.21Choosing a dominant balance between the term on the left and the second term on the right of (2.19*a*), we obtain that *α*=5/8 and matching coefficients we obtain the so-called -vertex solution [69]
2.22Finally, if *α*=1/2 then the toughness asymptote (2.16) is the so-called *k*-vertex solution *w*_{k} and matching the term on the left of (2.19*a*) with the first term on the right we find that .

*Generalized asymptotes.* From the vertex solutions described above, we see that for *s*≪1 there is a hierarchy of powers
2.23that establishes the region moving away from the fracture tip that will be occupied by the different solutions. For example closest to the tip the toughness solution (provided *K*′>0) will always be dominant. Then provided *C*′>0, as one moves further from the tip, the solution will transition to the leak-off asymptote characterized by the 5/8 power law. The length scale at which this transition occurs is determined by the relative magnitudes of the coefficients of the two solutions. Finally, the leak-off solution will transition to the viscous solution characterized by the 2/3 power, which prevails farthest from the tip. Thus, we can see that there are multiple physical processes present that can be active at vastly different length scales. If only two of these processes are present then the same ranking in (2.23) determines the regions occupied by the two processes relative to the tip.

Garagash [70] was the first to consider a two process solution involving the transition from the toughness to the viscous dominated regimes (*C*′=0) by solving the semi-infinite crack problem (2.19) with the boundary condition at infinity . The numerical solution to this connection problem involves the generalized asymptote shown in figure 3*a* in terms of which the asymptotic solution can be expressed in the form
2.24where ℓ_{mk} is the transition length scale. The *m* and *k* vertex solutions as well as the *m*–*k* transition solution have all been observed experimentally for penny-shaped fractures [28,71].

The corresponding viscous to leak-off transition (*K*′=0) involving the transition length scale has been considered by Adachi & Detournay [20]. Having established the two-process ‘edge solutions’ Garagash *et al.* [72] have provided the numerical solution for the generalized tip asymptotics for an HF in which all three processes are active (*K*′>0, *μ*′>0,*C*′>0). This generalized asymptote involves multiple length scales, which are characterized by the dimensionless parameter . Recently, a simpler non-singular formulation of the hypersingular integral equation (2.19) has made it possible to solve the multiscale problem using standard numerical techniques and yields an approximate formulation in terms of a separable first order ordinary differential equation (ODE) [73]. The closed form solution of this ODE yields not only all the vertex solutions but also the following approximate solution [73,74] in terms of a given function *g*_{δ}:
2.25that is able to capture the complete multiscale solution within a fraction of a per cent. This generalized asymptote has been used to model three-process multiscale planar [74] and multi-planar [75] HFs.

## 3. Discretization, coupled equations and the multiscale implicit level set algorithm scheme to locate the free boundary

In this section, we briefly describe the discretization processes on structured meshes appropriate for the DD and XFEM formulations described above. As the resulting coupled evolution equations are extremely stiff (for which explicit schemes require time-step restrictions Δ*t*≤*O*(Δ*x*^{3}) [31,61]) L-stable backward Euler implicit time stepping is often used to evolve the solution. Hence we also describe iterative schemes and preconditioners that have been developed to solve the fully coupled equations that need to be solved at each time step. Having established the importance of the tip asymptotic behaviour in governing the propagation of HF, we discuss specialized numerical schemes that have been developed to locate the free boundary and to accommodate this multiple length scale behaviour that can span 5–20 decades without having to resort to mesh refinement. In a recent comparative study [76], this class of schemes was shown to provide accurate solutions on a relatively coarse fixed structured mesh.

### (a) Discretization

#### (i) Displacement discontinuity formulation for planar fractures

We assume [32] that the fracture will grow within a rectangular region that has been tessellated into a fixed uniform rectangular mesh with dimensions Δ*x* and Δ*y* in the two coordinate directions (similar to the fracture shown in figure 1*a*). The fracture footprint is then covered by rectangular elements such that Constant DD elements are used for the elasticity computations [43] along with collocation at element centres to yield a fully populated linear system relating the element net pressures *p* to the element widths *w*
3.1in which, for a homogeneous elastic medium, the discrete Green’s function *G* assumes the convolution form
3.2Consistent with the constant DD discretization with collocation at element centres, the lubrication equation (2.8) is discretized via the finite-volume method to yield a five-node finite-difference stencil [32,42] and the backward Euler scheme to obtain
3.3where *S* is a vector of source/sink terms and *A*(*w*) is the difference operator defined by
3.4and are the fluxes along the vertical edges of the element defined as , where are the corresponding edge widths defined by . The fluxes and widths along the horizontal edges of the element are defined analogously. Zero-flux boundary conditions are implemented in tip elements by removing those terms associated with the element faces having zero boundary fluxes from the difference operator.

#### (ii) Extended finite-element method discretization for plane strain fractures

The fundamental idea of the XFEM is to represent cracks by augmenting the standard Lagrange shape functions in the elements intersecting and neighbouring the cracks by two types of specialized enrichment functions, namely: sign enrichment to represent the crack geometry
3.5in which is the signed distance function and is the extension of the crack to include any element that is cut by the crack; and tip enrichment to represent the singular behaviour of the form [*u*]=*As*^{α} at the tip in order to restore the *O*(*h*^{2}) convergence expected of bi-linear elements. The following basis functions for tip enrichment corresponding to an arbitrary power law 1/2≤*α*<1 have been derived [46] by determining the solution to the equilibrium equations (2.2) subject to power-law displacement jump conditions:
3.6and
3.7

The finite-dimensional Galerkin space *U* is defined by *U*=*H*^{1,h}×*H*^{1,h}×*H*^{1,h} and is spanned by the following shape functions:
3.8where , *N*_{i} are the standard piecewise bi-linear Lagrange basis functions on rectangular elements, and *a*_{i}, *b*_{i}, . Let *I* be the set of all nodes, *I*_{t} the set of tip enriched nodes that are within a prescribed distance *ρ* from the crack tips, *I*_{s} the set of sign enrichment nodes comprising all the nodes of elements cut by the crack excluding those already in *I*_{t}, so that , *I**_{s} is the set of all nodes of elements that are cut by the crack *and* that have at least one node in *I*_{s}, and *I**_{t} is the set of all nodes in elements that have at least one node in *I*_{t}. In addition, if tip enrichment is being used, the two blending functions and are introduced to blend the two enrichments.

The one-dimensional lubrication equation corresponding to (2.8) is discretized using a finite-volume formulation in which fluxes at element edges are obtained by differencing pressures evaluated at element centres [48].

### (b) Solving the coupled pressure–width equations

#### (i) Iteration schemes for the coupled equations

As the iterative schemes for solving the coupled equations for the DD and XFEM formulations are similar, we choose to discuss the solution of the coupled DD equations (3.1) and (3.3) for the volume conserving-equilibrating fluid pressure *p*_{f} and width *w* vectors at time *t*. In this section, we will assume that the fracture footprint is known (or the pressure and width fields for a trial footprint are being sought) and defer the discussion of determining the evolution of to the next subsection. Solving the coupled equations (3.1) and (3.3) for *p*_{f} and *w* as primary variables requires the inversion of the operator *A*(*w*), which is rank deficient owing to the Neumann boundary condition on *p*_{f} due to the zero flux boundary condition (2.15*c*) that holds in the zero lag case. In order to ensure a solution exists in this case it is necessary to impose the solvability condition
3.9where is the fluid volume leaked from the fracture over the current time step and (3.9) is obtained by integrating (2.7) in space over the region and in time over [*t*,*t*+Δ*t*], using the divergence theorem, and imposing the zero flux boundary condition (2.15*c*).

To avoid the complexity of having to impose the solvability condition (3.9), we use (3.1) to eliminate *p*_{f} from (3.3) to obtain the following evolution equation for the width *w* at time *t*:
3.10Once the width has been determined, the corresponding net and fluid pressures can be obtained directly from (3.1). We now outline a number of iterative schemes to solve these nonlinear equations.

*Fixed point iteration.* We represent the total change in width over the time step by Δ*w*=*w*−*w*_{t−Δt} and use the relation *w*=*w*_{t−Δt}+Δ*w* to rewrite (3.10) in the form
3.11in which the fixed point operator is defined as . Starting with an initial guess for *w* the updated width change Δ*w* is determined by inverting and a new estimate of *w* is obtained from the relation *w*=*w*_{t−Δt}+Δ*w*, which is in turn used to update *A*(*w*) on the right side of (3.11) and in . The process is repeated until the relative changes in *w* become less than a certain tolerance.

*Picard iteration.* The Picard iteration scheme also uses the fixed point operator to obtain a new estimate *w*_{j+1/2} for *w*, which is then relaxed using the parameter *γ*
3.12and
3.13Comparing (3.12) with (3.11), we observe that they are related in that each involves the inversion of the same operator with each iteration but differ in their right-hand sides.

If we define *δw*_{j}=*w*_{j+1}−*w*_{j} to be the change in *w* between two iterates (as opposed to Δ*w* the change in *w* over the time step [*t*,*t*+Δ*t*] used in the fixed point scheme) and combine (3.12) and (3.13), we observe that the Picard scheme can be rewritten in the form
3.14where *r*_{j}=*w*_{j}−*w*_{t−Δt}−Δ*tA*(*w*)(*Gw*_{j}+*σ*_{0})−*S* is the residual vector for (3.10) associated with *w*_{j}.

*Newton iteration.* The Newton scheme can be derived by substituting a perturbation *w*=*w*_{j}+*δw*_{j} to the current estimate *w*_{j} into (3.10), expanding the nonlinear terms using a Taylor series, and retaining only first order terms to yield the linear equation
3.15where is the Jacobian and *B*_{j}=*A*′(*w*)*p*_{f}.

#### (ii) Fast inversion of the operators and

The complexity involved with inverting the non-symmetric fixed point and Jacobian operators is that they both involve the matrix product *A*(*w*)*G* of the sparse matrix *A*(*w*) and the fully populated elasticity matrix *G*. For problems with a large number of degrees of freedom *N*, even storing the elasticity matrix *G* can be prohibitive, so iterative solvers for such non-symmetric systems such as GMRES or BiCGSTAB [77] are often preferred to direct solvers. The efficacy of these iterative techniques is dependent on the cost of performing a *G*-matrix vector product and the spectra of the operators and , which determines the number of iterations required. If the elastic medium is homogeneous the convolution structure of the matrix (3.2) implies that it can be represented by a bi-circulant Toeplitz matrix, which has a sparse representation in terms of the fast Fourier transform. This representation makes it possible to perform a matrix vector multiply by *G* in operations [78]. Alternative methods for performing a *G*-matrix vector product efficiently, which can be extended to unstructured meshes, are the fast multipole method (FMM) [79–81] or the use of H-matrices [82,83]. Both these methods involve approximations whose errors need to be traded against speed, and the feasibility of developing a FMM for a layered elastic material or the utility of the H-matrix approach if *G* needs to be updated frequently are questions that still need to be resolved. Finally, the number of iterations can be reduced significantly by implementing incomplete lower upper factorizations to approximations to and that are based on sparse approximations to *G* that only involve nearest-neighbour interactions [84].

The numerical linear algebra involved with XFEM formulations presents a wider variety of options than the DD formulation, which only involves unknowns on the fracture plane. One option is to couple the sparse symmetric stiffness matrix for the solid medium, involving many more degrees of freedom than the DD formulation, to the fluid flow process taking place in the fracture plane. While the sparse symmetric matrices might be desirable, is it efficient to couple the relatively few degrees of freedom needed to capture the fluid flow to the much larger number of degrees of freedom to capture the solid deformation? Another approach is to use the XFEM to determine a numerical Green’s function, which involves a much smaller dense matrix. The relative benefits of these two approaches is still an open question in this fully coupled nonlinear environment.

### (c) Locating the free boundary using the implicit level set algorithm

As there are an infinite number of volume preserving-equilibrating pairs (*w*,*p*_{f}) that are solutions to the nonlinear equations (3.1) and (3.3), each associated with a slightly different footprint , it is important to select the correct solution by determining the appropriate fracture boundary . Owing to the singular free boundary problem characteristic of HF propagating with zero lag, the local fluid velocity in the vicinity of the tip, which requires the numerical evaluation of a distinguished limit, cannot be used to locate the free boundary .

Indeed, the algorithm to locate the free boundary has to incorporate the singular tip asymptotics discussed in §2b(iv). We now describe the recently developed implicit level set algorithm (ILSA) [32,85] that uses the tip asymptotics to locate the free boundary. It is notionally convenient to separate the DD elements representing the fracture into two sets, namely: the tip elements and the channel elements . The tip elements are the partially filled elements coloured by light shading in figure 1*a*, while the channel elements comprise those elements surrounded by the tip elements that are completely filled with fluid. We also identify the ribbon of channel elements on the boundary of the channel region that share at least one side with a tip element (see the dark shaded elements in figure 1*a*), which we refer to as survey elements.

Steps in the ILSA procedure:

(i) *Determine equilibrium pressures and widths for the current footprint *

Given a trial footprint (initially taken to be that at the previous time step ) the coupled equations (3.10) are used to determine the width *w* (and hence *p*_{f}) that corresponds to the fluid flux into the well-bore over the current time step. In the first front iteration, the crack widths are treated as the primary variable over the whole fracture . In the second and subsequent front iterations, as estimates of the asymptotic tip widths are known (see step (iv) below), the fluid pressures are treated as the primary unknown in the tip region and the widths are treated as the primary unknown in the channel region .

(ii) *Use tip asymptote to find shortest distances from elements in to the new trial crack front *

Given the new width estimate *w* at each of the survey elements in the current trial fracture widths are used to invert the appropriate asymptotic relation from §2b(iv) to determine the shortest distance *s* from the centres of these elements to the free boundary . In particular, (2.19*c*) can be inverted to obtain *s* given *w*_{k}, (2.21) can be inverted to obtain *s* given *w*_{m}, (2.22) can be inverted to obtain *s* given , the generalized asymptote (2.24) can be inverted to find *s* given *w*_{mk}, and given the generalized asymptote (2.25) can be used to find *s*.

We observe that all the tip asymptotes (except for the toughness asymptote (2.19*c*)) involve the unknown front velocity *V* , which is so crucial to locating the free boundary. To resolve this issue, if *s* is the shortest distance from the current collocation point to the trial free boundary, then in each asymptote involving the front velocity explicitly we replace *V* by its difference quotient approximate
3.16The generalized *mk*-asymptote (2.24) has recently been implemented for plane strain fractures [48] and for planar fractures [85]. For each survey point in , the trial width *w* is used to find *s* from
The generalized -asymptote (2.25) has been implemented in planar [74] and multi-planar HF models [75]. For each survey point in the trial width *w* is used to determine *s* from

(iii) *Solve the eikonal equation to locate the new trial crack front *

The set of shortest distances from elements in to the new trial front establishes the following initial conditions:
3.17to construct the signed distance function by solving the eikonal equation
3.18The negative sign in the initial condition (3.17) enforces the sign convention that for all points (*x*,*y*) that lie within the fracture boundary curve , while points for which lie outside . Moreover, because the initial condition (3.17) defines the distance to the free boundary, the fracture boundary curve is defined by the level set .

(iv) *Set the tip volumes from the asymptote to impose the tip asymptotics in a weak form*

Having determined a new estimate for the location of the free boundary the volume of fluid in each tip element can be determined by integrating the applicable tip asymptote backward from the front to the opposite vertex of a triangular region shown in figure 1*b*
3.19where ℓ is the distance from the front to the farthest interior corner of a tip element and *θ* is the angle that the local outward normal to the front makes with the tip element edge. The volume over a triangular region can be used to construct the fracture volume within a partially filled rectangular element as follows:
3.20In order to evaluate these integrals only the first two moments of the tip asymptote are required.

The average width *w*^{t} in the tip element is set to
As the width in the tip element is determined from the asymptotics once the front position is set, the fluid pressure in the tip element now becomes the primary unknown. The equation for the fluid pressure comes from integrating the continuity equation (2.7) over the tip element and applying the zero flux boundary condition. This is precisely the procedure used in the finite-volume method and results in an expression for the change in fluid volume over the time step involving that should match the flux of fluid across the boundaries, which involves the fluid pressure at the centre of the tip element.

(v) *Convergence check*

If the iterations on the fracture front have converged then proceed to the nest time step, if not return to (i) with the new estimate of .

We observe that the concept of a partially filled element with tip asymptotics imposed in a weak sense is remarkably effective at representing a fracture with an arbitrary front on a structured mesh. The ability to capture the tip asymptotics in a weak form makes it possible to account for multiscale effects that act at length scales that are many orders of magnitude smaller than the computational mesh size. Modelling such multiscale behaviour directly would require prohibitively fine meshes. Given that the material parameters are assumed to be fixed, the tip velocity *V* is the primary determinant of the region in phase space in which a particular asymptote finds itself. As the ILSA scheme is able to autonomously determine the front velocity *V* from the asymptote, it is able to automatically adjust the point in phase space to determine the propagation regime of the fracture front in the neighbourhood of the survey points. In a recent study comparing the performance of a variety algorithms [76], the ILSA scheme was shown to provide accurate solutions while only requiring a relatively coarse mesh.

The ILSA scheme has also been implemented in XFEM HF models. In order to be able to impose different boundary conditions in the tip and channel regions and , respectively, it is necessary to use the mixed variational formulation (2.5). The so-called XFEM-*t* scheme [48] uses sign enrichment in the channel region and tip enrichment in the tip regions and exhibits the optimal *O*(*h*^{2}) convergence rate one would expect of bi-linear Lagrange elements in two dimensions. However, the coefficients in the stiffness matrix corresponding to tip elements have to be updated for each new trial fracture front position. As the singular integrations associated with this constant updating process can be computationally intensive, an alternative, sub-optimal so-called XFEM-*s* scheme [48] was developed. The XFEM-*s* scheme uses only sign enrichment on a region that both encompasses the crack and extends beyond the crack by linear extrapolation to the edge of the element in which the crack tip falls. Similar to the DD ILSA scheme the tip asymptote is imposed in a weak sense within this extended sign enriched region. Although the XFEM-*s* scheme yields sub-optimal *O*(*h*) convergence, it is considerably more efficient.

## 4. Numerical results

In this section, we present numerical solutions that illustrate the performance of the ILSA strategy for four different examples. In the first example, we consider a fracture that propagates in an elastic medium in which there are positive jumps in the confining stress field located symmetrically with respect to the fluid source and compare the ILSA solution with the results from a laboratory experiment. In the second example, we consider a fracture that propagates in an impermeable elastic medium in which *K*′>0 and the confining stress field drops discontinuously across two interfaces located symmetrically with respect to the fluid source. This example is chosen to illustrate the way in which the new ILSA scheme is able to use the *mk*-generalized asymtptote to capture the multiple length scales that are active at different points along the perimeter of the fracture, depending upon the magnitude of the local normal velocity. The third example provides an ILSA model of a buoyancy-driven *m*-vertex HF that evolves into a finger-like fracture when emanating from a point source and a planar front that breaks up into well-defined fingers when emanating from a line source. The fourth example provides an application of the XFEM to analyse the propagation of a two-dimensional plane strain HF in the vicinity of a free surface. The HF is initiated parallel to the free surface but as its length becomes comparable with the distance to the free surface the fracture starts to curve towards the free surface. The results are compared with those of OribiC—a DD-based HF code.

### (a) A symmetric stress barrier: comparison with a laboratory experiment

In this section, we compare the ILSA *m*-vertex solution with the results of a laboratory experiment [86] in which an HF was initiated at the centre of a 340×355×360 mm poly(methyl methacrylate) (PMMA) block whose surface was specially machined to match the displacement due to strip loads in order to establish a symmetric stress jump within the block. Within the block a layer of height *H*=50 mm is bounded by symmetric stress jumps Δ*σ*=4.3 MPa in *σ*_{0}, which is given by
The ILSA model assumes that the fracture grows within an infinite elastic medium subject to the above confining stress field *σ*_{0} and uses the following parameter values reported for the experiment:
As the fracture takes place between two unbonded impermeable blocks of PMMA, the fracture in the experiment propagates in the *m*-vertex viscous regime. In figure 2, the experimental (blue) and ILSA (black) blade-like fracture footprints and fracture width cross sections with the planes *y*=0 and *x*=0 are compared at a sequence of sample times. The numerical results show remarkable agreement with the experimental results.

### (b) A stress drop: distinct propagation regimes along the periphery

In this section, we illustrate the ILSA solution for an HF propagating in an infinite elastic medium subject to a confining stress field *σ*_{0} that drops discontinuously across two interfaces a distance *H*=20 m apart. The fracture initiates midway between these two interfaces and the toughness and viscous processes compete to determine the location of the advancing fracture front and exhibits a continuous variation in the propagation regime all along the fracture perimeter. The material parameters used in this example are
while the confining stress field is defined by

These results were generated using the ILSA scheme [85] in which the two-process generalized *mk*-edge asymptote (2.24) has been used. On the left part of figure 3, the generalized *mk*-edge asymptote [70] is represented by the solid line that changes colour from red to blue as the asymptote transitions from the *k* to *m* vertex solutions, which are denoted by red to blue dotted lines, respectively. On the right part of figure 3, the ILSA-*mk* fracture footprints at *t*=100, 400, 900, 1800, 3600 s are plotted while the interfaces across which the confining stress field jumps are indicated by the thick horizontal black lines, and the underlying computational mesh is indicated by the thin black squares. For each of these footprints the survey elements are shaded with a colour that corresponds to the length scale active within that element. Thus, we observe that for the earliest time step *t*=100 s most points along the fracture perimeter are advancing close to the viscosity dominated regime, while by time *t*=400 s all points along the perimeter are advancing in the intermediate regime in which . At subsequent time samples the fracture is propagating predominantly in the toughness regime in the top layer *y*>10 m in which the confining stress is the highest so that the fracture propagates much slower than in the other layers. Because the least resistance to the progress of the fracture is encountered in the bottom layer *y*<−10 *m* with the lowest confining stress, the fracture herniates into this region advancing much more rapidly and propagating in the intermediate regime close to the viscous asymptote. In the middle layer |*y*|<10 m with the intermediate confining stress, the fracture front develops an inflection point and the regime of propagation transitions from the intermediate asymptotic region at the bottom interface to the toughness regime at the top interface.

### (c) Implicit level set algorithm model of fingering associated with a buoyancy-driven crack

In this example, we consider the evolution of an HF under a confinement field that decreases linearly with increasing *y* according to
This is a model for a buoyancy-driven HF [2,4,5] in which the fracture finds it easier to grow into the region in which the confinement is steadily decreasing. The linear dependence of *σ*_{0} on *y* introduces an additional convective term
on the right side of (2.8), which can be seen by eliminating the fluid pressure *p*_{f} in favour of the net pressure. In this case, the HF propagates using the *m*-vertex asymptote within a plane perpendicular to the confining stress field *σ*_{0} in an infinite elastic medium and for which the following dimensionless parameters were used:

If fluid is injected from a point source the ILSA scheme using the *m*-vertex viscous asymptote yields a single finger-like fracture for which the fracture fronts are shown in figure 4*a*. If fluid is injected from a line source
then from figure 4*b* it can be seen that the ILSA *m*-vertex buoyancy-driven HF develops an unstable advancing front with well-developed fingers [87].

### (d) Extended finite-element method–displacement discontinuity comparison for a curving crack near a free surface

In this section, we illustrate an XFEM solution for an HF propagating in a state of plane strain close to a free surface [45]. The following dimensionless parameters were used: so that the dimensionless toughness , where

A fracture of half-length 18 dimensionless units was initiated parallel to the free surface a distance *H*=45.5 dimensionless units from it. A rectangular computational domain *L*_{x}=300 and *L*_{y}=121 of square finite elements of mesh size Δ*x*=Δ*y*=1 was combined with singly and doubly infinite elements of order 9 along the bottom and lateral two sides of the 300×121 rectangle to complete the half-space. A comparison between the XFEM results and those of OribiC, a DD HF code that incorporates Green’s function influences for a half-space, is shown in figure 5. We observe that as the length of the HF starts to approach *H* it becomes easier for the fracture to grow in the direction of the free surface. The OribiC and XFEM curved crack paths (figure 5*a*) show close agreement. In figure 5*b*, the XFEM fracture path is superimposed on the underlying FEM mesh. The OribiC and XFEM widths (figure 5*c*), plotted against the arc length *s* from the injection point along the fracture, show close agrement at all the sample times *τ*_{*}=[0.2,0.4,0.6]. The OribiC and XFEM pressures (figure 5*d*) are plotted against the arc length *s* from the injection point along the fracture until the boundary of the fluid front. These two solutions show close agrement. We observe that because the fracture is close to the free surface and the confinement *σ*_{0}=0 there is a non-zero fluid lag in this problem. The transition between the fluid-filled fracture and the stress free part of the fracture is responsible for the inflection point in the width field shown in figure 5*c*.

## 5. Conclusion

In this review paper, we have described continuum models of HF propagation and highlighted some of the unique challenges associated with these models, particularly when the fluid and fracture fronts coalesce. We have demonstrated that in this case the evolution of the fracture is a singular free-boundary problem in which the normal velocity of the free boundary is not easily evaluated numerically, but has to be determined by delicate asymptotic analysis. We have outlined this asymptotic analysis that highlights the way in which the many physical processes compete to yield an asymptotic solution that varies on multiple length scales. Precisely how to capture all these multiple length scales without having to resolve both the very finest and coarsest length scales was an open problem until the ILSA was developed in 2008 [32]. In this paper, we describe the class of numerical schemes that are derived from the ILSA paper. The ILSA strategy locates the free boundary using the tip asymptotics; once the free boundary has been found, the multiscale tip asymptotic behaviour is imposed in a weak sense. This two-step iterative process enables the free boundary to be located implicitly in a way that is consistent with the tip asymptotics. Moreover, the multiscale behaviour can be accounted for right down to the finest length scale while still computing on a relatively coarse-structured mesh. This review paper illustrates the ILSA strategy using both the DD formulation and the XFEM formulation for the fracture modelling. We also address some recently resolved practical issues, including new models for proppant transport able to model ‘tip screen-out’; efficient numerical schemes to solve the coupled nonlinear equations and fast methods to solve the resulting linear systems. Finally, numerical examples are provided to illustrate the performance of the method, including comparison with a symmetric stress jump laboratory experiment; a stress drop problem that exhibits a two-process-driven change in length scale all around the perimeter of the fracture; buoyancy-driven fingering; and an XFEM model of an HF curving towards a free surface.

There are a number of open problems still to be addressed, including:

(i)

*Accurate numerical models for fracture recession after shut-in*. To date, there are no closed-form solutions and tip asymptotics for recession of the fracture front due to shut-in, flow-back and leak-off. This becomes particularly delicate when the fracture closes on proppant. Closed-form solutions along with precise laboratory experiments would be very useful in developing stable numerical codes that are able to capture recession accurately.*Q*_{0}=0 and flow-back*Q*_{0}<0(ii)

*Plasticity in the solid medium and poroelastic effects for capturing the effects of long periods of injection*. With the advent of domain-based methods such as the XFEM, it will become easier to model not only heterogeneities in the material properties but also nonlinear effects such as plasticity in the vicinity of the fracture and the poroelastic effects of long injection times.(iii)

*Evolution, interaction and intersection of arbitrary fracture surfaces in three dimensions*. Capturing the propagation, interaction and possible intersection of multiple arbitrary fracture surfaces presents the numerical modeller with considerable challenges. From a continuum perspective, the most elegant representation of all these effects is the phase field approach [13–15] in which a single damage field variable is used to track and capture the evolving fractures in the solid three-dimensional medium. The primary advantage of both phase field and discrete models is that they are able to capture the evolution and interaction of complex fracture geometries with significantly less effort than the conventional continuum models in which cracks are modelled explicitly. On the other hand, phase-field models are extremely computationally intensive as they involve the solution of a non-convex optimization problem at each time step. Moreover, due to the discrete nature of lattice models and the smeared-out damage representation of cracks by the phase-field approach, it is not clear whether these methodologies will be able, without using prohibitive computing resources, to capture the complex multiscale behaviour characteristic of HF when multiple physical processes compete to determine their evolution. Fully coupled versions of these algorithms able to capture the full multiscale tip asymptotics would hold much promise. However, issues such as computation time and explicit versus implicit time stepping still need to be resolved.(iv)

*Implicit time stepping and advanced numerical linear algebra*. In the solution of the coupled equations there remain many open questions. For example, in the XFEM, is it efficient to couple the relatively few degrees of freedom needed to capture the fluid flow with the much larger number of degrees of freedom needed to capture the solid deformation? Are there special numerical linear algebra techniques that need to be developed to make this process optimal, as opposed to searching for a numerical Green’s function?(v)

*Using numerical HF models for monitoring and control*. With numerical models becoming more efficient is it possible to make better use of field measurements to improve the monitoring of propagating HF using tiltmeters or microseismic data, for example? There has been some degree of success using tiltmeter data along with the Kalman filter to monitor propagating HF using synthetic data [88,89] and field data with monitoring wells from mines [90].

## Competing interests

I declare I have no competing interests.

## Funding

The author gratefully acknowledges the grant support of NSERC (Canada).

## Acknowledgements

I thank Prof. E. Detournay (University of Minnesota) for convincing me of the crucial importance of the tip asymptotics in HF modelling, Prof. A. P. Bunger (University of Pittsburgh) for providing experimental results, and Prof. E. Dontsov (University of Houston) and Dr E. Gordeliy (SDR-Boston) for insightful feedback on this paper.

## Footnotes

One contribution of 12 to a theme issue ‘Energy and the subsurface’.

- Accepted July 22, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.