International Journal of Fracture (2005) 134:175–190 DOI 10.1007/s10704-005-0154-0

© Springer 2005

Toughness-dominated hydraulic fracture with leak-off ANDREW P. BUNGER1 , EMMANUEL DETOURNAY2,∗ and DMITRY I. GARAGASH3 1

CSIRO Petroleum, Melbourne, Vic., Australia Civil Engineering, University of Minnesota, 122 Civil Engineering, 500 Pillsbury Drive, Minneapolis, MN, 55455, U.S.A. 3 Clarkson University, Potsdam, NY, U.S.A. *Author for correspondence. (E-mail: [email protected]) 2

Received 20 April 2004; accepted in revised form 15 June 2005 Abstract. This paper considers the problem of a hydraulic fracture in which an incompressible Newtonian fluid is injected at a constant rate to drive a fracture in a permeable, infinite, brittle elastic solid. The two cases of a plane strain and a penny-shaped fracture are considered. The fluid pressure is assumed to be uniform and thus the lag between the fracture front and the fluid is taken to be zero. The validity of these assumptions is shown to depend on a parameter, which has the physical interpretation of a dimensionless fluid viscosity. It is shown that when the dimensionless viscosity is negligibly small, the problem depends only on a single parameter, a dimensionless time. Small and large time asymptotic solutions are derived which correspond to regimes dominated by storage of fluid in the fracture and infiltration of fluid into the rock, respectively. Evolution from the small to the large time asymptotic solution is obtained using a fourth order Runge–Kutta method. Key words: Asymptotic expansions, hydraulic fracturing, numerical methods.

1. Introduction The problem of fluid-driven fracture in a brittle medium is one with a number of applications in engineering and geophysics. For example, hydraulic fracturing has been used to stimulate about 70% of the gas wells and 50% of the oil wells drilled in North America in the last half-century (Economides and Nolte, 2000). Fluid-driven fracture is also an important mechanism by which magma is transported in the lithosphere. The economic and scientific importance of fluid-driven fracturing of rock has motivated many contributions since the pioneering work of Khristianovic and Zheltov (1955). This paper considers the particular problem of a plane strain, usually referred to as a KGD fracture (in honor of Khristianovic Geertsma and de Klerk), and a radially symmetric (penny-shaped) hydraulic fracture, with fluid leak-off under uniform fluid pressure. The problem of the KGD fracture under uniform pressure in a permeable elastic medium has previously been considered by Huang and Swadener (1993), who contributed the framework for the numerical scheme employed in the present study. While it is important to acknowledge this contribution and to note that many other contributions have been made which examine hydraulic fracture growth under some

176 A.P. Bunger et al. geometric and/or parametric idealizations (e.g. Cleary and Wong, 1985; Spence and Sharpe, 1985) much has been gained in the last decade by understanding hydraulic fracture as a problem that evolves between regimes in which the solution is dominated by certain dimensionless groups of parameters (Adachi, 2001; Savitski and Detournay, 2002; Detournay, 2004; Garagash and Detournay, 2005; Garagash, 2005a). For example, over a certain range of time, the behavior of a solution may be shown to depend only on a group of parameters which may be interpreted as a dimensionless viscosity while the material toughness and fluid leak-off have negligible effect. Using this approach, it can be shown that laboratory scale experiments are often dominated by material toughness with hardly any influence from the viscosity, even when the fracture is driven by a very viscous fluid. Formal treatment of the dominance of certain groups of parameters can lead to significant reduction in the complexity of the problem without resorting to ad hoc assumptions or heuristic arguments. This paper outlines the application of this approach to the particular cases described above in order to understand the validity of the uniform pressure assumption, which parameters govern the evolution of the hydraulic fracture, how the solution behaves at small and large time, and how the behavior evolves between the small and large time solution. The problem may be stated as follows. Consider a fracture that is driven by injecting an incompressible Newtonian fluid with viscosity µ at a constant rate Qo . The fracture is driven through a brittle elastic material that is characterized by Young’s modulus E, Poisson’s ratio ν, fracture toughness KI c , and which is subjected to far-field stress σo , perpendicular to the fracture plane. Loss of fluid into the solid material is characterized by the leak-off parameter Cl . The material parameters may alternately be given by E E = , 1 − ν2 



µ = 12µ,


1/2 2 K =4 KI c , π 

C  = 2Cl ,

