**2. Benchmarks**

#### **2.1. Plane-strain hydraulic fracture (KGD)**

The case of a plane strain hydraulic fracture driven by a Newtonian fluid under constant injection rate is also sometimes refered to as the KGD model (for Khristianovic [12], Geerstma and De klerk [13]). In the absence of leak-off the solution of the hydraulic fracture propagation is self-similar and depends on a single dimensionless number: i.e. a dimensionless toughness K

$$\mathcal{K} = \frac{K'}{E'^{3/4} \mu'^{1/4} \mathcal{Q}\_o^{1/4}}$$

where *Qo* is the volumetric injection rate per unit length in the out-of plane direction, *<sup>E</sup>*′ denotes the plane-strain elastic modulus, *<sup>µ</sup>*′ <sup>=</sup> <sup>12</sup>*µ<sup>f</sup>* is an equivalent viscosity (with *<sup>µ</sup><sup>f</sup>* the fluid viscosity), and *<sup>K</sup>*′ <sup>=</sup> <sup>√</sup>32/*<sup>π</sup> KIc* with *KIc* the fracture toughness (see [14] for more details). Equivalently a dimensionless viscosity M = K−<sup>4</sup> can be used. The complete solutions for the fracture length evolution, fracture width and net pressure have been obtained for the limiting cases of zero dimensionless toughness (equivalently infinite M) and zero dimensionless viscosity (infinite K) First order solutions for either small toughness or viscosity are also available [3, 4]. Semi-analytical solutions for any finite values of dimensionless toughness or viscosity are also available [5].

2 Effective and Sustainable Hydraulic Fracturing

of injection (zero leak-off).

propagation.

the light of different benchmarks.

dimensionless toughness K

**2.1. Plane-strain hydraulic fracture (KGD)**

**2. Benchmarks**

Viscosity M

Toughness K

**Propagation regimes Plane-strain Radial**

**Table 1.** List of available solutions for the propagation of hydraulic fractures driven by a Newtonian fluid under constant rate

In geoscience applications, hydraulic fractures propagate in a complex, often poorly characterized medium. Nevertheless, the description of the medium must be simplified in order to apply theoretical models. It is thus crucial that numerical implementations of such models for fracture growth be accurate such that differences from field observations can be

In the last ten years, a number of reference solutions (analytical and semi-analytical) have been obtained for propagating plane-strain [1–5] and radial hydraulic fractures [6, 7] (see Table 1). These solutions provide invaluable benchmarks for numerical simulators. We compare a number of simulators (2D and 3D) that use different propagation algorithms against these reference solutions for hydraulic fractures driven by a Newtonian fluid under a constant injection rate. For the sake of clarity, we do not address fluid leak-off in our discussion. Of particular interest is the accuracy of the different simulators in tracking the moving fracture front, particularly in the so-called viscosity-dominated regime of

An outstanding question relates to the convergence and robustness of numerical simulators with respect to the multiscale near-tip behavior of hydraulic fractures. The coupled lubrication (fluid flow) and elasticity equations are known to degenerate near the fracture tip, such that the solution of a semi-infinite fracture propagating at a constant velocity is characterized by a multiscale singular behavior near the tip [8–10]. The nature of the dominant singularity depends on the relative importance of two dissipative processes (viscous forces and fracture energy), as well as the reference length scale. Such a multiscale behavior near the fracture tip in turn governs the evolution of the velocity of a finite hydraulic fracture during injection. We discuss the degree to which a numerical simulator needs to include and resolve the near-tip behavior in order to accurately match reference solutions in

The case of a plane strain hydraulic fracture driven by a Newtonian fluid under constant injection rate is also sometimes refered to as the KGD model (for Khristianovic [12], Geerstma and De klerk [13]). In the absence of leak-off the solution of the hydraulic fracture propagation is self-similar and depends on a single dimensionless number: i.e. a

<sup>K</sup> <sup>=</sup> *<sup>K</sup>*′

*E*′3/4*µ*′1/4*Q*1/4

*o*

attributed to model assumptions rather than poor numerical solution.

(<sup>K</sup> <sup>=</sup> 0) [3] (with correction for small toughness) [6]

(K → <sup>∞</sup>) [4] (with correction for small viscosity) [6, 11]

We will restrict our comparisons to the case of relatively small dimensionless toughness (e.g. K < 1), which is known to be the more difficult condition to reproduce numerically. Therefore we express the solution in a so-called viscosity scaling. Because the solution is self-similar the time dependence can be obtained using dimensional analysis. We aim to compare the solutions provided by different numerical codes, which are typically developed in space-time. We thus introduce a dimensionless time *τ* = *t*/*tc* and a scaled coordinate *ξ* = *x*/ℓ, where *tc* is a characteristic time scale and ℓ is the fracture length. The fracture length, opening, and net pressure can be written as follows:

$$\begin{split} \ell &= \left(\frac{E'\Omega\_o^3 t^4}{\mu'}\right)^{1/6} \gamma\_m(\mathcal{K}) = \left(\frac{E'Q\_o^3 t\_c^4}{\mu'}\right)^{1/6} \underbrace{\tau^{2/3} \gamma\_m(\mathcal{K})}\_{\gamma(\tau,\mathcal{K})} \\ w &= \frac{Q\_o^{1/2} t^{1/3} \mu'^{1/6}}{E^{1/6}} \Omega\_m(\underline{\xi}, \mathcal{K}) = \frac{Q\_o^{1/2} t\_c^{1/3} \mu'^{1/6}}{E^{1/6}} \underbrace{\tau^{1/3} \Omega\_m(\underline{\xi}, \mathcal{K})}\_{\Omega(\underline{\xi}, \tau, \mathcal{K})} \\ p &= \frac{E'^{2/3} \mu'^{2/3}}{t\_c^{1/3}} \Pi\_m(\underline{\xi}, \mathcal{K}) = \frac{E'^{2/3} \mu'^{2/3}}{t^{1/3}} \underbrace{\tau^{-1/3} \Pi\_m(\underline{\xi}, \mathcal{K})}\_{\Pi(\underline{\xi}, \tau, \mathcal{K})} \end{split} \tag{1}$$

where we have highlighted the correspondence between results obtained using a time-based algorithm (say *γ*, Ω, Π) to the self-similar solution dimensionless solution F*m*(K) = {*γm*, Ω*m*, Π*m*}.

The dimensionless solution F*m*(K) = {*γm*, Ω*m*, Π*m*} for small toughness developed in [3] will be compared with the numerical solutions from different simulators. More precisely, for three small values of dimensionless toughness K = 0.01, 0.1, 0.5, we will focus on the comparisons of the dimensionless fracture length *γm*, opening profile Ω*m*(*ξ*) close to the fracture tip and error in the fracture volume etc. We are especially interested in the evolution of the error of the numerical solutions with respect to mesh sizes in the near-tip region of the fracture, where gradients are the largest. The solution, for small toughness (K < 1), is in fact governed by the hydraulic fracture viscosity tip asymptote: the opening behaves as *<sup>w</sup>* ∼ (ℓ − *<sup>x</sup>*)2/3 and the net pressure as *<sup>p</sup>* ∼ (ℓ − *<sup>x</sup>*)−1/3 close to the fracture tip (see [10] for details). The tip region affected by the asymptote actually extends to about 10 to 20 percent of the plane-strain fracture length for dimensionless toughnesses below 0.5.

10.5772/56212

(3)

859

*<sup>R</sup>* <sup>=</sup> *<sup>E</sup>*′1/9*Q*1/3

*<sup>p</sup>* <sup>=</sup> *<sup>E</sup>*′2/3*µ*′1/3

simulators to this zero-toughness/small time solution.

focusing mostly on fracture length versus time.

**2.3. Some practical numbers**

*<sup>w</sup>* <sup>=</sup> *<sup>Q</sup>*1/3

*<sup>o</sup> t* 4/9

1/9 *<sup>E</sup>*′2/9 <sup>Ω</sup>*m*0(*ρ*) =

*<sup>o</sup> <sup>µ</sup>*′2/9*<sup>t</sup>*

*<sup>µ</sup>*′1/9 *<sup>γ</sup>m*<sup>0</sup> <sup>=</sup> *<sup>E</sup>*′3*Qoµ*′

*<sup>t</sup>*1/3 <sup>Π</sup>*m*0(*ρ*) = *<sup>K</sup>*′<sup>3</sup>

where again, we have highlighted the correspondence between the zero-toughness/ M-vertex self-similar solution F*m*<sup>0</sup> = {*γm*0, Ω*m*0, Π*m*0}, which is independent of time and the dimensionless solution expressed as a function of the previously defined dimensionless time. The zero-toughness/M-vertex solution F*m*<sup>0</sup> has actually been found to correctly capture the propagation of hydraulic fractures up to a dimensionless toughness of K = 1, i.e. for dimensionless time *<sup>τ</sup>* ≤ <sup>1</sup> (*<sup>t</sup>* ≤ *tmk*) [6]. We will thus investigate the convergence of different

For dimensionless time above unity (*τ* > 1 (*t* > *tmk*)), the solution transitions from the viscosity dominated (early-time) to the toughness (large-time) dominated regime. The toughness dominated regime is reached for K ≈ 3.5 (*τ* ≈ 70000). Note that, for infinitely large dimensionless toughness (i.e. zero viscosity), the solution is also self-similar [6] and is also denoted as the K-vertex solution. We will also briefly investigate, for a subset of the simulators considered, the transition between these two regimes of propagation (M to K),

Although rock properties and stimulation practices vary, it is interesting to compute the scales and dimensionless numbers previously introduced. Let's assume a "tight" rock with the following realistic properties: a plane-strain Young's modulus of 40 GPa and a fracture toughness of 1.5 MPa.m1/2. First, for a hydraulic fracturing treatment using a highly viscous fluid (e.g. gel-like with *µ<sup>f</sup>* = 100 cPoise) at a practical rate of 10 Barrels per minute, we obtain a transition time scale *tmk* for a radial fracture of 4.2 10<sup>6</sup> seconds! The propagation of such a hydraulic fracture for a realistic injection duration (i.e. less than two hours) will always be in the viscosity dominated regime of propagation. Remember that the dimensionless toughness evolves as (*t*/*tmk*)1/9 for radial fractures. For the plane-strain geometry (where the injection rate is per meter in the out-of-plane direction), using the same parameters we obtain a dimensionless toughness of 0.26 clearly indicating a viscosity dominated propagation.

For a slick-water treatment, popular in shale-gas reservoirs, the injection rates are usually much higher in order to compensate for the low viscosity of water and to obtain a sufficiently wide fracture to accomodate proppant (see the scales in front of the opening *w* in Equations (1) and (3)). Using a value of 20 Barrels per minute (a realistic value for a single perforation

*<sup>K</sup>*′<sup>4</sup> *<sup>τ</sup>*4/9*γm*<sup>0</sup> *γ*(*τ*)

> *<sup>K</sup>*′ *<sup>τ</sup>*1/9Ω*m*0(*ρ*) Ω(*ρ*,*τ*)

> > Π(*ρ*,*τ*)

Accuracy and Convergence of Hydraulic Fracture Simulators

http://dx.doi.org/10.5772/56212

*<sup>E</sup>*′3/2*Qoµ*′ *<sup>τ</sup>*−1/3Π*m*0(*ρ*)

*E*′*Qoµ*′

#### **2.2. Radial hydraulic fracture**

The growth of a radial hydraulic fracture spans both the viscosity and toughness regimes of propagation [6, 11]. At early times, the perimeter and the opening of the fracture are small and most of the energy is spent in viscous flow, whereas at a later times, the fracture perimeter and opening are larger and the fracture energy required to extend the fracture dominates the energy required to drive the viscous fluid through the fracture. The radial solution is also dependent only on a dimensionless toughness which in this case is a function

of time. Introducing the characteristic time *tmk* = *<sup>µ</sup>*′5*Q*<sup>3</sup> *oE*′<sup>13</sup> *<sup>K</sup>*′<sup>18</sup> 1/2 , and the dimensionless time *τ* = *t*/*tmk* we have (see [6, 14] for more details):

$$\mathcal{K} = \tau^{1/9}$$

Solutions for the case of zero and infinite dimensionless toughness (i.e. small and large dimensionless time) have been obtained semi-analytically [6]. The complete transient solution can be obtained only numerically. A reference algorithm [7] based on an explicit moving mesh algorithm with proper matching of the multiscale HF tip asymptotics [10] will provide the baseline for the comparisons for intermediate times. The fracture radius *R*, width *w* and net pressure *p* can be written as:

$$\begin{aligned} R &= \frac{E'^3 Q\_o \mu'}{K'^4} \gamma(\tau) \\ w &= \frac{\sqrt{E' Q\_o \mu'}}{K'} \Omega(\rho, \tau) \\ p &= \frac{K'^3}{E'^{3/2} Q\_o^{1/2} \mu'^{1/2}} \Pi(\rho, \tau) \end{aligned} \tag{2}$$

where the dimensionless solution F = {*γ*, Ω, Π} depends only on dimensionless time *τ* = *t*/*tmk* and scaled position *ρ* = *r*/*R* along the fracture.

As before, we will pay particular attention to the case of small dimensionless toughness, i.e. early-time, which is the most challenging numerically. In the limit of zero-toughness/early-time (i.e. viscosity dominated propagation, here refereed as the M-vertex), the solution is self-similar and can be conveniently written in the following viscosity scaling:

$$\begin{aligned} R &= \frac{E^{1/9} Q\_o^{1/3} t^{4/9}}{\mu^{1/9}} \gamma\_{m0} = \frac{E^{3} Q\_o \mu'}{K'^4} \underbrace{\tau^{4/9} \gamma\_{m0}}\_{\gamma(\tau)} \\ w &= \frac{Q\_o^{1/3} \mu'^{2/9} t^{1/9}}{E^{2/9}} \Omega\_{m0}(\rho) = \frac{\sqrt{E^{1} Q\_o \mu'}}{K'} \underbrace{\tau^{1/9} \Omega\_{m0}(\rho)}\_{\Omega(\rho, \tau)} \\ p &= \frac{E^{2/3} \mu'^{1/3}}{t^{1/3}} \Pi\_{m0}(\rho) = \frac{K^3}{E'^{3/2} \sqrt{Q\_o \mu'}} \underbrace{\tau^{-1/3} \Pi\_{m0}(\rho)}\_{\Pi(\rho, \tau)} \end{aligned} \tag{3}$$

where again, we have highlighted the correspondence between the zero-toughness/ M-vertex self-similar solution F*m*<sup>0</sup> = {*γm*0, Ω*m*0, Π*m*0}, which is independent of time and the dimensionless solution expressed as a function of the previously defined dimensionless time.

The zero-toughness/M-vertex solution F*m*<sup>0</sup> has actually been found to correctly capture the propagation of hydraulic fractures up to a dimensionless toughness of K = 1, i.e. for dimensionless time *<sup>τ</sup>* ≤ <sup>1</sup> (*<sup>t</sup>* ≤ *tmk*) [6]. We will thus investigate the convergence of different simulators to this zero-toughness/small time solution.

For dimensionless time above unity (*τ* > 1 (*t* > *tmk*)), the solution transitions from the viscosity dominated (early-time) to the toughness (large-time) dominated regime. The toughness dominated regime is reached for K ≈ 3.5 (*τ* ≈ 70000). Note that, for infinitely large dimensionless toughness (i.e. zero viscosity), the solution is also self-similar [6] and is also denoted as the K-vertex solution. We will also briefly investigate, for a subset of the simulators considered, the transition between these two regimes of propagation (M to K), focusing mostly on fracture length versus time.

#### **2.3. Some practical numbers**

4 Effective and Sustainable Hydraulic Fracturing

**2.2. Radial hydraulic fracture**

of time. Introducing the characteristic time *tmk* =

*w* and net pressure *p* can be written as:

time *τ* = *t*/*tmk* we have (see [6, 14] for more details):

details). The tip region affected by the asymptote actually extends to about 10 to 20 percent

The growth of a radial hydraulic fracture spans both the viscosity and toughness regimes of propagation [6, 11]. At early times, the perimeter and the opening of the fracture are small and most of the energy is spent in viscous flow, whereas at a later times, the fracture perimeter and opening are larger and the fracture energy required to extend the fracture dominates the energy required to drive the viscous fluid through the fracture. The radial solution is also dependent only on a dimensionless toughness which in this case is a function

K = *τ*1/9

Solutions for the case of zero and infinite dimensionless toughness (i.e. small and large dimensionless time) have been obtained semi-analytically [6]. The complete transient solution can be obtained only numerically. A reference algorithm [7] based on an explicit moving mesh algorithm with proper matching of the multiscale HF tip asymptotics [10] will provide the baseline for the comparisons for intermediate times. The fracture radius *R*, width

*<sup>K</sup>*′<sup>4</sup> *<sup>γ</sup>*(*τ*)

*<sup>o</sup> <sup>µ</sup>*′1/2

where the dimensionless solution F = {*γ*, Ω, Π} depends only on dimensionless time *τ* =

As before, we will pay particular attention to the case of small dimensionless toughness, i.e. early-time, which is the most challenging numerically. In the limit of zero-toughness/early-time (i.e. viscosity dominated propagation, here refereed as the M-vertex), the solution is self-similar and can be conveniently written in the following

Π(*ρ*, *τ*)

*<sup>E</sup>*′*Qoµ*′

*<sup>R</sup>* <sup>=</sup> *<sup>E</sup>*′3*Qoµ*′

*<sup>p</sup>* <sup>=</sup> *<sup>K</sup>*′<sup>3</sup> *E*′3/2*Q*1/2

*w* =

*t*/*tmk* and scaled position *ρ* = *r*/*R* along the fracture.

viscosity scaling:

 *<sup>µ</sup>*′5*Q*<sup>3</sup> *oE*′<sup>13</sup> *K*′<sup>18</sup>

1/2

*<sup>K</sup>*′ <sup>Ω</sup>(*ρ*, *<sup>τ</sup>*) (2)

, and the dimensionless

of the plane-strain fracture length for dimensionless toughnesses below 0.5.

Although rock properties and stimulation practices vary, it is interesting to compute the scales and dimensionless numbers previously introduced. Let's assume a "tight" rock with the following realistic properties: a plane-strain Young's modulus of 40 GPa and a fracture toughness of 1.5 MPa.m1/2. First, for a hydraulic fracturing treatment using a highly viscous fluid (e.g. gel-like with *µ<sup>f</sup>* = 100 cPoise) at a practical rate of 10 Barrels per minute, we obtain a transition time scale *tmk* for a radial fracture of 4.2 10<sup>6</sup> seconds! The propagation of such a hydraulic fracture for a realistic injection duration (i.e. less than two hours) will always be in the viscosity dominated regime of propagation. Remember that the dimensionless toughness evolves as (*t*/*tmk*)1/9 for radial fractures. For the plane-strain geometry (where the injection rate is per meter in the out-of-plane direction), using the same parameters we obtain a dimensionless toughness of 0.26 clearly indicating a viscosity dominated propagation.