where E  is the so-called plane strain modulus and the other parameters are introduced to simplify the upcoming equations. Assume that the radius of the wellbore is small relative to the length of the fracture, thus fluid injection may be modeled as a point source at the origin. Furthermore, assume that the lag between the fracture front and the fluid is much smaller than the length of the fracture and may be taken as zero. This assumption can be shown to be justified for the regime of propagation to be considered (Garagash, 2005b). Solving this problem consists of obtaining solutions for the net pressure p(x, t), which equals the fluid pressure pf minus the farfield stress σo , the fracture aperture w(x, t), and the fracture length l(t). A sketch of the problem geometry is shown in Figure 1. 2. Mathematical formulation 2.1. Governing equations Consider a planar (KGD) and radially symmetric (penny-shaped) crack geometries characterized by the fracture half-length or fracture radius, l(t), respectively, and aperture w(x, t) where x is the spatial coordinate. The mathematical model begins with global conservation of incompressible fluid injected at the crack inlet

Toughness-dominated hydraulic fracture with leak-off 177

Figure 1. Geometry of the KGD (plane strain) hydraulic fracture.

Vinject (t) = Vcrack (t) + Vleak (t).

(1)

(Note that the “volume” in the planar geometry is per unit out-of-plane fracture extent and, consequently, has dimension of an area). The injected volume is assumed to be given. We consider the particular case of a constant injection rate, hence Vinject (t) = Qo t. The volume of the fluid stored in the crack and leaked into the porous solid is 

l(t)

Vcrack (t) = 2

w(x, t) (πx)m−1 dx,

t > 0,

(2)

0



t

Vleak (t) = 2

dt 0





l(t  )

g(x, t  ) (πx)m−1 dx,

t > 0,

(3)

0

where m = 1 and m = 2 for planar and radially symmetric geometries, respectively. Here g(x, t) is the volume of fluid leak-off per unit length of the fracture (Carter, 1957) C g(x, t) = √ , t − to (x)

t − to (x) > 0,

(4)

which is found by modeling leak-off as a one-dimensional diffusion process, where to (x) is the time it takes for the fluid front to reach point x. Substituting (4) into (3), reversing the order of integration, and integrating by parts yields the leak-off volume Vleak (t) = 2

 π m−1  2

0

t

C  l m (t  )  dt , √ t − t

t > 0.

(5)

178 A.P. Bunger et al. In addition to these global fluid conservation considerations, fluid in the crack is governed by the lubrication equation
 ∂w 1 ∂ m−1 3 ∂p + g =  m−1 x , t > 0, x < l(t), (6) w ∂t µx ∂x ∂x which is obtained by eliminating the fluid flux q between local mass conservation ∂w 1 ∂(x m−1 q) + m−1 + g = 0, ∂t x ∂x

t > 0, x < l(t),

and Poiseuille law q =−

w3 ∂p , µ ∂x

t > 0, x < l(t).

Further to these conditions on fluid flow, linear elastic fracture mechanics gives 
 l 8 x x w= , p(x  , t)x m−1 dx  , m = 1, 2, Gm (7) πE  l m−1 0 l l where G1 and G2 are known elastic kernels for planar and radial cracks, respectively, (Sneddon, 1951; Sneddon and Lowengrub, 1969; Barr, 1991)     1  1 − ξ 2 + 1 − ξ 2   G1 ξ, ξ = ln  , 2  1 − ξ 2 − 1 − ξ 2 

   −1 1−ξ 2 ξ 2 1  F sin , ,   2 2 ξ

1−ξ ξ 2   G2 ξ, ξ  =  1 F sin−1 1−ξ 2 , ξ , ξ 1−ξ 2 ξ 2

ξ > ξ , ξ <ξ

with F denoting the incomplete elliptic integral of the first kind (Abramowitz and Stegun 1965). Thus the field equations which govern this problem are given by (1), (2), (5), (6), and (7). Ignoring the initial conditions for now, the formulation is completed by addition of the following boundary conditions, which incorporate the assumption that the fracture is always propagating (K = KI c ), w∼

K (l − x)1/2 , E

q = 0,

x = l.

1−

x  1, l

(8) (9)

2.2. Scaling This problem is governed by two competing energy dissipation mechanisms and two competing fluid storage mechanisms (Detournay and Garagash, submitted). The energy dissipation mechanisms are associated with viscous flow (µ) and with creation of surface area in the solid material (KI c ). The two storage mechanisms are storage of fluid in the fracture (E  ) and leak-off of fluid into the permeable solid (Cl ).

Toughness-dominated hydraulic fracture with leak-off 179

Figure 2. Parametric space, after Detournay and Garagash (submitted).