For a slick-water treatment, popular in shale-gas reservoirs, the injection rates are usually much higher in order to compensate for the low viscosity of water and to obtain a sufficiently wide fracture to accomodate proppant (see the scales in front of the opening *w* in Equations (1) and (3)). Using a value of 20 Barrels per minute (a realistic value for a single perforation cluster) and the viscosity of water, the radial transition time scale *tmk* now reduces to two minutes. For radial fractures, the toughness dominated regime of propagation is obtained for dimensionless toughness above 3.5 (see [6]), which corresponds to *t* 75000 *tmk*, which translate for this particulate case to *t* 250 hours. This indicates that most of the duration of the treatment will take place in the transition from the viscosity to the toughness regimes of propagation for a radial fracture (assuming that no stress barriers affect the fracture geometry). For the plane-strain fracture geometry, we obtain a dimensionless toughness of 0.3 (assuming the same injection rate per meter in the out-of-plane direction), for which the propagation is still dominated by viscosity.

10.5772/56212

861

cohesive traction is zero) and hence there will be no cohesive traction contribution, but fluid pressure is acting on the open fracture surfaces. So a coupled fluid pressure-traction-separation relationship exists between the cohesive zone defined by the traction-separation law and the pressurised fracture as found from solving the lubrication equation with the constraint that all tractions acting on the entire fracture and the cohesive zone must be in equilibrium. In this cohesive finite element model [19], the irreversible bilinear traction-separation cohesive law is adopted. An incompressible Newtonian fluid is injected at the centre of the fracture at constant injection rate. There is no fluid leak-off through the impermeable surfaces of the fracture, so only flow in the fracture radius direction is modelled. The cohesive elements at the injection point are defined as initially open to allow entry of fluid, and so that the initial flow and fracture growth is possible. Infinite elements surrounding the finite domain, which contains a hydraulic fracture, have been used to model the far-field boundary. Further details of the finite element model can

Accuracy and Convergence of Hydraulic Fracture Simulators

http://dx.doi.org/10.5772/56212

This code, also based on a fixed mesh, simulates straight hydraulic fractures in two-dimensions (plane-strain and axisymmetric fractures) using a fully implicit scheme to solve for the coupling between the elasticity equation (discretized using the displacement discontinuity method), the fluid conservation (discretized using a finite volume scheme) and to locate the fracture front. An increment of fracture length is given and the corresponding time-step (to reach the new fracture length) is solved by satisfying the fracture propagation condition in the tip element in a weak form: i.e. the volume of the tip element is enforced to be equal to the LEFM square-root asympote (The algorithm is similar to the one described in [20], see also [21]). The HF viscosity tip asymptote [8] can also be used for the case of low fracture toughness, its performance will be compared to the LEFM asymptote. Results obtained using the LEFM and the viscosity asymptotes will be denoted as 1DPlanarHF\_lefm, and 1DPlanarHF\_m respectively. All the results presented here use a grid with a constant element size (i.e. without any refinements), a

EMMA is an Explicit Moving Mesh Algorithm for radial geometry, which embeds the proper multiscale tip asymptotes of the hydraulic fracture depending on its velocity (see [7] for more details). It is extremely accurate and the moving mesh nature of the algorithm allows it to span more than ten orders of magnitude in dimensionless time. It notably provides a good solution for the transition between the viscosity and toughness

This algorithm [22] models the evolution a hydraulic fracture with an arbirarily shaped boundary that is assumed to propagate in a plane (a planar fracture in a 3D elastic medium), which is typically perpendicular to the minimum principal stress direction. The three dimensional elastic equilibrium equations are discretized using the dispacement discontinuity boundary integral method in which the fracture within the plane is represented by constant width rectangular elements that are collocated at element centres.

re-coarsening of the mesh during the simulation is possible.

dominated regime of propagation for a radial fracture.

**3.2. Three dimensional codes**

1. The Implicit Level Set Algorithm (ILSA)

be found in [19].

3. 1DPlanarHF

4. EMMA

These examples show the importance of the viscosity dominated regime of propagation for oil and gas hydraulic fracturing applications. Numerical simulators therefore need to be able to capture this regime of propagation, which is a difficult task especially if the algorithm relies solely on the linear elastic fracture mechanics propagation condition that manifests itself at a length scale near the fracture tip that is much smaller than the modeling lengthscale in the viscosity dominated case [10, 15].

#### **3. Simulators tested**

Two classes of simulators have been tested: codes simulating a two-dimensional configuration (plane-strain or/and axisymmetry) where the fracture is a one dimensional geometrical object, and three dimensional codes simulating planar fractures (which are two-dimensional objects in three dimensions). We now briefly describe the algorithms used by these different simulators.

#### **3.1. Two dimensional codes**

1. MineHF2D

This simulator (see [16] for more details) handles the propagation of both straight and curved hydraulic fractures in plane-strain. The algorithm is based on a fixed grid. It uses the displacement discontinuity method to solve the elastic equations coupled to a finite difference scheme for the fluid flow within the fracture. This algorithm includes the presence of a fluid lag at the fracture tip; for the simulated case reported here, a large confining stress value was used to minimize the fluid lag (see [17, 18] for more discussion on the effect of fluid lag). Explicit time-stepping is used and a volume of fluid method locates the fluid front. The fracture propagation criterion is based on the linear elastic fracture mechanics asymptote. The stress intensity factors are obtained using the displacement method with an adjusting factor of 0.88. A mesh with variable element sizes was used with refinement toward the fracture tip.

2. FEM\_Cohesive

This code is based on a finite element model and the pore pressure cohesive element implemented in Abaqus (Abaqus 6-10.2, 2010). It can handle both plane-strain and axi-symmetric configurations (e.g. plane-strain and radial hydraulic fractures). In this model, a pre-defined surface made up of elements that support the cohesive zone traction-separation calculation is embedded in the rock and the hydraulic fracture grows along this pre-defined surface. The fracture process zone (unbroken cohesive zone) is defined within the separating surfaces where the surface tractions are nonzero. The fracture is fully filled with fluid in the fully damaged cohesive zone (where the

cohesive traction is zero) and hence there will be no cohesive traction contribution, but fluid pressure is acting on the open fracture surfaces. So a coupled fluid pressure-traction-separation relationship exists between the cohesive zone defined by the traction-separation law and the pressurised fracture as found from solving the lubrication equation with the constraint that all tractions acting on the entire fracture and the cohesive zone must be in equilibrium. In this cohesive finite element model [19], the irreversible bilinear traction-separation cohesive law is adopted. An incompressible Newtonian fluid is injected at the centre of the fracture at constant injection rate. There is no fluid leak-off through the impermeable surfaces of the fracture, so only flow in the fracture radius direction is modelled. The cohesive elements at the injection point are defined as initially open to allow entry of fluid, and so that the initial flow and fracture growth is possible. Infinite elements surrounding the finite domain, which contains a hydraulic fracture, have been used to model the far-field boundary. Further details of the finite element model can be found in [19].

#### 3. 1DPlanarHF

6 Effective and Sustainable Hydraulic Fracturing

the propagation is still dominated by viscosity.

in the viscosity dominated case [10, 15].

**3. Simulators tested**

by these different simulators.

**3.1. Two dimensional codes**

was used with refinement toward the fracture tip.

1. MineHF2D

2. FEM\_Cohesive

cluster) and the viscosity of water, the radial transition time scale *tmk* now reduces to two minutes. For radial fractures, the toughness dominated regime of propagation is obtained for dimensionless toughness above 3.5 (see [6]), which corresponds to *t* 75000 *tmk*, which translate for this particulate case to *t* 250 hours. This indicates that most of the duration of the treatment will take place in the transition from the viscosity to the toughness regimes of propagation for a radial fracture (assuming that no stress barriers affect the fracture geometry). For the plane-strain fracture geometry, we obtain a dimensionless toughness of 0.3 (assuming the same injection rate per meter in the out-of-plane direction), for which

These examples show the importance of the viscosity dominated regime of propagation for oil and gas hydraulic fracturing applications. Numerical simulators therefore need to be able to capture this regime of propagation, which is a difficult task especially if the algorithm relies solely on the linear elastic fracture mechanics propagation condition that manifests itself at a length scale near the fracture tip that is much smaller than the modeling lengthscale

Two classes of simulators have been tested: codes simulating a two-dimensional configuration (plane-strain or/and axisymmetry) where the fracture is a one dimensional geometrical object, and three dimensional codes simulating planar fractures (which are two-dimensional objects in three dimensions). We now briefly describe the algorithms used

This simulator (see [16] for more details) handles the propagation of both straight and curved hydraulic fractures in plane-strain. The algorithm is based on a fixed grid. It uses the displacement discontinuity method to solve the elastic equations coupled to a finite difference scheme for the fluid flow within the fracture. This algorithm includes the presence of a fluid lag at the fracture tip; for the simulated case reported here, a large confining stress value was used to minimize the fluid lag (see [17, 18] for more discussion on the effect of fluid lag). Explicit time-stepping is used and a volume of fluid method locates the fluid front. The fracture propagation criterion is based on the linear elastic fracture mechanics asymptote. The stress intensity factors are obtained using the displacement method with an adjusting factor of 0.88. A mesh with variable element sizes