These competing processes are illustrated by the parametric space shown in Figure 2. The base of the rectangle is the storage-dominated regime, the top edge is the leak-off dominated regime, the left edge is the viscosity-dominated regime, and the right edge is the toughness-dominated regime. Under injection of a Newtonian fluid at constant rate, it can be shown that the regime of hydraulic fracture propagation evolves in time from the storage edge to the leak-off edge along a trajectory determined by dimensionless parameter Mk , a monomial grouping power of µ , E  , K  , and Qo , but independent of time. For the KGD fracture, trajectories are parallel to the K K˜ edge, thus the toughness-dominated regime (Mk  1) means simply that this trajectory is very near the K K˜ edge. For the radial fracture, it may be shown that the fracture evolves not only from storage to leak-off, but also from viscosity to toughness (Savitski and Detournay, 2002). Hence, there are two independent timescales, the ratio of which is given by Mk . The toughness-dominated regime corresponds to the case where evolution from viscosity to toughness occurs on a much shorter timescale than evolution from storage to leak-off. Here we consider the case where the trajectory passes sufficiently close to the K vertex such that an “early-time” solution, valid in the neighborhood of K, is an asymptotic solution over some range of time. Comparison of relative sizes of parameters first requires appropriate scaling of the problem. We begin with the general scaling w = εL(ξ, τ ),

p = εE  (ξ, τ ),

l = γ (τ )L,

(10)

where t = t∗ τ,

x = l(t)ξ.

(11)

Appropriate scalings along the toughness edge of the parametric space for the KGD and radial fractures are summarized in the Table 1. Under this scaling the dimensionless formulation of the problem becomes  1 1 τ = γ m (τ ) (ξ, τ ) (πξ )m−1 dξ 2 0  π m−1  τ γ m (τ  ) + dτ  , τ > 0, (12) √  2 τ −τ 0

180 A.P. Bunger et al. Table 1. Small parameter ε, lengthscale L, timescale t∗ and dimensionless parameter Mk for the planar and radial crack geometry. Geometry

ε C 2 Qo

KGD



Radial

C 2 K  2 Qo E  2

1/3

L 

2



2/3

K  Qo E C 2 K  Qo E C 2

t∗ K  4 Q2o E 4 C 6



K  4 Qo E 4 C 5

2/3

Mk   3 µ QKo E 4   4  11 1/3 µ C KE 14 Qo

1 ∂ ξ dγ ∂ − +√ ∂τ γ dτ ∂ξ τ − τo (x)
 ∂ 1 m−1 3 ∂  ξ , τ > 0, |ξ | < 1, = Mk γ m ξ m−1 ∂ξ ∂ξ    8γ 1 Gm ξ, ξ  (ξ  , τ )ξ m−1 dξ  , τ > 0, |ξ | ≤ 1, = π 0  ∼ γ 1/2 (1 − ξ )1/2 , 1 3 ∂  = 0, γ ∂ξ

1 − ξ  1,

(13) (14) (15)

ξ = 1.

(16)

In the case where Mk → 0, solution of (13) subject to the boundary condition (16) requires that = const, thus showing the equivalence between considering the toughness edge solution (Mk = 0) and the assumption of uniform internal pressure distribution. Note that the scaling is such that
 ξ dγ ∂ 1 m ∂ − +√ γ = O(1) ∂τ γ dτ ∂ξ τ − τo (x) for all τ > 0, thus the deviation from constant pressure distribution is of the same order as Mk . Furthermore, the fluid lag cannot exist in the absence of a pressure gradient, thus the zero-lag assumption is valid for small values of Mk (Garagash, 2000, 2005b). Finally, note that for Mk  1, under the assumption that the crack is always propagating, the opening and pressure are given as functions of fracture length by linear elastic fracture mechanics (14)–(15),  = 2−1/2 γ 1/2 (1 − ξ 2 )1/2 ,

= (π/2)m−1 2−5/2 γ −1/2 .

(17)

The global volume balance Equation (12) may be interpreted as an evolution equation governing the position of the fracture front γ as a function of dimensionless time τ . Substituting (17) into (12), evaluating the spatial integrals and changing variables in the temporal integral yields the governing equation  π m−1 √ 1 τ = Am γ m+1/2 + τ 2 2



1 0

γ m (ητ ) dη, 1−η

where A1 = 2−5/2 π and A2 = 3−1 2−1/2 π.

m = 1, 2,

(18)