This code is based on a finite element model and the pore pressure cohesive element implemented in Abaqus (Abaqus 6-10.2, 2010). It can handle both plane-strain and axi-symmetric configurations (e.g. plane-strain and radial hydraulic fractures). In this model, a pre-defined surface made up of elements that support the cohesive zone traction-separation calculation is embedded in the rock and the hydraulic fracture grows along this pre-defined surface. The fracture process zone (unbroken cohesive zone) is defined within the separating surfaces where the surface tractions are nonzero. The fracture is fully filled with fluid in the fully damaged cohesive zone (where the This code, also based on a fixed mesh, simulates straight hydraulic fractures in two-dimensions (plane-strain and axisymmetric fractures) using a fully implicit scheme to solve for the coupling between the elasticity equation (discretized using the displacement discontinuity method), the fluid conservation (discretized using a finite volume scheme) and to locate the fracture front. An increment of fracture length is given and the corresponding time-step (to reach the new fracture length) is solved by satisfying the fracture propagation condition in the tip element in a weak form: i.e. the volume of the tip element is enforced to be equal to the LEFM square-root asympote (The algorithm is similar to the one described in [20], see also [21]). The HF viscosity tip asymptote [8] can also be used for the case of low fracture toughness, its performance will be compared to the LEFM asymptote. Results obtained using the LEFM and the viscosity asymptotes will be denoted as 1DPlanarHF\_lefm, and 1DPlanarHF\_m respectively. All the results presented here use a grid with a constant element size (i.e. without any refinements), a re-coarsening of the mesh during the simulation is possible.

#### 4. EMMA

EMMA is an Explicit Moving Mesh Algorithm for radial geometry, which embeds the proper multiscale tip asymptotes of the hydraulic fracture depending on its velocity (see [7] for more details). It is extremely accurate and the moving mesh nature of the algorithm allows it to span more than ten orders of magnitude in dimensionless time. It notably provides a good solution for the transition between the viscosity and toughness dominated regime of propagation for a radial fracture.

#### **3.2. Three dimensional codes**

#### 1. The Implicit Level Set Algorithm (ILSA)

This algorithm [22] models the evolution a hydraulic fracture with an arbirarily shaped boundary that is assumed to propagate in a plane (a planar fracture in a 3D elastic medium), which is typically perpendicular to the minimum principal stress direction. The three dimensional elastic equilibrium equations are discretized using the dispacement discontinuity boundary integral method in which the fracture within the plane is represented by constant width rectangular elements that are collocated at element centres. The Reynolds lubrication equation, expressing the conservations of mass of the viscous fluid contained within the crack surfaces, is discretized using a finite volume method also defined with respect to quantities sampled at the centres of the rectangular elements. At the periphery of the fracture, which may not conform to the structured rectangular mesh, the boundary is represented using a concept of partially filled tip elements that are used to define average fracture widths, which are also sampled at element centres. The distinguishing feature of this algorithm is its ability to locate the fracture free boundary using the asymptotic behavior of the hydraulic fracture width that is applicable at a particular point on the fracture perimeter. The free boundary is located by the following iterative process: given an initial guess for the fracture boundary *∂S*, determine the corresponding trial fracture width *w* and fluid pressure field *pf* ; in the ribbon of elements that are completely filled with fluid and, which share at least one side with a partially filled tip element, use the trial width values to estimate the distance to the free boundary by inverting the applicable tip asymptotic behavior [10]; use these estimates of the distance to the free boundary as initial conditions for the eikonal equation |∇*T*(*x*, *y*)| = 1, whose level set curve *T*(*x*, *y*) = 0 is the free boundary. The fracture boundary is then moved to the curve *T*(*x*, *y*) = 0 and the iterative process is repeated until convergence is achieved. The algorithm uses the multi-scale hydraulic fracture tip asymptotics solution [10] and thus automatically captures the different type of propagation regimes with relatively coarse mesh.

10.5772/56212

**Plane-strain Radial**

(K ≪ 1)

10−<sup>1</sup> < *τ* < 104 (.5 < K < 3.5)

http://dx.doi.org/10.5772/56212

Accuracy and Convergence of Hydraulic Fracture Simulators

◆◆◆◆◆◆◆◆ ◆◆◆◆◆◆◆ ◆◆◆◆◆◆ ◆◆◆◆◆ ◆◆◆◆ ◆◆◆ ◆◆◆ ◆◆◆ ◆◆ ◆◆ ◆◆

● MineHF2D

Simulators

■ FEM\_Cohesive

▲ 1DPlanarHF\_m

◆ 1DPlanarHF\_lefm

●●●●●● ● ●●●●●●●●●● ●●●●●●●●● ●●●●●●● ●●●●●● ●●●●● ●●●● ●●

▲▲▲▲▲▲▲▲▲▲▲ ▲▲▲▲▲▲▲▲ ▲▲▲▲▲▲ ▲▲▲▲▲ ▲▲▲▲ ▲▲▲▲ ▲▲▲ ▲▲ ▲▲ ▲▲

■■

◆

▲

863

<sup>K</sup> <sup>=</sup> 0.01 <sup>K</sup> <sup>=</sup> 0.1 *<sup>τ</sup>* <sup>≪</sup> <sup>1</sup>

MineHF2D n.a n.a

1DPlanarHF (\_lefm only) (\_lefm only) ILSA n.a n.a

■■■■ ■■■■ ■■■■■■■ ■■■■■■ ■■■■■ ■■■■■ ■■■■ ■■■■ ■■■■ ■■■ ■■■ ■■■ ■■■ ■■ ■■ ■■ ■■ ■■ ■■ ■■■■■ ■ ■ ■

0.001 0.005 0.010 0.050 0.100 0.500

The solution for a plane-strain hydraulic fracture driven by the injection of a Newtonian fluid at a constant rate is self-similar (i.e. evolves as a power-law of time). Our comparisons here focus on the case of viscosity dominated fractures (K < 0.5, spefically K = 0.1, 0.001) which are the most difficult to simulate numerically. The simulators tested for that configuration

Figure 1 displays, for the case K = 0.01 the dimensionless fracture opening profiles from the tip of the fracture obtained with the different simulators (at the last step of their simulations) as well as both the analytical solution and the viscosity HF tip asymptote [10]. Figure 2 is similar but for the case K = 0.1 . We first observe that the viscosity tip asymptote covers a region about 10 to 20 percent of the fracture from its tip. The different simulators provide

◆ ◆

▲

◆◆◆

▲ ▲

▲▲▲

FEM\_Cohesive

HFLattice n.a n.a

●●●● ●●●● ●●●●●●●●● ●●●● ●●●● ●●●●●●● ●●●●●● ●●●●●● ●●●●● ●●●●● ●●●● ●●●● ●●●● ●●● ●●● ●●● ●● ●● ●● ●● ●● ●●●●●● ● ● ●

▲

0.02

**4. Results and discussion**

0.05

0.10

0.20

0.50

1.00

◆

◆

**Figure 1.** Dimensionless opening Ω*<sup>m</sup>* from the fracture tip in log-log scale; plane-strain fracture K = 0.01.

are (Table 2): MineHF2D, FEM\_Cohesive, 1DPlanarHF\_lefm and 1DPlanarHF\_m.

**Table 2.** The benchmarks tested () for the different simulators.

Solution

m-tip asymptote

**4.1. The plane strain small toughness benchmark**

For this paper, a simplified version of the algorithm was also designed to only use the LEFM asymptote (hereafter denoted as ILSA\_lefm) for comparisons with other algorithms (MineHF2D, 1DPlanarHF etc.). In this version, we adapted the ILSA code to damp the front advance by rescaling the level set function *T*(*x*, *y*) so that the the maximum distance between any point in the ribbon elements and the damped free boundary is no more than three element lengths. This sequence of damped front positions enables the trial widths to be relaxed until fracture width profile presents a close approximation to the viscosity dominated solution, in spite of the fact that the tip elements are, by the nature of the ILSA\_lefm algorithm, locked into the LEFM asymptote.

2. HFLattice:

This code [23] simulates fracture propagation without limitation of shape, direction or number of fractures, as well as slip and opening along pre-existing joints. A 3D lattice formulation is used for simulation of deformation and fracturing. The lattice is a quasi-random assembly of nodes connected by non-linear shear and normal springs. The lattice resolution is given by the average node spacing. Newton's law of motion (for translation and rotation) is solved at the nodes using an explicit central difference scheme. The normal force in the spring is tested and a micro-crack is formed when breakage is detected (spring strength is adjusted to give the correct rock strength). A macro-fracture that develops in intact rock is thus characterized by an assembly of micro-cracks. Fluid flow and storage are based on a network of fluid nodes, located at broken springs or springs intersected by pre-existing joints, connected by pipes. The fluid network is updated continuously as new fracturing occurs. An explicit fluid pressure scheme is used to solve for fracture and matrix flow. The mechanical and flow models are fully coupled: fracture permeability depends on aperture (i.e. deformation of the mechanical components), fluid pressure affects deformation and strength of the solid model, and deformation of the solid model affects fluid pressures. In the algorithm, the lattice springs carry total forces, which affects force balance and motion. Also, effective stress is considered for joint slip or opening.


**Table 2.** The benchmarks tested () for the different simulators.

8 Effective and Sustainable Hydraulic Fracturing

relatively coarse mesh.

considered for joint slip or opening.

2. HFLattice:

The Reynolds lubrication equation, expressing the conservations of mass of the viscous fluid contained within the crack surfaces, is discretized using a finite volume method also defined with respect to quantities sampled at the centres of the rectangular elements. At the periphery of the fracture, which may not conform to the structured rectangular mesh, the boundary is represented using a concept of partially filled tip elements that are used to define average fracture widths, which are also sampled at element centres. The distinguishing feature of this algorithm is its ability to locate the fracture free boundary using the asymptotic behavior of the hydraulic fracture width that is applicable at a particular point on the fracture perimeter. The free boundary is located by the following iterative process: given an initial guess for the fracture boundary *∂S*, determine the corresponding trial fracture width *w* and fluid pressure field *pf* ; in the ribbon of elements that are completely filled with fluid and, which share at least one side with a partially filled tip element, use the trial width values to estimate the distance to the free boundary by inverting the applicable tip asymptotic behavior [10]; use these estimates of the distance to the free boundary as initial conditions for the eikonal equation |∇*T*(*x*, *y*)| = 1, whose level set curve *T*(*x*, *y*) = 0 is the free boundary. The fracture boundary is then moved to the curve *T*(*x*, *y*) = 0 and the iterative process is repeated until convergence is achieved. The algorithm uses the multi-scale hydraulic fracture tip asymptotics solution [10] and thus automatically captures the different type of propagation regimes with

For this paper, a simplified version of the algorithm was also designed to only use the LEFM asymptote (hereafter denoted as ILSA\_lefm) for comparisons with other algorithms (MineHF2D, 1DPlanarHF etc.). In this version, we adapted the ILSA code to damp the front advance by rescaling the level set function *T*(*x*, *y*) so that the the maximum distance between any point in the ribbon elements and the damped free boundary is no more than three element lengths. This sequence of damped front positions enables the trial widths to be relaxed until fracture width profile presents a close approximation to the viscosity dominated solution, in spite of the fact that the tip elements are, by the nature of the ILSA\_lefm algorithm, locked into the LEFM asymptote.

This code [23] simulates fracture propagation without limitation of shape, direction or number of fractures, as well as slip and opening along pre-existing joints. A 3D lattice formulation is used for simulation of deformation and fracturing. The lattice is a quasi-random assembly of nodes connected by non-linear shear and normal springs. The lattice resolution is given by the average node spacing. Newton's law of motion (for translation and rotation) is solved at the nodes using an explicit central difference scheme. The normal force in the spring is tested and a micro-crack is formed when breakage is detected (spring strength is adjusted to give the correct rock strength). A macro-fracture that develops in intact rock is thus characterized by an assembly of micro-cracks. Fluid flow and storage are based on a network of fluid nodes, located at broken springs or springs intersected by pre-existing joints, connected by pipes. The fluid network is updated continuously as new fracturing occurs. An explicit fluid pressure scheme is used to solve for fracture and matrix flow. The mechanical and flow models are fully coupled: fracture permeability depends on aperture (i.e. deformation of the mechanical components), fluid pressure affects deformation and strength of the solid model, and deformation of the solid model affects fluid pressures. In the algorithm, the lattice springs carry total forces, which affects force balance and motion. Also, effective stress is

**Figure 1.** Dimensionless opening Ω*<sup>m</sup>* from the fracture tip in log-log scale; plane-strain fracture K = 0.01.

#### **4. Results and discussion**

#### **4.1. The plane strain small toughness benchmark**

The solution for a plane-strain hydraulic fracture driven by the injection of a Newtonian fluid at a constant rate is self-similar (i.e. evolves as a power-law of time). Our comparisons here focus on the case of viscosity dominated fractures (K < 0.5, spefically K = 0.1, 0.001) which are the most difficult to simulate numerically. The simulators tested for that configuration are (Table 2): MineHF2D, FEM\_Cohesive, 1DPlanarHF\_lefm and 1DPlanarHF\_m.

Figure 1 displays, for the case K = 0.01 the dimensionless fracture opening profiles from the tip of the fracture obtained with the different simulators (at the last step of their simulations) as well as both the analytical solution and the viscosity HF tip asymptote [10]. Figure 2 is similar but for the case K = 0.1 . We first observe that the viscosity tip asymptote covers a region about 10 to 20 percent of the fracture from its tip. The different simulators provide 10 Effective and Sustainable Hydraulic Fracturing

● ●

● ●

1.0

Relative

 error

 on

2.0

1.5

3.0

5.0

in (%)

(in%)

1.0

Relative

 error

 on

2.0 3.0

1.5

10.0

in (%)

(in%)

15.0

5.0

7.0

● ●

**Figure 3.** Relative error in the fracture length as a function of the ratio mesh-size over fracture length K = 0.01.

■ ■

●

**Figure 4.** Relative error in the fracture length as a function of the ratio mesh-size over fracture length K = 0.1.

■■■■■■ ◆

■ ■■■

■ ■■ ■■■ ■

Figures 3 and 5 for K = 0.01 to Figures 4 and 6 for K = 0.1.

●

1 × 10- <sup>4</sup> 5 × 10- <sup>4</sup> 0.001 0.005 0.010 0.050 0.100

Similar observations can be made for the case of a dimesionless toughness K = 0.1, although all the algorithms using a fracture energy propagation condition perform slightly better due to the higher value of dimensionless toughness. In other other words, similar relative errors are obtained for larger values of *h*/ℓ (i.e. fewer elements), as can be seen by comparing

1 × 10- <sup>4</sup> 5 × 10- 40.001 0.005 0.010 0.050 0.100

■

◆ ◆◆ ◆◆ ◆◆◆ ◆◆◆◆◆◆ ◆◆◆◆ ◆◆◆◆ ◆◆◆◆◆ ◆◆◆◆◆◆

■ ■ ■ ■■ ■ ■■ ■■■ ■■

■■ ■■■

■

■

◆

◆ ◆◆ ◆◆ ◆◆◆ ◆◆◆◆◆◆

▲▲▲▲▲

Simulators

■ FEM\_Cohesive

▲ 1DPlanarHF\_m

◆ 1DPlanarHF\_lefm

● MineHF2D

▲▲▲▲

▲▲▲ ▲▲▲ ▲▲▲▲ ▲▲▲▲▲

Simulators

http://dx.doi.org/10.5772/56212

Accuracy and Convergence of Hydraulic Fracture Simulators

■ FEM\_Cohesive

▲ 1DPlanarHF\_m

◆ 1DPlanarHF\_lefm

● MineHF2D

▲▲▲▲ ▲▲▲▲▲ ▲▲▲▲▲

◆ ◆ ◆ ◆

10.5772/56212

865

**Figure 2.** Dimensionless opening Ω*<sup>m</sup>* from the fracture tip in log-log scale; plane-strain fracture K = 0.1.

width estimates that all correctly fall on the analytical solution "away" from the fracture tip. The distances from the tip at which the simulators recover the analytical solution appear to depend on both the mesh-size and the type of propagation condition used. For algorithms using the linear elastic fracture mechanics (LEFM) propagation condition (opening as a square root of the distance from the tip), this recovery distance from the tip is larger for coarser mesh sizes. The algorithm using a cohesive zone model appears to need significantly more refinement. In contrast, the algorithm using the viscosity HF tip asymptote (i.e. 1DPlanarHF\_m) is able to capture the fracture opening exactly all the way to the tip and with a much coarser mesh than was used for the other computations.

This dependence of the convergence toward the exact solution on the mesh size and propagation condition can be further observed in Figures 3-5, which display the rate of convergence for the fracture length and fracture volume as a function of the ratio of the mesh size over fracture length (i.e. the inverse of the number of elements to discretize the frature for uniform mesh). All simulators converge correctly toward the analytical solution but at very different computational costs. We can see that for algorithms using the LEFM condition or a cohesive zone model, the mesh size required to reach the same level of accuracy is about 20 times smaller than for the algorithm that uses the correct HF viscosity tip asymptote. In the case K = 0.01, 1DPlanarHF\_lefm needs about 400 elements (*h*/ ∼ .0025) to obtain a relative error of about 1 and 3 percent in the fracture length and fracture volume respectively, while smaller relative errors are already obtained when using 20 elements (*h*/ ∼ .05) for 1DPlanarHF\_m. The cost is even greater when a cohesive zone model is used: about 3000 elements (*h*/ <sup>∼</sup> <sup>3</sup> <sup>×</sup> <sup>10</sup>−4) are needed to reach a relative error of 1.5 and three percent in the fracture length and fracture volume respectively.

**Figure 3.** Relative error in the fracture length as a function of the ratio mesh-size over fracture length K = 0.01.

10 Effective and Sustainable Hydraulic Fracturing





 -





 



width estimates that all correctly fall on the analytical solution "away" from the fracture tip. The distances from the tip at which the simulators recover the analytical solution appear to depend on both the mesh-size and the type of propagation condition used. For algorithms using the linear elastic fracture mechanics (LEFM) propagation condition (opening as a square root of the distance from the tip), this recovery distance from the tip is larger for coarser mesh sizes. The algorithm using a cohesive zone model appears to need significantly more refinement. In contrast, the algorithm using the viscosity HF tip asymptote (i.e. 1DPlanarHF\_m) is able to capture the fracture opening exactly all the way to the tip and

This dependence of the convergence toward the exact solution on the mesh size and propagation condition can be further observed in Figures 3-5, which display the rate of convergence for the fracture length and fracture volume as a function of the ratio of the mesh size over fracture length (i.e. the inverse of the number of elements to discretize the frature for uniform mesh). All simulators converge correctly toward the analytical solution but at very different computational costs. We can see that for algorithms using the LEFM condition or a cohesive zone model, the mesh size required to reach the same level of accuracy is about 20 times smaller than for the algorithm that uses the correct HF viscosity tip asymptote. In the case K = 0.01, 1DPlanarHF\_lefm needs about 400 elements (*h*/ ∼ .0025) to obtain a relative error of about 1 and 3 percent in the fracture length and fracture volume respectively, while smaller relative errors are already obtained when using 20 elements (*h*/ ∼ .05) for 1DPlanarHF\_m. The cost is even greater when a cohesive zone model is used: about 3000 elements (*h*/ <sup>∼</sup> <sup>3</sup> <sup>×</sup> <sup>10</sup>−4) are needed to reach a relative error of 1.5 and three percent in the

 

 











**Figure 2.** Dimensionless opening Ω*<sup>m</sup>* from the fracture tip in log-log scale; plane-strain fracture K = 0.1.

with a much coarser mesh than was used for the other computations.

fracture length and fracture volume respectively.



 




 



(in%)

(in%)


**Figure 4.** Relative error in the fracture length as a function of the ratio mesh-size over fracture length K = 0.1.

Similar observations can be made for the case of a dimesionless toughness K = 0.1, although all the algorithms using a fracture energy propagation condition perform slightly better due to the higher value of dimensionless toughness. In other other words, similar relative errors are obtained for larger values of *h*/ℓ (i.e. fewer elements), as can be seen by comparing Figures 3 and 5 for K = 0.01 to Figures 4 and 6 for K = 0.1.

867

**4.2. The radial hydraulic fracture benchmark**

version ILSA\_lefm and the HFLattice model (Table 2).

**Table 3.** Range of dimensionless time of the simulations for the radial benchmark

also exhibits convergence as *h*/*R* decreases.

Owing to its previously argued importance in practice, we focus on early-time, where the relevant analytical solution corresponds to the case of zero-toughness. The simulators tested for this geometry and regime of propagation are the two-dimensional codes under axi-symmetry; FEM\_Cohesive and 1DPlanarHF\_lefm, and the 3D codes ILSA, the simplified

Accuracy and Convergence of Hydraulic Fracture Simulators

http://dx.doi.org/10.5772/56212

The different simulators have been run for different ranges of dimensionless time all within the viscosity dominated regime (*τ* < 1), and are described in Table 3. The HFLattice

> **FEM\_Cohesive 1DPlanarHF\_lefm ILSA/ILSA\_lefm** *<sup>τ</sup>* ∈ [10−<sup>10</sup> − <sup>2</sup> × <sup>10</sup>−8] *<sup>τ</sup>* ∈ [10−<sup>7</sup> − <sup>10</sup>−1] *<sup>τ</sup>* ∈ [10−<sup>18</sup> − <sup>10</sup>−17]

simulator has been run for the specific case of zero toughness (in fact zero tensile strength as

The convergence toward the zero-toughness solution as a function of the mesh size can be observed in Figure 7 (convergence of the fracture radius). Similar trends to the plane-strain case can be observed. The convergence requires a much finer mesh for the simulators using a cohesive zone model (axi-symmetric FEM\_Cohesive) and the LEFM propagation condition (axisymmetric 1DPlanarHF and ILSA\_lefm). ILSA, the only simulator using the appropriate HF tip asymptote, achieves the same accuracy with a much coarser mesh compared to all the other simulators. In particular, the version ILSA\_lefm, which uses the LEFM asymptote, needs about an order of magnitude finer mesh compared to ILSA for the same relative error. The FEM\_Cohesive algorithm requires a ratio *h*/*R* about two to three orders of magnitude smaller than ILSA for the same relative error. The HFLattice model, although less accurate,

The openings at the last time step of the simulation (refer to Table 3 for the corresponding dimensionless time) in the tip coordinate system are compared to the zero-toughness analytical solution in Figure 8. The FEM\_Cohesive algorithm captures the solution away from the fracture tip well, i.e at distance larger than 3-5 % of the fracture radius from the tip. Closer to the tip, the opening from the cohesive zone algorithm appears to slightly overestimate the fracture opening. The results of the algorithms using the LEFM propagation condition (1DPlanarHF\_lefm, ILSA\_lefm) converge toward the analytical opening as their mesh gets finer. On the other hand, ILSA exactly matches the analytical solution for the opening with a relatively coarse mesh (the last element of ILSA corresponds to a partially fractured element for which the fracture width is also function of the fracture front location

A comparison of the net pressure profile obtained by the different simulators and the analytical solution is displayed on Figure 9. Similar trends to the fracture opening can be observed. One has to note that the HFLattice simulator approximates the point-source by a finite volume, hence the discrepancy on the net pressure with respect to the point-source

*4.2.1. Viscosity dominated regime*

per the formulation).

within the element).

solution close to the fracture inlet.

**Figure 5.** Relative error in the dimensionless fracture volume as a function of mesh size over fracture length - Plane-strain fracture K = 0.01.

**Figure 6.** Relative error in the dimensionless fracture volume as a function of mesh size over fracture length - Plane-strain fracture K = 0.1.

### **4.2. The radial hydraulic fracture benchmark**

#### *4.2.1. Viscosity dominated regime*

12 Effective and Sustainable Hydraulic Fracturing

● ●

● ●

■ ■

■ ■

1.0

5.0

2.0

Relative

fracture K = 0.1.

 error on

Volume (%)

1.5

3.0

2.0

Relative

fracture K = 0.01.

 error on

1.5

3.0

10.0

Volume (%)

20.0

15.0

5.0

7.0

■ ■

■ ■

●●

1 × 10- <sup>4</sup> 5 × 10- 40.001 0.005 0.010 0.050 0.100

**Figure 5.** Relative error in the dimensionless fracture volume as a function of mesh size over fracture length - Plane-strain

●

●

1 × 10- <sup>4</sup> 5 × 10- <sup>4</sup> 0.001 0.005 0.010 0.050 0.100

**Figure 6.** Relative error in the dimensionless fracture volume as a function of mesh size over fracture length - Plane-strain