Toughness-dominated hydraulic fracture with leak-off 181 3. Asymptotic behavior It is useful to examine the small and large time asymptotic behavior, which corresponds to the behavior of the solution near the storage/toughness (K) and leak˜ vertices, respectively, in Figure 2. We assume that for small and off/toughness (K) large values of τ an asymptotic solution exists of the form τ 1:

γ (τ ) = τ β

n 

γk i τ iα + O(τ (n+1)α+β ),

(19)

i=0

τ 1:

γ (τ ) = τ

β˜

n 

˜

˜ β γk˜ i τ i α˜ + O(τ (n+1)α+ ).

(20)

i=0 ˜

β in the asymptotic expansions (19) and (20), The leading terms γk0 τ β and γk0 ˜ τ respectively, correspond to the solution in the storage-dominated (K-vertex) and leak˜ off dominated (K-vertex) regimes, respectively. In the global fluid balance Equation (18), these solutions correspond to balancing of the dimensionless injected fluid volume (left hand side) by the first (storage in the crack) or the second (leak-off) terms on the right hand side, respectively. This yields time-exponents β and β˜

β=

1 , m + 1/2

β˜ =

1 , 2m

and coefficients 1  m+1/2
1 γk0 = , 2Am

(m = 1, 2).

γk0 ˜ =

21−1/m , π

(m = 1, 2).

Next-order terms in (19) and (20) (with i
1) govern the departure of the solution from the corresponding vertex. Time-exponents α and α˜ are determined from the next-order terms in (18) with (19) and (20), respectively, α=

1 m − 1/2 , 2 m + 1/2

α˜ = −

m − 1/2 , 2m

(m = 1, 2).

The rest of coefficients γk i and γk˜ i (with i
1) are determined by solving at each order according to regular perturbation theory (see Appendix B for some details). Results are given in Table 2 and plots are given along with the numerical results in Figures 3–6. The asymptotic solutions (19)–(20) imply the existence of a particular scaling at each vertex in which the dimensionless fracture length tends to a finite limit. The problem is naturally scaled according to the K-scaling at small time and by the ˜ K-scaling for large time (see Appendix A for further detail). These two scalings are given by γ γk = β , Ck = τ α , τ  1, (21) τ γk˜ =

γ τ β˜

,

Sk˜ = τ α˜ ,

τ  1.

(22)

182 A.P. Bunger et al. Table 2. Values required for first five terms of asymptotic series solutions (19) and (20). Formal expressions for the γki ’s and γki˜ ’s are given in Appendix A. τ 1

β α γk0 γk1 γk2 γk3 γk4

τ 1

KGD

Radial

2/3 1/6 0.9324 −1.714 2.196 −1.863 0.7093

2/5 3/10 0.8546 −1.110 1.562 −1.772 1.337

β˜ α˜ γk0 ˜ γk1 ˜ γk2 ˜ γk3 ˜ γk4 ˜

KGD

Radial

1/2 −1/4 0.3183 −5.706 ∗ 10−2 1.341 ∗ 10−2 −3.131 ∗ 10−3 6.368 ∗ 10−4

1/4 −3/8 0.4502 −3.824 ∗ 10−2 4.685 ∗ 10−3 −3.322 ∗ 10−4 −

1

Penny-shaped 0.1

γ

KGD 0.01

0.001 0.001

0.01

0.1

τ

1

Figure 3. Plot of fracture length (solid line) versus time along with small and large time asymptotic solutions. Note that solutions have been shown only over the last 3 decades to make the figure clearer.

Re-scaling in this way allows examination of the limiting small and large time behavior on the time scale (Ck , Sk˜ ) with which the problem is evolving. 4. Numerical solution 4.1. Algorithm It is advantageous to obtain a governing equation in the form of a first order ordinary differential equation, owing to the suite of available tools for solving such problems. This may be done by differentiating (18) to give, γ  (τ ) =

1 (2m + 1)Am γ m−1/2 (τ )



 1 − 2π m−1 0

τ

 γ m−1 (θ)γ  (θ) dθ , √ τ −θ

(23)

Toughness-dominated hydraulic fracture with leak-off 183 10

KGD

Π

1 Penny-shaped

0.1 0.001

0.01

τ

0.1

1

Figure 4. Plot of pressure (solid line) versus time along with small and large time asymptotic solutions.

Figure 5. Plot of fracture length in storage/toughness (K vertex) scaling (solid line) along with small time asymptotic solution (dashed) and limit as Ck → 0 (dotted).