◆

◆

◆ ◆ ◆◆ ◆◆◆◆◆

◆ ◆ ◆◆ ◆◆ ◆◆ ◆◆ ◆◆◆ ◆◆◆ ◆◆◆ ◆◆◆◆ ◆◆◆◆ ◆◆◆◆◆◆◆

> ▲▲ ▲▲▲

▲▲▲▲ ▲

Simulators

■ FEM\_Cohesive

▲ 1DPlanarHF\_m

◆ 1DPlanarHF\_lefm

Simulators

■ FEM\_Cohesive

▲ 1DPlanarHF\_m

◆ 1DPlanarHF\_lefm

● MineHF2D

● MineHF2D

Owing to its previously argued importance in practice, we focus on early-time, where the relevant analytical solution corresponds to the case of zero-toughness. The simulators tested for this geometry and regime of propagation are the two-dimensional codes under axi-symmetry; FEM\_Cohesive and 1DPlanarHF\_lefm, and the 3D codes ILSA, the simplified version ILSA\_lefm and the HFLattice model (Table 2).

The different simulators have been run for different ranges of dimensionless time all within the viscosity dominated regime (*τ* < 1), and are described in Table 3. The HFLattice


**Table 3.** Range of dimensionless time of the simulations for the radial benchmark

simulator has been run for the specific case of zero toughness (in fact zero tensile strength as per the formulation).

The convergence toward the zero-toughness solution as a function of the mesh size can be observed in Figure 7 (convergence of the fracture radius). Similar trends to the plane-strain case can be observed. The convergence requires a much finer mesh for the simulators using a cohesive zone model (axi-symmetric FEM\_Cohesive) and the LEFM propagation condition (axisymmetric 1DPlanarHF and ILSA\_lefm). ILSA, the only simulator using the appropriate HF tip asymptote, achieves the same accuracy with a much coarser mesh compared to all the other simulators. In particular, the version ILSA\_lefm, which uses the LEFM asymptote, needs about an order of magnitude finer mesh compared to ILSA for the same relative error. The FEM\_Cohesive algorithm requires a ratio *h*/*R* about two to three orders of magnitude smaller than ILSA for the same relative error. The HFLattice model, although less accurate, also exhibits convergence as *h*/*R* decreases.

The openings at the last time step of the simulation (refer to Table 3 for the corresponding dimensionless time) in the tip coordinate system are compared to the zero-toughness analytical solution in Figure 8. The FEM\_Cohesive algorithm captures the solution away from the fracture tip well, i.e at distance larger than 3-5 % of the fracture radius from the tip. Closer to the tip, the opening from the cohesive zone algorithm appears to slightly overestimate the fracture opening. The results of the algorithms using the LEFM propagation condition (1DPlanarHF\_lefm, ILSA\_lefm) converge toward the analytical opening as their mesh gets finer. On the other hand, ILSA exactly matches the analytical solution for the opening with a relatively coarse mesh (the last element of ILSA corresponds to a partially fractured element for which the fracture width is also function of the fracture front location within the element).

A comparison of the net pressure profile obtained by the different simulators and the analytical solution is displayed on Figure 9. Similar trends to the fracture opening can be observed. One has to note that the HFLattice simulator approximates the point-source by a finite volume, hence the discrepancy on the net pressure with respect to the point-source solution close to the fracture inlet.

(in%)

10.5772/56212

●●●● ●● ●● ●● ● ●● ●● ● ● ●

▼▼▼ ▼▼ ▼▼ ▼▼ ▼▼ ▼ ▼ ▼

◆◆◆◆◆ ◆◆◆◆◆◆◆◆ ◆◆◆◆◆◆ ◆◆◆◆◆ ◆◆◆◆ ◆◆◆ ◆◆◆ ◆◆◆ ◆◆ ◆◆

869

▲▲▲▲▲▲▲▲ ▲▲▲ ▲▲ ▲▲ ▲▲ ▲ ▲ ▲

▼▼▼▼▼▼▼▼▼▼▼▼▼▼

■■■■ ■■■ ■■■■ ■■ ■■■■ ■ ■ ■ ■

http://dx.doi.org/10.5772/56212

Simulators

● ●

Accuracy and Convergence of Hydraulic Fracture Simulators

● HFLattice

▲ ILSA\_lefm

▼ ILSA

■ FEM\_Cohesive

◆ 1DPlanarHF\_lefm

■

◆◆ ◆◆◆◆ ◆

▲

●

■

●

◆

0.001 0.005 0.010 0.050 0.100 0.500 1.000

**Figure 8.** Fracture opening in the tip coordinates system (log-log scale) for a radial hydraulic fractures propagating in the viscosity dominated regime. Results from the different simulators for a dimensionless time *τ* < 10−<sup>6</sup> . Note that for better

Finally, we investigate the performance of a subset of the algorithms on the transition of the solution toward the toughness dominated solution. Such a transition typically happens between *τ* = 1 (K = 1) and *τ* ∼ 70000 (K = 3.4). The Explicit Moving Mesh Algorithm (EMMA), the 1DPlanarHF\_lefm and the ILSA codes are compared, focusing on the evolution of the dimensionless fracture radius with time. Figure 10 display the results. The fracture radius have been averaged over the fracture footprint for the results of ILSA (3D code), while

We clearly see that the different algorithms capture the transition between the viscosity and the toughness propagation regimes extremely well. They are virtually indistinguishable.

In this paper, we have investigated the performance of a number of hydraulic fracture propagation algorithms based on different propagation conditions: LEFM, cohesive zone model/tensile strength and algorithms accounting for the multi-scale nature of hydraulic fracture propagation in the near-tip region. This exercise was made possible thanks to the existence of analytical solutions for both geometries of hydraulic fractures. All the simulators investigated here are able to capture the analytical solutions of the different configurations

▼

■■■■■■ ■■■■■■ ■■■■■ ■■■■■ ■■■■■ ■■■■ ■■■■ ■■■■ ■■■■ ■■■ ■■■ ■■■ ■■■ ■■ ■■ ■■ ■■ ■■ ■ ■ ■ ■ ■ ■ ■

clarity of the plot, all the mesh points are not displayed for 1 − *ρ* > 0.04 for FEM\_Cohesive.

*4.2.2. Transition toward the toughness dominated regime*

the other codes are axi-symmetric by assumption.

tested, but at widely different computational costs.

0.02

**5. Conclusions**

0.05

0.10

0.20

0.50

1.00 M-solution

◆

■

◆ ◆

▼

▲

**Figure 7.** Evolution of the relative error in the fracture radius as a function of the mesh size for the different simulators radial fracture in the viscosity dominated regime (i.e. zero toughness / early time, *τ* < 1). All simulations displayed here are for *τ* < 10−6.

It is interesting to investigate more closely how an algorithm using the LEFM asymptote (e.g. ILSA\_lefm) is able to converge to the zero-toughness analytical solution. Consider the case of modeling a radial hydraulic fracture starting at a very small initial time (*τ* = <sup>10</sup><sup>−</sup>18), which corresponds to a dimensionless toughness K = 0.01, and is therefore very close to the M-Vertex (zero-toughness solution). If the LEFM asymptote is used in ILSA at this time, the asymptote would dictate that the fracture front needs to be advanced by roughly 104 element lengths! As already mentioned, to circumvent this problem, the ILSA\_lefm code damps the front advance by rescaling the level set function *T*(*x*, *y*) so that the the maximum distance between any point in the ribbon elements and the damped free boundary is no more than three element lengths. This sequence of damped front positions enables the trial widths to be relaxed until the fracture width profile presents a close approximation to the viscosity dominated solution, in spite of the fact that the tip elements are locked into the LEFM. Indeed, the algorithm settles on a solution in which the LEFM tip widths are very small and contribute very little to the net volume of the fracture, while over the remainder of the fracture away from the leading edge, the widths are locked into the viscous asymptote as dictated by the conservation of volume. Thus the damped front advances continue until the conservation of fluid volume dictates that it should stop at which time it approximates the viscous asymptote. Results from the standard ILSA code with the correct viscous asymptote and damped ILSA\_lefm code with the LEFM asymptote indicates that in order to obtain similar relative errors compared to the analytical solution, a mesh size that is an order of magnitude smaller is needed for ILSA\_lefm compared to the standard ILSA code (see Figure 7).

**Figure 8.** Fracture opening in the tip coordinates system (log-log scale) for a radial hydraulic fractures propagating in the viscosity dominated regime. Results from the different simulators for a dimensionless time *τ* < 10−<sup>6</sup> . Note that for better clarity of the plot, all the mesh points are not displayed for 1 − *ρ* > 0.04 for FEM\_Cohesive.

#### *4.2.2. Transition toward the toughness dominated regime*

Finally, we investigate the performance of a subset of the algorithms on the transition of the solution toward the toughness dominated solution. Such a transition typically happens between *τ* = 1 (K = 1) and *τ* ∼ 70000 (K = 3.4). The Explicit Moving Mesh Algorithm (EMMA), the 1DPlanarHF\_lefm and the ILSA codes are compared, focusing on the evolution of the dimensionless fracture radius with time. Figure 10 display the results. The fracture radius have been averaged over the fracture footprint for the results of ILSA (3D code), while the other codes are axi-symmetric by assumption.

We clearly see that the different algorithms capture the transition between the viscosity and the toughness propagation regimes extremely well. They are virtually indistinguishable.

#### **5. Conclusions**

14 Effective and Sustainable Hydraulic Fracturing

■ ■ ■ ■ ■■ ■■ ■■■■ ■■■

Simulators

● HFLattice

▲ ILSA\_lefm

1 × 10 - <sup>4</sup> 5 × 10 - <sup>4</sup> 0.001 0.005 0.010 0.050 0.100

**Figure 7.** Evolution of the relative error in the fracture radius as a function of the mesh size for the different simulators radial fracture in the viscosity dominated regime (i.e. zero toughness / early time, *τ* < 1). All simulations displayed here are for

It is interesting to investigate more closely how an algorithm using the LEFM asymptote (e.g. ILSA\_lefm) is able to converge to the zero-toughness analytical solution. Consider the case of modeling a radial hydraulic fracture starting at a very small initial time (*τ* = <sup>10</sup><sup>−</sup>18), which corresponds to a dimensionless toughness K = 0.01, and is therefore very close to the M-Vertex (zero-toughness solution). If the LEFM asymptote is used in ILSA at this time, the asymptote would dictate that the fracture front needs to be advanced by roughly 104 element lengths! As already mentioned, to circumvent this problem, the ILSA\_lefm code damps the front advance by rescaling the level set function *T*(*x*, *y*) so that the the maximum distance between any point in the ribbon elements and the damped free boundary is no more than three element lengths. This sequence of damped front positions enables the trial widths to be relaxed until the fracture width profile presents a close approximation to the viscosity dominated solution, in spite of the fact that the tip elements are locked into the LEFM. Indeed, the algorithm settles on a solution in which the LEFM tip widths are very small and contribute very little to the net volume of the fracture, while over the remainder of the fracture away from the leading edge, the widths are locked into the viscous asymptote as dictated by the conservation of volume. Thus the damped front advances continue until the conservation of fluid volume dictates that it should stop at which time it approximates the viscous asymptote. Results from the standard ILSA code with the correct viscous asymptote and damped ILSA\_lefm code with the LEFM asymptote indicates that in order to obtain similar relative errors compared to the analytical solution, a mesh size that is an order of magnitude smaller is needed for ILSA\_lefm compared to the standard ILSA code

▼ ILSA

■ FEM\_Cohesive ◆ 1DPlanarHF\_lefm

■■ ■ ■■ ■■ ■■■ ■■ ■■■ ■ ■ ■ ■ ■ ■ ■■■■ ■ ■■ ■ ■■ ■ ■ ■ ■■ ■

■

■ ■

■■

0.01

0.1

Relative

*τ* < 10−6.

(see Figure 7).

 error

 on

in (%)

(in%)

1

10

■ ■

■■ ■ ■■

■ ■ ■ ■

■ ■ ■ ■ ■

■ ■■■ ■ ■ ■ ■ ■ ■ ■ ■ ■■ ■ ■■■ ■ ■ ■ ■■■■ ■ ■■ ■ ■■ ■■■■ ■■ ■ ■

● ● ●

▲▲▲ ▼

▲ ▲▲▲▲ ▲▲▲ ▲▲ ▲ ▲ ▲▲▲ ▲▲ ▲▲▲▲▲▲ ▲▲▲▲ ▲▲ ▲▲▲ ▲▲ ▲▲▲

▼ ▼▼▼▼▼ ▼▼▼▼▼▼ ▼▼ ▼▼ ▼▼ ▼▼ ▼▼▼▼ ▼▼ ▼ ▼ ▼▼▼ ▼▼

◆◆◆◆◆◆◆◆◆◆◆ ◆

▼ ▼ ▼ ▼▼ ▼

> In this paper, we have investigated the performance of a number of hydraulic fracture propagation algorithms based on different propagation conditions: LEFM, cohesive zone model/tensile strength and algorithms accounting for the multi-scale nature of hydraulic fracture propagation in the near-tip region. This exercise was made possible thanks to the existence of analytical solutions for both geometries of hydraulic fractures. All the simulators investigated here are able to capture the analytical solutions of the different configurations tested, but at widely different computational costs.

871

Algorithms based on the Linear Elastic Fracture Mechanics propagation condition requires a fine mesh (*h*/ℓ <sup>10</sup>−<sup>2</sup> − <sup>10</sup><sup>−</sup>3) in order to capture viscosity dominated hydraulic fracture propagation. A fine mesh is needed for these algorithms to capture the viscosity opening asymptote in order to properly match the fracture volume. Cohesive zone models, which model the fracture process zone, require even finer meshes. This is due to the fact that the cohesive zone length-scale is even smaller than that of the region of influence of the linear elastic fracture mechanics (LEFM) square-root near-tip asymptote. By contrast, when the algorithms use the appropriate multi-scale hydraulic fracture asymptote in the near-tip region, the exact solution can be matched with a very coarse mesh (i.e. *<sup>h</sup>*/ℓ ≈ <sup>10</sup><sup>−</sup>1). Extremely efficient and fast propagation algorithms can thus be developed with even better accuracy than algorithms based on the classical LEFM propagation condition. Computational cost and accuracy may not be the only concern when developing a simulator. Algorithm flexibility may also be important. We hope that the study reported in this paper can help in

Accuracy and Convergence of Hydraulic Fracture Simulators

http://dx.doi.org/10.5772/56212

Finally, we also would like to advocate that the different analytical reference solutions used in this paper be used as a minimal series of benchmarks for any hydraulic fracturing simulator.

Brice Lecampion1,⋆, Anthony Peirce2, Emmanuel Detournay3, Xi Zhang4, Zuorong Chen4, Andrew Bunger4, Christine Detournay5, John Napier, Safdar Abbas1, Dmitry Garagash6 and

[1] R. Carbonell, J. Desroches, and E. Detournay. A comparison between a semi-analytical and a numerical solution of a two-dimensional hydraulic fracture. *Int. J. Sol. Struct.*,

[2] J. Adachi and E. Detournay. Self-similar solution of a plane-strain fracture driven by a power-law fluid. *International Journal for Numerical and Analytical Methods in*