where the integral in (18) is first integrated by parts and then differentiated according to Liebnitz’s rule. Special care must be given to the numerical integration owing to the singular nature of the kernel near 0 and τ . A quadrature scheme based on analytical treatment of the singular parts of the integral and Simpson’s quadrature rule over the regular parts of the integral is derived in Appendix B. In this manner we obtain a governing equation of the form  γi = F (τi , γi , γi−1 , γj
(24)

This history-dependent first order ordinary differential equation is solved using a fourth order Runge–Kutta method, which is a modification of Huang and Swadener (1993). The method employs the recurrence

184 A.P. Bunger et al.

Figure 6. Plot of solution for fracture length in leak-off/toughness (K˜ vertex) scaling (solid line) along with large time asymptotic solution (dashed) and limit as Sk˜ → 0 (dotted).

γi+1 = γi +

h (V1 + 2V2 + 2V3 + V4 ) 6

(25)

with  V1 = F (τi , γi , γi−1 , γj
(26)

V2 = F (τi + 21 h, γi + 21 hV1 , V1 , γj
i , γj
i ),

(27)

V3 = F (τi + 21 h, γi + 21 hV3 , V3 , γj
i , γj
i ),

(28)

V4 = F (τi + h, γi + hV3 , V3 , γj
i , γj
i ).

(29)

Initial conditions are obtained from the small time asymptotic behavior, which is used to obtain γ at the first two time nodes and γ  at the first time node. 4.2. Results Figures 3–4 show the numerical solution, where pressure is obtained from γ via (17). These plots also include the small and large time asymptotic solutions and it is clear that the solution evolves between these asymptotics. Notice that the transient solution matches well with the large-time asymptote even while τ < 1. We owe this to the convergence properties of the convolution integral, which results in each term in the asymptotic series decaying much more rapidly than its respective power of τ . This is easily seen by inspecting the coefficients for the τ  1 series given in Table 2. One must re-scale according to (21) and (22) to properly examine early and late time behavior. The results for fracture length at these two scalings are shown in Figures 5 and 6. When properly scaled the evolution solution approaches the appropriate finite limit at small and large time, that is as the evolution parameters, Ck and Sk˜ respectively, become small.

Toughness-dominated hydraulic fracture with leak-off 185 5. Conclusions The problem of the hydraulic fracture along the toughness edge is equivalent to the problem under the assumption of uniform fluid pressure. The validity of the uniform pressure assumption is properly understood by considering the parameters governing the mathematical model, rather than relying on various heuristic arguments often employed in the literature. It is shown that under the uniform pressure assumption the solution for fracture length, opening, and pressure can be tabulated once for all cases owing to the fact that the solution depends only on one parameter, a dimensionless time. Evolution of the solution between the small and large time asymptotic solutions is obtained numerically. We show that when properly scaled the solution approaches a known finite limit at early and late time. This comprises a strong test of the accuracy of the numerical method. Indeed, little can be said about the accuracy of the numerical results without knowledge of the small and large time asymptotic forms. Physical insight may be gained by examining the asymptotic solutions in dimensional form. The first term in the asymptotic solution for the KGD fracture is given by

l∼

    2 E  Qo 2/3 t 2/3 , 

πK Qo 1/2 t , π C

E 4 C 6 t K 4 Q2o 4 6 E C t K 4 Q2o

 1,  1.

(30)

Similarly for the penny-shaped fracture  2/5   3E√ Qo t 2/5 , π 2K  l∼  1/2   1 2Qo t 1/4 , π C

 

E 4 C 5 K 4 Qo

2/3

 5 2/3

E 4 C K 4 Qo

t  1, t  1.

(31)

Examination of (30)–(31) emphasizes that the solution evolves from a storage-dominated regime associated with E  at small time to a leak-off dominated regime associated with C  at large time. It is easily seen that the asymptotic solutions give very accurate approximation for τ < 0.01 and τ > 0.1. However, even more may be said for the quality of approximation. Figures 3 and 4 show that the small and large time asymptotes give nearly the same result in the region 0.01 < τ < 0.1. In fact, it may be shown that the match is exact for the penny-shaped fracture. Hence, the asymptotic solutions are not only useful for small and large time, but taken together, they yield approximate solutions valid up to a reasonable degree of accuracy for all time. The validity of this model is tied to whether the fracture is propagating such that the geometrical assumptions of plane strain or radial growth are satisfied. Although there may be some field cases for which these geometries provide intermediate asymptotic solutions over some period of time, the value of these solutions does not hinge on whether or not the geometry is realistic. Rather, these solutions provide strong benchmarks for numerical hydraulic fracturing simulators which must be able to produce correct results in these idealized cases if their predictions for more complicated hydraulic fracture growth are to be believed.

186 A.P. Bunger et al. Acknowledgements Support has been provided by a Graduate Assistance in Areas of National Need fellowship from the United States Department of Education and the University of Minnesota Department of Civil Engineering. Further support has been provided by CSIRO Division of Petroleum Resources in Melbourne, Australia. Appendix A. Perturbation Method This appendix details the application of regular perturbation theory to find asymptotic solutions to the problem at hand. This discussion follows the approach described by Lin and Segel (1988). We wish to find a solution to (18) for small values of τ α , where α > 0 for τ  1, or τ α˜ , where α˜ < 0 for τ  1. Letting λ = τ α (τ  1) or λ = τ α˜ (τ  1) and letting F be related directly to γ , we wish to find solutions in the form of the Maclaurin expansion   n  λn+1 dn+1 F (ξ ) i F(n) (λ) = , ai λ + dλn+1 (n + 1)! i=0

1 di F (0) ai ≡ , i! dλi

0 < ξ < λ.

(A.1)

The prerequisite of finding a solution in this way is proper scaling, which means that each term is comprised of a dimensionless factor of unit order and a dimensionless constant which closely estimates the term’s order of magnitude Lin and Segel (1988). Inspection of (18) shows clearly that they are not properly scaled for τ  1 or τ  1. At this point, the proper scaling of γ is unknown, however we assume ˜ that such a scaling exists of the form γk = γ /τ β (τ  1) or γk˜ = γ /τ β (τ  1) and we assume that (A.1) exists for F = γk or F = γk˜ , respectively. Substitution into (A.1) leads directly to (19) with ai = γki (τ  1) or to ( 20) with ai = γki˜ (τ  1). For the KGD case (m = 1) substitution of the small time asymptote (19) into (18) gives  τ π  3/2 1/2 = 5/2 τ 3β/2 γk0 + τ 3β/2+α ( 23 γk1 γk0 ) + O(τ 3β/2+2α ) + 2 2 1 + τ β+1/2 ηβ γk0 + τ β+1/2+α ηβ+α γk1 0  dη +O(τ β+1/2+2α ηβ+2α ) . (A.2) 1−η This equation must be “balanced” by appropriate selection of α and β, that is, the powers of τ must be such that a solution to (A.2) exists and may be found by solving at each order of τ independently, beginning with the term of largest order. Balancing the first terms, namely, the left hand side of (A.2) (the dimensionless injected fluid volume) with the first term in the first square brackets on the right hand side (leading storage term) gives the value of β = 2/3. Similarly, balancing of the next order

Toughness-dominated hydraulic fracture with leak-off 187 terms, namely, the second term in the first square bracket (next-order storage term) with the first term under the integral (first-order leak-off term) gives the appropriate value of α = 1/6. Time-exponents β˜ = 1/2 and α˜ = −1/4 in the large time asymptotic expansion (20) are obtained in analogous fashion from an equation analogous to (A.2), except that in this case the leading term in the expansion is provided by the balancing of the leak-off integral with the injected fluid volume, while the next-order terms are related to the storage in the crack. Solving of (A.2) and of an analogous equation for the large time asymptote term by term gives the values for the γki and γki˜ , respectively. Numerical values of the first six coefficients for the KGD fracture are recorded in Table 2. The corresponding analytical expressions are given below for the first four coefficients (omitting the lengthy expressions for the high-order coefficients for brevity): for τ  1 γk0 =

2 π 2/3

γk1 = − γk2 = − γk3 =

(A.3)

,

16 I2/3 , 3π 4/3 32 I2/3  9π 2

(A.4)  I2/3 − 4I5/6 ,

(A.5)

 512 I2/3  −8I5/6 + I2/3 (2 − I2/3 + 3I5/6 ) , 8/3 81π

(A.6)

1 , π

(A.7)

for τ  1 γk0 ˜ =

1 , γk1 ˜ =− √ 4 2πI1/4 3 , 128 I1/4 √ 3 π(1 + 3I1/4 ) γk3 , √ ˜ =− 2 1024 2I−1/4 I1/4

γk2 ˜ =

(A.8) (A.9) (A.10)

where Iϕ denotes integrals which are solved in terms of the Gamma function (Spiegel and Liu, 1999),  1 ηϕ [ϕ + 1][1/2] Iϕ = . (A.11) dη = [ϕ + 3/2] 1−η 0 This method may be used for the penny-shaped fracture in a completely analogous way. Values of α, β, α, ˜ and β˜ obtained by balancing, are reported in Table 2, along with the numerical values of the expansion coefficients. Interestingly, the asymptotic method adopted here limits the large time asymptotic expansion for the penny-shape crack to four terms (and to six terms for KGD crack). The limitation comes from the convergence properties of the leak-off integral near τ = 0. Numerical values of

188 A.P. Bunger et al. the first six coefficients for small-time asymptote and the first four coefficients for large-time asymptote for the penny-shape fracture are recorded in Table 2 . The corresponding analytical expressions are given below (omitting the lengthy expressions for the high-order coefficients in the τ  1 case for brevity): for τ  1,


2/5 3 γk0 = , √ π 2 √
1/5 3 2I4/5 3 , γk1 = − √ 5 π 2 γk2 =

9I4/5 (8I11/10 − 3I4/5 ), 50

γk3 =

23/5 314/5 π 1/5 I4/5 2 (−2I4/5 − 8I11/10 I7/5 + 125 +I4/5 (6I11/10 + I7/5 )),

and for τ  1, √ 2 γk0 , ˜ = π 21/4 γk1 = − , ˜ 3π 3/2 I1/8

(A.12)

(A.13)

(A.14)

(A.15)

(A.16) (A.17)

γk2 ˜ =−

I−1/4 − 5I1/8 , 2 18π 2 I−1/4 I1/8

(A.18)

γk3 ˜ =−

4I−5/8 (I−1/4 − 5I1/8 ) + 5I1/8 (I−1/4 + 10I1/8 ) . 2 216π 5/2 21/4 I−5/8 I−1/4 I1/8

(A.19)

B. Numerical method This appendix outlines the quadrature scheme used to evaluate the governing integro-differential Equations (18). While this covers specifically the method used for the KGD case (m = 1), the method used for the penny-shaped case is analogous. First, the integral is split into its regular and singular parts. Here we anticipate that the routine will begin at the small initial value of time τo and the time domain is to be discretized uniformly with spacing h. The governing equation may be written as 
 τ −h  τ   τo   γ (θ) 25/2 γ = + + 1−2 dθ . √ 3πγ 1/2 (τ ) τ −θ τo τ −h 0 

(B.B.1)

The kernel is smooth over the domain of the first integral. Letting Fij be a matrix that incorporates the denominator of the kernel and which is constructed according to Simpson’s quadrature rule, we may approximate the first integral as

Toughness-dominated hydraulic fracture with leak-off 189 

τ −h

τo

i−1  γ  (θ) Fij γj . dθ ≈ √ τ −θ j =1

(B.B.2)

The second and third parts of the integral in (B.B.1) must be treated in a special way due to the singular nature of the kernel at θ = 0 (owing to singular γ  ) and at θ = τ . The second part may be integrated analytically by changing variables and using a linear approximation of γ  between τ − h and τ . The result is  τ √  1 (1 − χ )γ  (τ − h) + χ γ  (τ ) γ  (θ) dχ dθ ≈ h √ τ −θ 1−χ τ −h 0 4√  2√  hγ (τ − h) + hγ (τ ). (B.B.3) ≈ 3 3 The third part of the integral is evaluated by approximating γ  by γˆe , where γˆe is the truncated asymptotic series for small time from (19). Furthermore, to enable analytical evaluation of the singular part, we write the integral as
 τ  τ −h  τ    τo  γˆ (θ) γ (θ) − − dθ ≈ dθ. √ √e τ −θ τ −θ 0 τo τ −h 0 When the integral is split in this way, the second part is regular and may be integrated using Simpson’s quadrature rule. The third part may be approximated in the same manner as (B.B.3). Upon an appropriate change of variables, the first part is easily evaluated analytically using the early time asymptote (19) and solutions to integrals of the form (A.11). The discretized form of the governing integro-differential equation may thus be obtained by collecting all the terms of γ  at the current time step on the left hand side. Using the early time asymptote to evaluate the integral near τ = 0 avoids large numerical integration errors that would otherwise arise from the singularity in γ  . However, for this method to work with some prescribed accuracy ε requires that τo be sufficiently small. We make the following order of magnitude approximation. Taking τo = h, we require element size h to be such that h5/6 < ε and h3/2 < ε for the KGD and penny-shaped cases, respectively. This corresponds to the relative size of the next order term when analytical integration is carried out on five terms of the early time asymptotic series. Hence, if 100 elements are used, we would expect ε = 0.01 precision for approximately τ < 0.5 for KGD and τ < 4 for the penny-shaped fracture. The integral in (B.B.1) is a convolution involving γ  (and γ for the penny-shaped case) at all previous time steps. Hence, the number of points to be considered grows each time the routine advances. This growing convolution is a problem that is often encountered when modeling diffusive processes and it can lead to huge computation time if one wishes to compute results over several orders of magnitude in time. This problem was overcome by devising the routine such that the number of points in the convolution is always kept between n/2 and n, where n is some fixed number. Upon reaching n points, the grid spacing is doubled and the process repeated. For example, the results presented for the penny-shaped fracture are computed by keeping between 50 and 100 points in the convolution, using an initial step size of 10−6 , and doubling

190 A.P. Bunger et al. 18 times. In this way the solution was computed for about seven orders of magnitude in just a few seconds on a desktop computer. References Abramowitz, M. and Stegun, I.A. (eds.) (1965). Handbook of Mathematical Functions, Dover. Adachi, J. (2001). Fluid-driven Fracture in Permeable Rock, Ph.D. thesis, University of Minnesota, Minneapolis, MN. Barr, D.T. (1991). Leading-Edge Analysis for Correct Simulation of Interface Separation and Hydraulic Fracturing, Ph.D. thesis, Massachusetts Institute of Technology, Cambridge, MA. Carter, E.D. (1957). Optimum fluid characteristics for fracture extension. In: Drilling and Production Practices. (edited by Howard, G.C. and Fast, C.R.), American Petroleum Institute, Tulsa OK, pp. 261–270. Cleary, M. and Wong, S. (1985). Numerical simulation of unsteady fluid flow and propagation of a circular hydraulic fracture. International Journal of Numerical and Analytical Methods in Geomechanics 9, 1–14. Detournay, E. (2004). Propagation regimes of fluid-driven fractures in impermeable rocks. International Journal of Geomechanics 4(1), 1–11. Economides, M.J. and Nolte. K.G. (2000). Reservoir Stimulation, John Wiley & Sons, Chichester, UK, 3rd Edition. Garagash, D.I. (2000). The tip region of a fluid-driven fracture in an elastic medium. ASME Journal of Applied Mechanics 67, 183–192. Garagash, D.I and Detournay, E. (2000). Hydraulic fracture propagation in elastic rock with large toughness. In: Rock Around the Rim – Proc. 4th North American Rock Mechanics Symposium. (edited by Girard, J., Liebman, M., Breeds, C. and Doe, T.), Balkema, pp. 221–228. Garagash, D.I. and Detournay, E. (2005). Plane strain propagation of a hydraulic fracture: SmallToughness solution. ASME Journal of Applied Mechanics In press. Garagash, D.I. (2005a). Plane-Strain Propagation of a Fluid-Driven Fracture during Injection and Shut-in: Asymptotics of Large Toughness. Engineering Fracture Mechanics. In press. Garagash, D.I. (2005b). Propagation of a Plane-Strain Fluid-Driven Fracture with a Fluid Lag: Early Time Solution. International Journal of Solids and Structure, (Submitted May 2005) Geertsma, J. and de Klerk, F. (1969). A rapid method of predicting width and extent of hydraulic induced fractures. Journal of Petroleum Technology 246, 1571–1581, (SPE 2458). Huang, N.C. and Swadener, J.G. (1993). An investigation of hydraulic fracturing, part I: one-phase flow model. Theorertical and Applied Fracture Mechanics 18, 89–102. Khristianovic, S.A. and Zheltov, Y.P. (1955). Formation of vertical fractures by means of highly viscous fluids. Proceedings 4th World Petroleum Congress, Vol. II, Rome, pp. 579–586. Lin, C.C. and Segel, L.A. (1988). Mathematics Applied to Deterministic Problems in the Natural Sciences, SIAM, Philadelphia. Savitski, A. and Detournay, E. (2002). Propagation of a penny-shaped fluid-driven fracture in an impermeable rock: asymptotic solutions. International Journal of Solids and Structures 39, 6311– 6337. Sneddon, I.N. (1951). Fourier Transforms, McGraw-Hill, New York, NY. Sneddon, I.N. and Lowengrub, M. (1969). Crack Problems in the Classical Theory of Elasticity, John Wiley & Sons, New York, NY. Spence, D. and Sharpe, P. (1985). Self-similar solution for elastodynamic cavity flow. Proceedins of the Royal Society of London: Series A 400, 289–313. Spiegel, M.R. and Liu, J. (1999). Mathematical Handbook of Formulas and Tables, 2nd ed., McGrawHill, New York.