[3] D.I. Garagash and E. Detournay. Plane-strain propagation of a fluid-driven fracture: Small toughness solution. *ASME J. Appl. Mech.*, 72:916–928, November 2005.

[4] D. I. Garagash. Plane-strain propagation of a fluid-driven fracture during injection and shut-in: Asymptotics of large toughness. *Engineering Fracture Mechanics*, 73(4):456–481,

2 Department of Mathematics, University of British Columbia, Vancouver, Canada 3 Department of Civil Engineering, University of Minnesota, Minneapolis, USA 4 CSIRO Earth Science and Resource Engineering, Clayton, Victoria, Australia

6 Department of Civil and Resource Engineering, Dalhousie University, Canada

making an educated choice of algorithm.

<sup>⋆</sup> Address all correspondence to: BLecampion@slb.com

1 Schlumberger Doll Research, Cambridge, USA

5 Itasca Consulting Group, Minneapolis, USA

*Geomechanics*, 26(6):579–604, May 2002.

36(31):4869–4888, 1999.

March 2006.

**Author details**

Peter Cundall5

**References**

**Figure 9.** Dimensionless net pressure for a radial hydraulic fracture in the viscosity dominated regime (i.e. zero-toughness / early-time). Comparisons of the simulators (for a dimensionless time *τ* < 10−6) with the zero-toughness solution. The results of FEM\_Cohesive and 1DPlanarHF\_lefm are plotted only every 50 and 2 mesh points respectively for clarity.

**Figure 10.** Radial hydraulic fracture - Transition between the viscosity dominated regime at early time to Toughness dominated regime at large time. Comparisons between different simulators and the zero (blue) and infinite toughness (red) solution.

Algorithms based on the Linear Elastic Fracture Mechanics propagation condition requires a fine mesh (*h*/ℓ <sup>10</sup>−<sup>2</sup> − <sup>10</sup><sup>−</sup>3) in order to capture viscosity dominated hydraulic fracture propagation. A fine mesh is needed for these algorithms to capture the viscosity opening asymptote in order to properly match the fracture volume. Cohesive zone models, which model the fracture process zone, require even finer meshes. This is due to the fact that the cohesive zone length-scale is even smaller than that of the region of influence of the linear elastic fracture mechanics (LEFM) square-root near-tip asymptote. By contrast, when the algorithms use the appropriate multi-scale hydraulic fracture asymptote in the near-tip region, the exact solution can be matched with a very coarse mesh (i.e. *<sup>h</sup>*/ℓ ≈ <sup>10</sup><sup>−</sup>1). Extremely efficient and fast propagation algorithms can thus be developed with even better accuracy than algorithms based on the classical LEFM propagation condition. Computational cost and accuracy may not be the only concern when developing a simulator. Algorithm flexibility may also be important. We hope that the study reported in this paper can help in making an educated choice of algorithm.

Finally, we also would like to advocate that the different analytical reference solutions used in this paper be used as a minimal series of benchmarks for any hydraulic fracturing simulator.

#### **Author details**

16 Effective and Sustainable Hydraulic Fracturing

▲

Simulators

● HFLattice

▲ ILSA\_lefm

▼ ILSA

0

5

10

15

M solution

20


0.0

0.5

1.0

■ FEM\_Cohesive

◆ 1DPlanarHF\_lefm

◆

▲

▼ ▼

● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●

▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼▼ ▼ ▼

◆

◆◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆

0.0 0.2 0.4 0.6 0.8 1.0 - 1.0

**Figure 9.** Dimensionless net pressure for a radial hydraulic fracture in the viscosity dominated regime (i.e. zero-toughness / early-time). Comparisons of the simulators (for a dimensionless time *τ* < 10−6) with the zero-toughness solution. The results

K solution

0 500 1000 1500 2000 2500 3000 3500

**Figure 10.** Radial hydraulic fracture - Transition between the viscosity dominated regime at early time to Toughness dominated regime at large time. Comparisons between different simulators and the zero (blue) and infinite toughness (red) solution.

of FEM\_Cohesive and 1DPlanarHF\_lefm are plotted only every 50 and 2 mesh points respectively for clarity.

■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■

▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲▲

M-solution

●

◆ ◆

▲

▼ ▼ ▼ ▼

▼ ◆

Simulators ● 1DPlanarHF\_lefm

ILSA

EMMA

■ ■

▲ ▲

●

◆

▼

■

Brice Lecampion1,⋆, Anthony Peirce2, Emmanuel Detournay3, Xi Zhang4, Zuorong Chen4, Andrew Bunger4, Christine Detournay5, John Napier, Safdar Abbas1, Dmitry Garagash6 and Peter Cundall5

<sup>⋆</sup> Address all correspondence to: BLecampion@slb.com

1 Schlumberger Doll Research, Cambridge, USA


5 Itasca Consulting Group, Minneapolis, USA

6 Department of Civil and Resource Engineering, Dalhousie University, Canada

### **References**


[5] J. Hu and D. I. Garagash. Plane-Strain Propagation of a Fluid-Driven Crack in a Permeable Rock with Fracture Toughness. *Journal of Engineering Mechanics*, 136(9):1152, 2010.

10.5772/56212

873

[20] E. Gordeliy and Emmanuel Detournay. A fixed grid algorithm for simulating the propagation of a shallow hydraulic fracture with a fluid lag. *International Journal for*

Accuracy and Convergence of Hydraulic Fracture Simulators

http://dx.doi.org/10.5772/56212

[21] B. Lecampion. Hydraulic fracture initiation from an open-hole: Wellbore size, pressurization rate and fluid-solid coupling effects. In *46th U.S. Rock*

[22] A. P. Peirce and E. Detournay. An implicit level set method for modeling hydraulically driven fractures. *Computer Methods in Applied Mechanics and Engineering*,

[23] B. Damjanac, C. Detournay, P. Cundall, and Varun. Three-Dimensional Numerical Model of Hydraulic Fracturing in Fractured Rock Masses. In A. Bunger, R.G. Jeffrey, and J. McLennan, editors, *Effective and Sustainable Hydraulic Fracturing*, Brisbane, Australia,

*Mechanics/Geomechanics Symposium*, number ARMA 2012-601, June 24–27 2012.

197(33-40):2858–2885, June 2008.

20–22 May 2013. Intech.

*Numerical and Analytical Methods in Geomechanics*, 35(5):602–629, April 2011.


[20] E. Gordeliy and Emmanuel Detournay. A fixed grid algorithm for simulating the propagation of a shallow hydraulic fracture with a fluid lag. *International Journal for Numerical and Analytical Methods in Geomechanics*, 35(5):602–629, April 2011.

18 Effective and Sustainable Hydraulic Fracturing

39(26):6311–6337, December 2002.

medium. *J. Appl. Mech.*, 67:183–192, 2000.

fracturing offset rocks". *J. Geoph. Res.*, 81:6292, 1976.

*International Journal of Geomechanics*, 4(1):35, 2004.

early-time solution. *Int. J. Sol. Struct.*, 43:5811–5835, 2006.

*of Petroleum Science and Engineering*, 88–89:136–144, 2012.

University of Minnesota, 2003.

2010.

1955.

29:1317–1340, 2005.

[5] J. Hu and D. I. Garagash. Plane-Strain Propagation of a Fluid-Driven Crack in a Permeable Rock with Fracture Toughness. *Journal of Engineering Mechanics*, 136(9):1152,

[6] A. Savitski and E. Detournay. Propagation of a penny-shaped fluid-driven fracture in an impermeable rock: asymptotic solutions. *International Journal of Solids and Structures*,

[7] M. Madyarova. Fluid-driven penny-shaped fracture in permeable rock. Master's thesis,

[8] J. Desroches, E. Detournay, B. Lenoach, P. Papanastasiou, JRA Pearson, M. Thiercelin, and A. Cheng. The crack tip region in hydraulic fracturing. *Proceedings of the Royal*

[9] D.I. Garagash and E. Detournay. The tip region of a fluid-driven fracture in an elastic

[10] D. I. Garagash, E. Detournay, and J. Adachi. Multiscale tip asymptotics in hydraulic fracture with leak-off. *Journal of Fluid Mechanics*, 669:260–297, February 2011.

[11] H. Abé, L. M. Keer, and T. Mura. "growth rate of a penny-shaped crack in hydraulic

[12] S. Khristianovic and Y. Zheltov. Formation of vertical fractures by means of highly viscous fluids. In *Proc., 4th World Petroleum Congress*, volume II, pages 579–586, Rome,

[13] J. Geertsma and F. De Klerk. A rapid method of predicting width and extent of hydraulically induced fractures. *Journal of Petroleum Technology*, 21(12):1571–1581, 1969.

[14] E. Detournay. Propagation regimes of fluid-driven fractures in impermeable rocks.

[15] D. I. Garagash. Scaling of Physical Processes in Fluid-Driven Fracture: Perspective from the Tip. In F.M. Borodich, editor, *IUTAM Symposium on Scaling in Solid Mechanics*,

[16] X. Zhang, R. G. Jeffrey, and E. Detournay. Propagation of a fluid-driven fracture parallel to the free surface of an elastic half plane. *Int. J. Numer. Anal. Meth. Geomech.*,

[17] D. Garagash. Propagation of a plane-strain fluid-driven fracture with a fluid lag:

[18] B. Lecampion and E. Detournay. An implicit algorithm for the propagation of a hydraulic fracture with a fluid lag. *Comp. Meth. Appl. Mech. Engrg.*, 196:4863–4880, 2007.

[19] Z.R. Chen. Finite element modelling of viscosity-dominated hydraulic fractures. *Journal*

volume 10 of *Iutam Bookseries*, pages 91–100, Dordrecht, 2009. Springer.

*Society of London. Series A: Mathematical and Physical Sciences*, 447(1929):39, 1994.


**Section 14**

**Mining and Measurement**
