1

Use of a speed equation
for numerical simulation of hydraulic fractures
Alexander M. Linkov
Institute for Problems of Mechanical Engineering, 61, Bol’shoy pr. V. O., Saint Petersburg, 199178, Russia
Presently: Rzeszow University of Technology, ul. Powstancow Warszawy 8, Rzeszow, 35-959, Poland,
e-mail: linkoval@prz.edu.pl

Abstract. This paper treats the propagation of a hydraulically driven crack. We explicitly write the
local speed equation, which facilitates using the theory of propagating interfaces. It is shown that when
neglecting the lag between the liquid front and the crack tip, the lubrication PDE yields that a solution
satisfies the speed equation identically. This implies that for zero or small lag, the boundary value
problem appears ill-posed when solved numerically. We suggest - regularization, which consists in
employing the speed equation together with a prescribed BC on the front to obtain a new BC
formulated at a small distance behind the front rather than on the front itself. It is shown that regularization provides accurate and stable results with reasonable time expense. It is also shown that
the speed equation gives a key to proper choice of unknown functions when solving a hydraulic
fracture problem numerically.
Keywords: hydraulic fracturing, numerical simulation, speed equation, ill-posed problem,
regularization
1. Introduction
Hydraulic fracturing is a technique used extensively to increase the surface to or from which a fluid
flows in a rock mass. It is applied for various engineering purposes such as stimulation of oil and gas
reservoir recovery, increasing heat production of geothermal reservoirs, measurement of in-situ
stresses, control of caving in the roof of coal and ore excavations, enhancing efficiency of CO2
sequestration and isolation of toxic substances in rocks. In natural conditions, a similar process occurs
when a pressurized melted substance fractures impermeable rock leading to the formation of veins of
mineral deposits. Beginning with researchers such as Khristianovich & Zheltov (1955), Carter (1957),
Perkins & Kern (1961), Geertsma & de Klerk (1969), Howard & Fast (1970), Nordgren (1970),
Spence &Sharp (1985), Nolte (1988) numerous studies have been published on the theory and
numerical simulation of hydraulic fracturing (see, e. g., Desroches et al. 1994; Lenoach 1995;
Garagash & Detournay 2000; Adachi & Detournay 2002; Detournay et.al. 2002; Savitski & Detournay
2002; Jamamoto et al. 2004; Pierce & Siebrits 2005; Garagash 2006; Adachi et al. 2007; Mitchel et al.
2007; Kovalyshen & Detournay 2009; Kovalyshen 2010; Hu & Garagash 2010; Garagash et al. 2011
and detailed reviews in many of them). The review by Adachi et al. (2007) is specially organized to
give a comprehensive report on the computer simulation; it actually covers the present state of the art
as well. Thus there is no need to dwell on the historical background, detailed analysis of the processes
in a near-tip zone resulting in particular asymptotics and regimes of flow, general equations and
general approaches used to date. Being interested in computer aided simulation of hydraulic fracturing,
in this paper we would rather focus on a key area that needs to be addressed for further progress in
numerical simulation (Adachi et al. 2007, p. 754); this is the need “to dramatically speed up”
simulators of fracture propagation. Naturally, reasonable accuracy of results should be guaranteed. The
goal cannot be reached without clear understanding of underlying computational difficulties which
strongly influence the accuracy and stability of numerical results and robustness of procedures.
The paper addresses this issue. It presents in detail the results of brief communications by the
author (Linkov 2011a, b). Our prime objective is to delineate and to overcome the computational
difficulty caused, in essence, by strong non-linearity of the lubrication equation and by the moving
boundary. In contrast with the cited previous publications which employed the global form of the mass

2
balance to trace liquid front propagation in time, we explicitly write and use the local speed equation
(SE). The speed function entering the SE may serve to employ methods developed in the theory of
propagating surfaces (Sethian 1999). The SE gives also a key to proper choice of unknown functions,
which are analytical up to the front. We show that at points where the lag between the liquid front and
the crack tip is zero, the lubrication equation yields that its solution identically satisfies the SE. This
implies that for zero or small lag, the problem will appear ill-posed when solved numerically for a
fixed front at a time step, and consequently it requires appropriate regularization to have accurate and
stable numerical results. We suggest a method of regularization that employs the very source of the
difficulty to overcome it. The SE and a BC at the front are used together to derive a new BC
formulated at a small distance
behind the front rather than on the front itself. This leads to regularization which provides accurate and stable numerical results. The Nordgren model serves to
clearly display the computational features discussed. It is also used to obtain benchmarks with five
correct significant digits, at least.
2. Global mass balance, speed equation and ill-posed problem for hydraulically driven fracture
Reynolds transport theorem (e. g. Crowe et al. 2009), applied to the mass of an arbitrary volume of a
medium in a narrow channel between closely located boundaries, after averaging over the channel
width, reads:
dM
w
w
dS
wvn dL ,
(2.1)
dt
t
t
Sm
L
where dM / dt is the external mass coming into or out of the considered volume per unit time, S m is
the middle surface of the volume, w is the width (opening) of the channel, L is its contour, ρ is the
mass density averaged over the width, vn is the normal to L component of the particle velocity also
averaged over the width; the normal n is assumed to be in the plane tangent to S at a considered point.
Applying (2.1) to the total volume of the medium with the middle surface S t and contour Lt at a
moment t we have:
dM e
w
dS
w vn dL .
(2.2)
dt
t
t
St
Lt
Here, and henceforth, the star denotes that a value refers to the medium front. When writing (2.2)
we assume that there is no significant sucking or vaporization through the front. Then the speed V* of
the front propagation coincides with the normal component vn of the particle velocity. Thus we have
the key equation:
dxn *
,
(2.3)
V* vn
dt
where xn* is the normal to the front component of a position vector of a point on the front; in global
coordinates, x x * (t ) is a parametric equation of the front contour Lt . As a rule, the tangent
component of the particle velocity is small as compared with the normal component at the front. In this
case, the front moves with the speed exactly equal to the velocity of fluid particles comprising it, and
we have V* vn* v * , where v * is the vector of the particle velocity averaged over the front width.
For incompressible homogeneous liquid, ρ = const, and (2.2) may be written as the equation of the
liquid volume balance:
dVe
w
dS
w vn dL ,
(2.4)
dt St t
Lt
where Ve

Me / .

3
Comment 1. For 1-D case, the volume conservation equation (2.4) yields

x (t )
*

w( x, t )dx Ve (t ) .

0

Then for a „rigid-wall‟ channel ( w w(x) , w / t 0 ), integration gives the front location as a
function of time x* (t ) f (Ve (t )) . It is easily seen from the last expression that if the width w(x)
decreases fast enough with growing x, then for a prescribed influx Ve (t ) , the front coordinate turns to
infinity at a finite time t* ( x* (t* )
). A solution does not exist for t t* . For instance, in the case of
the constant influx rate q0 ( Ve (t ) q0t ) and exponentially decreasing width ( w a exp( x) , a, α >
0), the front turns to infinity as t
t* a /( q0 ) . There is no solution for t t* . This clearly indicates
that problems involving flow of incompressible liquid in a thin channel are quite tricky. A solution
may not exist or it might be difficult to find the solution numerically, especially when the rigidity of
channel walls is very high (according to Pierce & Siebrits (2005), high rigidity leads to a stiff system
of ordinary differential equations, when solving a boundary value problem by finite differences).
By definition of the flux through the channel width (e.g. Batchelor 1967), we have q wv . Then
in (2.4) w vn
qn is the normal to L component of the flux through the channel width at the liquid
front, and equation (2.4) may be written as
q n* (x * )
.
(2.5)
V* v n (x * )
w* (x * )
This is the speed equation (SE). Its right hand side (r. h. s) defines the so-called speed function.
Emphasize that the SE is general. It is not influenced by viscous properties of a particular
incompressible liquid and by the presence of leak-off or influxes through the walls of a thin channel.
In some cases, for instance, when a liquid flows in a fracture without a lag between the liquid front
and the crack tip, both the flux and opening turn to zero at the liquid front. Then (2.5) takes the limit
form:
q ( x)
.
(2.6)
vn (x* ) lim n
x x w(x)
It is highly apparent that even when w* ( x * ) = 0 and qn* (x* ) = 0, the limit on the r. h. s. of (2.6)
should be finite to exclude the front propagation with infinite velocity. This suggests using the particle
velocity as an unknown in numerical calculations, because it is non-singular in entire flow region
including the front. What is also beneficial, the velocity is non-zero except for flows with return
points.
In formulations of problems for hydraulic fracturing, the average particle velocity
v(x) q(x) / w(x) does not enter equations. Rather, the total flux q through a cross-section is used.
Specifically, the divergence theorem applied to (2.1) in the case of incompressible liquid yields the
continuity equation in terms of the flux q and opening:
w
div q qe ,
(2.7)
t
where qe is the prescribed intensity of distributed external sources of liquid, q is the flux vector,
defined in the tangent plane to St at a considered point. The Poiseuille law, used for a flow of
incompressible liquid in a narrow channel, connects the flux q with the pressure gradient
(2.8)
q
D(w, p)gradp .
Herein, D is a prescribed function or operator, gradient is also confined to the coordinates in the
tangent plane. Substitution of (2.8) into (2.7) gives the lubrication partial differential equation (PDE) at
points of liquid:
w
(2.9)
div D( w, p)gradp qe 0 .
t

Some initial (normally zero-opening) condition is assumed to account for the presence of the time
derivative in (2.9). Being of the second order and elliptic in spatial derivatives, equation (2.9) requires

4
a boundary condition (BC) on the liquid contour Ll. Normally it is the condition of the prescribed flux
q0 at a part Lq and of the prescribed pressure p0 at the remaining part Lp of the contour Ll:
(2.10)
qn (x* ) q0 (x* ) x* Lq ;
p(x* ) p0 (x* ) x* L p .
We see that neither partial differential equations (PDE) (2.7)-(2.9), nor BC (2.10) involve the
particle velocity. The latter enters only the SE (2.5) or its limit form (2.6) for points on the front. At
these points, it is defined by formula (2.3). Since according to (2.8) qn
D(w, p) p / n for any
direction n in a tangent plane, the flux entering the numerator on the r. h. s. of (2.5) is
p
qn
D( w, p )
.
(2.11)
nx x
*

Then for the hydraulic fracturing, the SE (2.5) is specified as:
1
p
vn
D( w, p )
,
w* (x * )
nx x

(2.12)

*

where n is the outward normal to Lt in the tangent plane at a considered point of the liquid front. The
r. h. s. of (2.12) specifies the speed function, which is the basic concept of the theory of propagating
interfaces (Sethian 1999). It may serve to employ level set methods and fast marching methods.
In hydraulic fracture problems, the opening is unknown. To have a complete system, the
lubrication PDE (2.9) is complemented with an equation of solid mechanics connecting the opening w
and pressure p:
(2.13)
A(w, p) 0 .
As a rule, the operator A in (2.13) is prescribed by using the theory of linear elasticity. In addition,
to let the fracture propagate, we need a fracture criterion (Otherwise, the liquid front reaches the crack
contour and stops.). Commonly, at points of the crack contour, authors impose the condition of linear
fracture mechanics:
(2.14)
K I K Ic ,
where K I is the stress intensity factor (SIF), K Ic its critical value.
At the crack contour Lc , the opening is set zero:
(2.15)
w(xc ) 0 .
In general, for physical reasons, excluding negative pressure, which tends to minus infinity at
points of the front, where the latter coincides with the crack contour, the liquid surface St is solely a
part of the crack surface Sc so that the liquid contour Lt is within the crack contour Lc . Thus, in
general there is a lag between Lt and Lc . In this case, the second of the conditions (2.10) is prescribed
at the boundary of the propagating liquid front:
(2.16)
p(x* ) p0 (x* ) .
In many cases, the change of the pressure p0 (x* ) along the liquid front is small. Consequently,
the tangential derivative p / of the pressure is small as compared with the normal derivative
p / n , and (2.8) implies that the tangential component of the flux is small as compared with its
normal component. Since v q / w , it means that the tangential component of the fluid velocity at the
front may be neglected. Then, as mentioned above, the front velocity equals to the particle velocity
itself. This case will be under discussion further on. We see that the existence of a lag, whatever small
it is, has important physical implications for fracture propagation.
Normally the lag is small and accounting for it strongly complicates a problem, while neglecting it
may be justified (e. g., Garagash & Detournay 2000). For these reasons, many papers on hydraulic
fracturing (e. g., Spence & Sharp 1985; Lenoach 1995; Garagash & Detournay 2000; Adachi &
Detournay 2002; Detournay et. al. 2002; Savitski & Detournay 2002; Jamamoto et al. 2004; Pierce &
Siebrits 2005; Adachi et al. 2007; Mitchel et al. 2007; Hu & Garagash 2010; Garagash et al. 2011)
assume that the lag is zero. Then at all points of the propagating liquid front, coinciding in this case

5
with the crack contour, the flux is zero: qn* (x* ) 0 . The latter condition is met in view of (2.15),
because the operator D(w, p) in (2.11) is such that D(0, p) 0 . Still, to satisfy both (2.15) and (2.14),
the elasticity equations require specific asymptotic behavior of the opening with the coefficient of the
asymptotic proportional to the SIF, when the latter is not zero. The SIF depends on the pressure (see, e.
g. Spence & Sharp 1985). As a result, at points of the front we have the boundary condition (2.14)
which, similar to (2.16), involves the pressure.
The fact that the lag is small and it is commonly neglected has significant consequences for
numerical calculations. To show it, consider the problem in the local Cartesian coordinates x1 ' O1 ' x2 '
with the origin O1 ' at a point x* at the liquid front and the axis x1 ' opposite to the direction of
external normal to the front at this point. We employ aforementioned advantages of using the particle
velocity v q / w . In the local system, for points close to the front, the normal component of the
1
p
velocity is v
. Then the lubrication PDE (2.9) at points near the front takes the form:
D( w, p)
w
x1 '
v
ln w
ln w
(2.17)
(v v* )
qe 0 ,
x1 '
x1 '
t x ' const
1

where the partial time derivative is evaluated under constant x1 ' . Using ln w serves to account for an
0 ) at the front. In
arbitrary power asymptotic behavior of the opening w( x1 ' , t ) C (t )(x1 ' ) (
particular, in the case of Newtonian liquid, for zero lag, the exponent α = 2/3 when in (2.14) K Ic 0
and α = ½ when K Ic 0 (Spence & Sharp 1985). In the case of non-Newtonian liquid, formulae for α
are given by Adachi & Detournay (2002). For the Nordgren problem discussed below, α = 1/3. Note
that for a flow with Carter‟s leak-off, so-called intermediate asymptotics may appear (e.g., Lenoach
1995; Kovalyshen & Detournay 2010). These asymptotics manifest themselves at some distance from
a crack contour rather than at the contour itself. For this reason, we shall not use them below in
equations involving a point on the front.
When the opening has the power asymptotic near the front, it is reasonable, in addition to the

particle velocity, to use also the variable y w1 / , which is linear in x1 ' near the front. The PDE
(2.17) then becomes
v v* y
v
1 y
(2.18)
qe 0 .
x1 '
y
x1 '
y t x ' const
1
Note that the derivative v / x1 ' and, under the supposed asymptotic, the derivatives y / x1 ' , the
multiplier (v v * ) / y and the term (1 / y) y / t in (2.18) are finite at the liquid front. Therefore, the
D( w, p )
grad p and y w1 / present a proper choice of unknown functions for
w
problems of hydraulic fracturing.
Write (2.18) as
v
y
y
(2.19)
y
(v v* )
yqe 0 .
x1 '
x1 '
t x1' const

variables v

Equation (2.19) implies that v v* at any point, where the opening w, and consequently y w1 / ,
is zero. Such are points at the front for zero lag. This means that when neglecting the lag, the SE (2.12)
is satisfied identically by a solution of the PDE (2.18). Obviously, the same holds for a solution of the
starting PDE (2.9).
We see that for zero lag, when solving the boundary value (BV) problem for (2.9), one implicitly
has satisfied the SE (2.12) additional to the prescribed BC of zero opening (2.15). Note that for the
mentioned power asymptotic of the opening, the SE may be re-written in terms of the normal

6
derivative of y

w1/

as

y
x1 ' x

Bv * , where B is a function of time only. Then at a point of the

1 0

front, we have satisfied two equations: y ( x* )

0 and

y
n

Bv* . Recall that the operator
x x*

div D( w, p)grad p is elliptic and requires only one BC, whereas actually there are two BC.
Consequently, we have a Cauchy problem for the elliptic operator. As known (e. g. Lavrent‟ev &
Savel‟ev 1999), such a problem is ill-posed in the Hadamard sense (1902). To have accurate and stable
results when solving an ill-posed problem numerically, one needs its proper regularization (e. g.
Tychonoff 1963; Lavrent‟ev & Savel‟ev 1999).
In the case when the lag is not neglected, we come to similar conclusions if the lag is small
enough. In this case, at a point of the liquid front, we have the BC of (2.16) type. As the lag is small,
the opening w* (x* ) at the liquid front is small, as well. Then from (2.19) it follows that for a solution
of the lubrication PDF, the SE is met approximately. Therefore, at points of the liquid front, in addition
1
p
D( w, p )
vn
to the BC (2.16) for the pressure, we have the approximate equation
w* (x * )
nx x
*

for its normal derivative. Hence, the problem of solving the lubrication PDF will appear ill-posed in
numerical calculations having the accuracy less than that of the approximate equation. It is reasonable
to have a method of regularization, which removes computational difficulties caused by the discussed
feature.
Comment 2. If a problem is self-similar and the lag is neglected, then integration of the lubrication
equation in the automodel coordinate (e. g. Spence & Sharp 1985; Adachi & Detournay 2002) from the
liquid front removes the difficulty: it is sufficient to seek the solution by taking into account for
asymptotic representation of the opening and pressure near the liquid front. In papers by Spence &
Sharp (1985) and Adachi & Detournay (2002), such representations meet both conditions (2.15) and
(2.6), which uniquely define the coefficients of the asymptotics at one point (the crack tip). In fact, the
authors solve an initial value (Cauchy) well-posed problem. The boundary condition of the prescribed
flux at another point (the inlet) is not used: the corresponding influx is found after obtaining the
solution of the Cauchy problem. Similar approach is applicable to the Nordgren problem. It serves us
to obtain the benchmarks in Sec. 4. Unfortunately, this method cannot be applied in a general case
when a self-similar formulation is not available.
3. Nordgren problem. Evidence of ill-posed problem
Consider the Nordgren (1972) model to see unambiguously that the BV problem is ill posed, to find a
proper means for its regularization and to obtain accurate numerical results, which may serve as
benchmarks. The analysis also confirms that the front velocity does satisfy the asymptotic equation
(2.6) despite that both the opening w(x ) and the flux qn (x ) are zero at the front.
3.1. Problem formulation
Recall the assumptions of the Nordgren (1972) problem. Similar to the Perkins-Kern (1961) model, it
is assumed that a vertical fracture of a height h (Fig. 1) is in plane-strain conditions in vertical cross
sections perpendicular to the fracture plane. The cross section is elliptical and the maximal opening w
decreases along the fracture. Nordgren‟s improvement of the model includes finding the fracture
length x (t ) as a part of the solution. Nordgren also accounts for the fluid loss due to leak-off. The
corresponding term is actually a prescribed function of time; we neglect it to not overload the analysis.
In this case, the continuity equation (2.7) reads q / x w / t 0 , where w is the average opening in
a vertical cross section, q is the flux through a cross section divided by the prescribed height h. The

7
liquid is assumed Newtonian with the dynamic viscosity
kl

. Then in (2.8), D( w, p)

kl w3 , where

1 /( 2 ) in the case of an elliptic cross section, considered by Nordgren (1972); for an arbitrary

thin plane channel, the Poiseuille value kl 1 /(12 ) is often used (e. g., Garagash & Detournay 2000;
Savitski & Detournay 2002; Garagash 2006; Mitchel et al. 2007; Hu & Garagash 2010). Thus the
equation (2.8) becomes:
p
q
kl w3
.
(3.1)
x
The dependence (2.13) between the average opening and pressure is taken in the simplest form
p kr w ,
(3.2)
found from the solution of a plane strain elasticity problem for a crack of the height h;
is the Poisson‟s ratio. Therefore, for a
k r (2 / h) E /(1 2 ) , E is the rock elasticity modulus,
non-negative opening, the pressure is non-negative behind the front. With lag neglected, the condition
(2.15) in view of (3.2) implies that the pressure becomes zero at the liquid front. The opening (as well
as the pressure) should then be positive behind the front and zero ahead of it. Under these assumptions,
there is no need in the fracture criterion (2.14).
In view of (3.1) and (3.2), the continuity equation (2.9) becomes the Nordgren PDE:
2 4
1
w
w
kl k r
0.
(3.3)
2
4
t
x
It is solved under the initial condition of zero opening
(3.4)
w(x,0) 0
for any x along the prospect path of the fracture.
The BC for the partial differential equation (PDF) (3.3) includes the condition of the prescribed
influx q0 at the fracture inlet x = 0:
(3.5)
q(0, t ) q0
and the condition that there is no lag between the crack tip and the liquid front:
w( x , t ) 0 .
(3.6)
The solution should be such that the opening is positive behind the front and zero ahead of it:
w( x, t ) 0 0 x x , w( x, t ) 0 x x .
(3.7)
The Nordgren problem consists in finding the solution of PDE (3.3) under the zero-opening initial
condition (3.4) and the BC (3.5), (3.6). The solution should comply with (3.7).
3.2. Speed equation, self-similar problem formulation, clear evidence that the problem is ill-posed
Nordgren used the conditions (3.7) to find the front propagation, rather than the global mass balance
commonly employed for this purpose (e. g., Howard & Fast 1970; Spence & Sharp 1985; Adachi &
Detournay 2002; Savitski & Detournay 2002; Jamamoto et al. 2004; Garagash 2006; Adachi et al.
2007; Mitchel et al. 2007; Hu & Garagash 2010; Garagash et al. 2011). For our purposes, we employ
the SE (2.12). In the case considered, it becomes:
1
w3
v
kl k r
,
(3.8)
3
x x x
where according to (2.3) v dx / dt .
Introduce dimensionless variables:
x
x
t
xd
, xd
, td
,v d
xn
xn
tn

dx d
dtd

v
, wd
vn

w
, pd
wn

p
, qd
pn

q
, q0d
qn

q0
qn

where xn (kl k r )1 / 5 qn3 / 5t n3 / 5 , v n xn / tn , wn qnt n / xn , pn 4k r wn , and t n , qn are arbitrary
scales of the time and flux respectively. In terms of dimensionless values, the equations (3.1), (3.2)

8

pd
, pd
xd
become, respectively,
read qd

wd 3

wd4
. Then the PDE (3.3) and the BC (3.5), (3.6)
xd

4wd ; hence, qd
2

w4
x2
w4
x

w
t

0,

(3.9)

q0 ,

(3.10)

x 0

w( x , t ) 0 .
(3.11)
From this point on, we omit the subscript d at variables and consider only dimensionless values. The
homogeneous conditions (3.4) and (3.7) do not change their form. Nordgren (1972) solved the problem
by finite differences not using the SE.
In dimensionless variables, the SE (3.8) takes the form:
4 w3
v
.
(3.12)
3 x x x

In compliance with the said in Sec. 2, the SE (3.12) gives a key to proper choice of the unknown
function. Indeed, to have the front velocity finite, the partial derivative w3 / x should not be singular
at the liquid front. This yields that w 3 is an analytical function of x and it may be represented by a
x at any instant t: w3 ( x, t )

power series in x*

j 0

a j (t )( x*

x) j . The zero-opening condition (3.11)

gives a0 (t ) 0 , while the SE (3.12) gives a1 (t ) 0.75v* (t ) . Then at the vicinity of the crack tip we
have
w( x, t )

w3 ( x, t )
a11/ 3 ( x*

a1 (t )( x*

x) O ( x*

x)1/ 3 O ( x*

x) 2

what

means

that

the

opening

behaves

as

x) 2 / 3 . Therefore, near the crack tip, in contrast with w3 , the

derivative w / x of w is singular as ( x*

x)

2/3

.

We see that using w 3 avoids unfavorable asymptotic behavior of w, while the SE (3.12) governs
the linear asymptotic behavior of w 3 . Hence, it is reasonable to use w3 as the unknown function,
rather than w as used by Nordgren (1972) or w 4 : in contrast with the derivatives of w3 , the
2 4
2 3
w
w
4 w w3
w
w
derivatives
and
are singular at the front x x . In terms of w3
2
x
3 x x
x
x2
PDE (3.9) becomes:
2

2

w3
1
w3
1 w3
0.
(3.13)
x
x 2 3w3
4 w3 t
For further discussion we use the fact that the problem is self-similar which serves to reduce PDF
(3.13) to an ordinary differential equation (ODE). We express the variables x and w via automodel

variables

and

function y( )

as x
3

t 4 / 5 , w( x)

t1 / 5 ( xt 4 / 5 ) . Then (3.13) becomes ODE with the unknown

( ):
d2y

where a( y, dy / d , )

(dy / d

dy

3
20

0,
d
0.6 ) /(3 y) . The BC (3.10) and (3.11) read:
2

a( y, dy / d , )

dy

0.75
0

y( )

q0
,
3 y (0)

0,

(3.14)

(3.15)
(3.16)

9
and the SE (3.12) becomes:

dy

0.6 .

(3.17)

Re-write (3.14) by using the expression for a( y, dy / d , ) as

d2y

1 dy
0.6
3

dy

3
y
20

(3.18)
0,
d 2
In limit
, for a solution, satisfying the BC (3.16), ODE (3.18) turns into the SE (3.17).
Hence, for ODE (3.14), at the point
, we have imposed not only the BC (3.16) for unknown
function y , but also the BC (3.17) for its derivative dy/ d . Note that equations (3.16), (3.17) imply
that the factor a in (3.14) is finite at the liquid front: lim a( y, dy / d , )
1 /(3 ) .

y

It is easy to check by direct substitution that if y1 ( 1 ) is the solution of the problem (3.14)-(3.17)
for q0 q01 with
y1 ( 2 k ) / k is the solution of the problem (3.14)-(3.17)
1 , then y 2 ( 2 )
k 5 / 6q01 with

for q02
C

( q0 )

0.6

/

and C0

2

y ( 0) /

1/
2

k ; herein, k is an arbitrary positive number. This implies that

are constants not dependent on the prescribed influx q 0 . As
(q0 ) / C , it is a matter of convenience to prescribe q 0 or . A particular value of q 0 or
may be also taken as convenient. Indeed, with the solution y1 ( 1 ) for q0 q01 , we find the solution
0.6

for any q0 : y( ) y1 ( k ) / k , where k (q01 / q0 ) 6 / 5 ,
1 / k ).
1/ k (
Let us fix . According to (3.16), (3.17), at the point , we have prescribed both the function y
and its derivative dy / d . Thus, for the ODE of the second order (3.14) we have a Cauchy (initial
value) problem. Naturally, its solution defines y(0) and dy / d
0 and consequently the flux q 0 at

0 . Hence, even a small error when prescribing q 0 in (3.15), excludes the existence of the solution
of the BV problem (3.14)-(3.16). Therefore, by Hadamard (1902) definition (see also Tychonoff 1963;
Lavrent‟ev & Savel‟ev 1999), the BV problem (3.14)-(3.16) is ill-posed. It cannot be solved without a
proper regularization (Tychonoff 1963; Lavrent‟ev & Savel‟ev 1999). To make conclusions on the
accuracy of numerical results obtained without and with regularization, it is reasonable to obtain
benchmarks.
4. Benchmark solution
The initial value (Cauchy) problem (3.14), (3.16), (3.17) is well-posed. Thus its solution provides the
needed benchmarks. To solve the system we transformed the problem (3.14), (3.16), (3.17) to the
equivalent problem in two unknowns Y1 ( ) y ( ) and Y2 ( ) dy / d . This yields the system of two
ODE:
dY1
Y2
d
(4.1)
dY2
3
a(Y1 , Y2 , )Y2
d
20
under the Cauchy conditions at the point , corresponding to (3.16) and (3.17),
Y1 ( ) 0 , Y2 ( )
0.6 .
(4.2)
The Cauchy problem (4.1), (4.2) is solved by using the fourth order Runge-Kutta scheme (see, e.g.
Epperson 2002). Calculations were preformed with double precision. For certainty, we set
= 1= 1

10
(calculation with other values of
gave the same results for C and C0 to the seventh significant
digit including). The integration step
was changed from 10-2 to 10-5, and the number of steps was
consequently changed from 102 to 105 to reach the inlet point
0 . The number of iterations for the
nonlinear factor a( y, dy / d , ) was changed from 20 to 1000. We could see that the step 10-5 and the
number of iterations 50 are sufficient to guarantee at least six correct digits. The values of the
constants C and C0 are evaluated to the accuracy of seven digits:
C 0.7570913 , C0 0.5820636.
For values depending on , we shall use the subscript 1 when they correspond to
= 1 = 1.
Thus we have:
q01 C*5 / 3 0.6288984 , 1 (0) 3 C0 =0.8349418.
Values for an arbitrary flux q0 may be obtained as

q0

0.6

/C

1.3208446q0

0.6

( 0)

,

3

C0

2

1.0051356 q00.4 .

(4.3)

1.0073486 ,
For the value q0 2 / , used by Nordgren (1972), equations (4.3) give
(0) 0.8390285 against the values given by this author to the accuracy of about one percent:
1.01 , (0) 0.83. The values of 1 3 y1 and d 13 / d 1 are presented in Table 1 with five
correct digits.
Values of ( ) and d 3 / d for an arbitrary prescribed flux q0 may be obtained from those in

the

Table

1

3

as

y1 (

k)/k ,

d

3

/d

3
1

(d

/ d 1) / k

k

1

with

k

(q01 / q0 )1.2

0.5731872 / q10.2 .

Comment 3. Table 1 shows that the derivative d

3

/ d , defining the particle velocity, is nearly

constant along the entire liquid being close to its limit value at the front
3

3

approximately linear in . Next, from the BC (3.16) we obtain
( ) / (0) 1
approximate analytical solution of the Nordgren problem is given by the equation:
w( x, t ) / w(0, t ) (1 x / x )1/ 3
0.2

3

. This implies that
/

is

. Hence, the
(4.4)

0.8

with w(0, t ) t
are found
(0) and x (t )
t ; for a given flux q0 , the values of (0) and
from (4.3). As clear from Table 1, the error of (4.4) does not exceed one percent. The graph
corresponding to the approximate solution (4.4) is indistinguishable from that given by Nordgren
(1972). Naturally, the asymptotic of the solution (4.4) agrees with the predicted asymptotic behavior
a11 / 3 ( x*

w( x, t )

, and consequently, w 3 is proportional to x

proportional to
The
(d

3
1

accurate

/ d 1) /

k

1

0.6
1.3208446q0 ,

data
k

with k

x.

(0) 1.0051356 q00.4 ,

( )

,

( ) and d

3

/d

0.2

, w( x, t )

d

3

/d

are known, we can find the front location

x (t ) , the front velocity v (t ) , the opening w( x, t ) and the particle velocity

0.8 t

and

0.5731872 / q10.2 serve to estimate errors when solving the Nordgren

problem as a BV problem. When

v (t )

3 is

x)1 / 3 near the crack tip because the conditions (4.2) guarantee that

t 0.2 ( xt

0.8

) and ( x, t )

4 0.2 d
t
3
d

( x, t ) as x (t )

3
xt 0.8

.

5. Straightforward solving self-similar BV problem. Method of regularization
5.1. Straightforward integration by finite differences

t 0.8 ,

11
Forget for a while about the SE and all the said on its influence on a BV problem. Let us see what
happens when solving the BV problem (3.14)-(3.16) in a common way by finite differences.
We performed hundreds of numerical experiments with various numbers of nodal points and
iterations and different values of the prescribed influx q01 at the inlet. Finite difference approximations
of second order for d 2 y / d 2 and dy/ d were combined with iterations for a( y, dy / d , ) . Up to
100 000 nodal points and up to 1500 iterations were used in attempts to reach the accuracy of three
correct digits, at least. The attempts failed: by no means could we have more than two correct digits.
Moreover, the results always strongly deteriorate near the liquid front. The numerical results clearly
demonstrate that the BV problem (3.14)-(3.16) is ill-posed. It cannot be solved accurately without
regularization.
As illustration, the dashed line in Fig. 2 presents a typical graph of d 3 / d , obtained under the
Nordgren boundary value q01 2 / . For comparison, the benchmark values of d 3 / d , calculated
by using the Table 1, are shown by the solid line with markers. Obviously, the results strongly
deteriorate near the liquid front
(in the considered example, the benchmark value of
equals
1.0073486).
Comment 4. Using the variable 3 ( ) , which is linear near the liquid front, removes a suggestion
that the deterioration is caused by singularity of d / d and d 2 / d 2 at the point
.
Comment 5. It is worth noting that the accuracy of two correct digits was obtained at points not too
close to the front even when using a rough mesh with a hundred or even only ten nodes. This indicates
that using a rough mesh may serve to regularize a problem when high accuracy is not needed.
5.2.

- regularization.

The numerical experiments evidently confirm that the considered ill-posed BV problem (3.14)(3.16) cannot be solved accurately without regularization. A regularization method is suggested by the
conditions (3.16), (3.17). Indeed, we may use them together to get the approximate equation
y 0.6 (
) near the front. Hence, instead of prescribing the BC (3.16) at the liquid front
,
where it is implicitly complimented by the SE (3.17), we may impose the boundary condition, which
combines (3.16) and (3.17) at a point
(1 ) at a small relative distance
1
/ from the
front:
0 .6 2 .
(5.1)
The BV problem (3.14), (3.15), (5.1) is well-posed and may be solved by finite differences.
Numerical implementation of this approach shows that with
= 10-3, 10-4 the results for the step
/ = 10-3, 10-4, 10-5, 10-6 coincide with those of the benchmark solution. The time expense is
fractions of a second. The results are stable if and
are not simultaneously too small (both and
-5
are greater than 10 ).
As could be expected, the results deteriorate when both the regularization parameter and the step
become too small. Specifically, when
=
= 10-6, the results are completely wrong. Actually,
in this case, to the accuracy of computer arithmetic, the problem is solved without regularization.
We could also see that with growing step
, the accuracy decreases and for a coarse mesh it
actually does not depend on the regularization parameter. In particular, for a quite coarse mesh with
the step
= 0.1, the accuracy is about one percent, and the results stay the same to this accuracy for
any
from 10-2 to 10-9.
The essence of the suggested regularization consists in using the SE together with a prescribed BC
to formulate a BC at a small distance behind the liquid front rather than on the front itself. We call
such an approach -regularization. The next section contains its extension to the cases when a selfsimilar formulation is not available or is not used.
y(

)

12
6. Straightforward solution of starting BV problem. Regularization
6.1. Straightforward integration by time steps with finite differences on a time step
Forget again about the SE and its influence on numerical solution of a BV problem. Try to solve the
starting Nordgren problem by common finite differences. Nordgren (1972) used straightforward
numerical integration of the problem (3.9)-(3.11) under the zero-opening initial condition with the
conditions that opening is positive behind the liquid front and zero ahead of it. This author applied
Crank-Nicolson finite difference scheme to approximate PDF (3.9) and to meet the BC (3.10), (3.11).
The resulting non-linear tridiagonal system was linearized by employing linear approximation of w 4 .
Nordgren (1972) does not include details of calculations on the initialization, the time step, the number
of nodes in spatial discretization, the number of iterations, stability of numerical results and expected
accuracy. To obtain knowledge on these issues, we also solved the system (3.9)-(3.11) in a
straightforward way by using the Crank-Nicolson scheme. The results are as follows.
Actually performing 20 iterations to account for the non-linear term w 4 is sufficient to reproduce
four digits of the fracture opening, except for close vicinity of the liquid front (Increasing the number
to 100 iterations does not improve the solution for all tested time and spatial steps.). For various time
-2
-3
-4
-2
-3
-4
steps ( t
10 , 10 , 10 ) and different spatial steps ( x 10 , 10 , 10 ) taken in various
combinations, the results are stable along the main part of the interval [0, x (t ) ], but deteriorate and
are unreliable in close vicinity of the front ( 1 x / x (t ) < 0.001). This yields changes in the third digit
of
and (0) , calculated by using x (t ) and w(0, t ) . The asymptotic behavior (4.4), as x x (t ) ,
is reproduced near the front except for its close vicinity. The results coincide with those given by
Nordgren (1972) to the accuracy of two significant digits accepted in his work.
In all the calculations, by no means could we have a correct third digit. Similar to self-similar
solution, fine meshes did not improve the accuracy as compared with a rough mesh having the step
x / x = 0.01. The results clearly show that the problem cannot be solved accurately without
regularization.
Comment 6. As mentioned in Sec. 3, the variable w has a singular partial derivative w / x at the
liquid front. To remove the influence of the singularity, we also solved the problem by using w3 as an
unknown function, because according to the SE (3.12) its spatial derivative is not singular. The
conclusions when using w3 are the same as those above. Again, by no means could we have reliable a
third digit, and results strongly deteriorated at a close vicinity of the liquid front. This shows that the
inaccuracy is caused not by singularity of the derivative w / x at the liquid front. It is caused by the
fact the BV problem, when solved by common finite differences, appears ill-posed.
6.2. Reformulation of PDE to form appropriate for using

- regularization. Numerical results

Extension of - regularization to solve PDE requires the combined use of the BC (3.11) on the front
with the SE (3.12) to impose a BC at a small relative distance from the front. The distance being
relative, we need to count it in the local system with the origin at the front. Hence it is reasonable to
( x x) / x from the front. The relative distance from the inlet is
introduce the relative distance
/ t x const of a function ( x, t ) , which enters PDE
1
x / x (t ) . The partial time derivative
and is evaluated under constant x, should be transformed into the partial time derivative

/ t

const

( x (t ), t ) evaluated under constant . Omitting routine details of the change
of the function ( , t )
of variables, we have the transformation:
v (t )
,
(6.1)
t x const
t const
x (t )

where v (t ) dx / dt . When using the variable

and transformation (6.1), equation (3.13) becomes:

13
2

Y
2

where Y ( , t )

A(Y , Y /

w3 ( x (t ), t ) , A(Y , Y /

Y

,x v )

,x v )

B(Y , x )

Y/

Y
t

0.75 x v
3Y

0,

, B(Y , x )

(6.2)

x2
.
4Y

The BC (3.10), (3.11) in new variables read:

43Y Y
3 x ς
Y ( ,t)

q0 ,

(6.3)

ς 0
1

0.

(6.4)

The SE (3.12) takes the form:
Y

0.75 x v .

(6.5)

1

Prove firstly, that the problem (6.2)-(6.4) is ill posed because like the self-similar formulation, the
SE (6.5) is met identically by a solution of PDF (6.2) under the BC (6.4) of zero opening at the liquid
front. Indeed, re-write (6.2) by using the expression for A(Y , Y / , x v ) as
2

1 Y
Y
Y
(6.6)
0.75x v
0.75x 2
0.
3
t
x ), for a solution, satisfying the BC (6.4), PDE (6.6) turns into the SE (6.5).
In limit
1 (x
Hence, for PDE (6.2), at the point
1 , we have imposed not only the BC (6.4) for unknown function
Y , but also the BC (6.5) for its spatial derivative dY / d . Therefore, we have two rather than one BC
at
1 and the problem appears ill-posed. Consequently, the starting problem (3.9)-(3.11) is ill-posed,
as well, what explains the failure to solve it to the accuracy greater than two correct digits. Note that
equations (6.4), (6.5) imply that the factor A(Y , Y / , x v ) in (6.2) is finite at the liquid front
despite its denominator 3Y turns to zero: lim A(Y , dY / d , x*v* )
1/ 3 .
Y

Y
2

The regularization of the problem (6.2)-(6.4) follows the line used for the self-similar formulation.
Like the self-similar formulation, the BC (6.4) and the SE (6.5) yield the approximate equation near
the liquid front
1:
Y ( , t ) 0.75 x (t )v (t )(1 ) ,
(6.6)
which defines the asymptotic behavior of the solution when
1 . As mentioned, the non-linear
multiplier A(Y , Y / , x v ) in PDF (6.2) is non-singular at the liquid front
1 . The factor
B (Y , x ) is finite except for the liquid front (
1 ), because the opening is positive behind the front.
Hence, similar to (5.1), we may impose the BC at the relative distance from the liquid front:
(6.7)
Y ( , t ) 0.75x (t )v (t ) ,
where
1 . In contrast with the problem (6.2)-(6.4), the problem (6.2), (6.3), (6.7) does not
involve an additional BC. We may expect that it is well-posed and provides the needed regularization.
Extensive numerical tests confirm the expectation.
We solved the problem (6.2), (6.3), (6.7) by using the Crank-Nicolson scheme and iterations for
non-linear multipliers A(Y , Y / , x v ) , B (Y , x ) at a time step. The velocity v (t ) is also iterated
by using the equation following from (6.6):
Y
3
xv .
(6.8)
4
The condition (6.8) expresses (with an accepted tolerance) the continuity of the particle velocity at
the point
. After completing iterations, the final value of v (t ) shows the new coordinate of the
t ) x (t ) v (t ) t . The value x (t
t ) is used on the next time step.
liquid front x (t

14
Initialization. In the considered problem there is no characteristic time and length. Thus, the
starting equation x (0) 0 should be adjusted to this uncertainty. The adjustment concerns the
initialization of time stepping. At t = 0 the liquid front coincides with the inlet. The opening is also
zero. Then for „small‟ time, to have the factor A(Y , Y / , x v ) finite in the entire liquid, it should
x v f ( ) with f ( )
be Y /
0.75 as
1 . Hence for „small‟ time, we have Y x v f ( )
with f (1) 0 . Substitution of Y and Y /
into the BC (6.3), using v dx / dt and integration of
0.6

the resulting ODE yield x (t )

Ct

0.8

, where C 1.25

0.75q0

0.8
3

f (0) f (0)

is a constant. Then

v (t ) 0.8Ct 0.2 and Y 0.8C 2t 0.6 f ( ) . Insertion of these equations into PDF (6.2) gives ODE (3.14)
/ . Similarly, the BC (6.3), (6.4) and the SE (6.5) turn into the
with y 0.8 f and
corresponding equations (3.15), (3.16) and (3.17). The constant C actually represents
,
corresponding to the prescribed flux q0 .
We see that when using time stepping, initialization requires solving the self-similar problem. The
latter being ill-posed, the solution is obtained by - regularization as explained in Sec. 5.2. Having its

solution, the initial data for an arbitrary chosen initialization time t 0 are found as x (t 0 )

t 00.8 ,

2 0.6
0.8 t 0 0.2 , Y ( , t0 )
t0 ( * ) . For certainty, in the following calculations we set q0 =
1, t0 = 0.01. The same regularization parameter
and the same spatial step
were used for both
the initialization and for each of the time steps.
Numerical results. Two objectives were sought in numerical tests. First, we wanted to check the
efficiency of -regularization, that is, its accuracy, stability and robustness. To this end, we used small
relative distance , small spatial step
, small time step t , and a large number of time steps.
Secondly, we checked if the beneficial features of coarse meshes, observed for self-similar
formulation, hold in time steps. Exploratory calculations have shown that fifty iterations in non-linear
terms at a time step are sufficient to reproduce seven significant digits. Thus, in further tests the
number of iterations was set equal to fifty. The benchmark solution served to evaluate the accuracy. In
particular, for q0 = 1, the benchmark values of the front position and the front velocity are

v (t 0 )

x (t ) 1.3208446 t 0.8 and v (t ) 1.0566757 t 0.2 , respectively.
Results for fine meshes and small time steps. Obviously, the accuracy at the first time steps
depends on the accuracy of the data obtained at the initialization stage. We could see that = 0.0001
and
= 0.01 provide the accuracy of 0.026% for . To keep the accuracy on this level for x (t )
and v (t ) at first time steps, the step t should be notably less than t 0 . By taking t = 0.01 t 0 = 10-4,
we had the accuracy of some 0.03% for both x (t ) and v (t ) at first hundred steps. With further
growth of the time, the accuracy of v (t ) stayed on the same level, while the relative error of x (t )
decreased. In particular, at the step m = 1000 (t = 0.11), it is 0.0077%; at the step m = 20000 (t = 2.01)
it is 0.0043%. The decrease of the relative error in x for large time is related to the growth of the

absolute value of x (Recall that x (t )
t 0.8 .). The time expense on a conventional laptop for 20000
steps with 50 iterations at a step is near 15 s. Note that in view of the growth of x (t ) in time, there is
no need to have the time step constant. By using exponentially growing time steps, the number of time
steps may be drastically reduced without loss of the accuracy. For instance, with the same time
expense we could reach time t = 36128 ( x = 5889.4) with no loss of the accuracy. The number of
iterations could be reduced as well.
There were no signs of instability in these and many other specially designed experiments. We
10 4 and
could see that for 10 2
0.01, the results are accurate and stable in a wide range of

15
values of the time step and for very large number (up to 100000) of steps. Therefore - regularization
(6.7) is quite efficient.
Results for coarse mesh and large number of time steps. For the coarse mesh
= 0.0(9) and =
0.1 at the initialization stage, we have one percent accuracy except for close vicinity of the liquid front.
There is no need to take time steps notably less than the initial time. When taking t = 0.25 t0 =
0.0025, we had x with the error not exceeding 1.2% at any time instant for 20000 equal steps. Again,
the relative error decreased with growing x . Using exponentially growing time step serves to obtain
x (t ) and v (t ) with the accuracy of one percent even when x (t ) becomes of order 10000 (recall that
the starting value is x (t0 ) 0.033170). As the number of nodes is small (only 11), the time expense is
fractions of a second. Again, there were no signs of instability. Moreover, increasing the initial time
step 3-10-fold, influenced the accuracy only of the first steps, while with growing time, the accuracy
kept to the mentioned level of about 1%. Hence, use of a coarse mesh efficiently maintains this
accuracy.
7. Extension of the regularization method to 2D fracture propagation
According to the rationale presented in the preceding section, it appears that the strategy of using regularization when tracing 2-D hydrofracture propagation is as follows. At each point of the liquid
front we introduce the local coordinate system moving with the front in the direction normal to the
front. We transform field equations to this system. The prescribed BC at the front and the SE
associated with front points are transformed accordingly. They are combined to formulate a new BC at
a small distance behind the front. Below, we follow this path.
We introduce the local Cartesian coordinates x1 ' O1 ' x2 ' with the origin O1 ' at a point x* at the
liquid front and the axis x1 ' opposite to the direction of external normal to the front at this point. As
the flux q* at the front is co-linear to the normal, its tangential component is zero. For the point x* ,
equation (2.11) written in the system x1 ' O1 ' x2 ' becomes a scalar equation qn*

D ( w, p )

p
x1'

,
x1 ' 0

and the SE (2.12) reads;
1
p
D( w, p)
w*
x1'

v

where v
obtain:

,

(7.1)

x1 ' 0

dxn* / dt is the absolute value of the front velocity. Combining (7.1) with the BC (2.16) we
p

1
D( w, p)dp v* x1 ' .
p0 w
Equation (7.2) serves to impose the BC at a small distance x1 '
p

(7.2)
from the front:

1
(7.3)
D( w, p)dp v* ,
p0 w
where p
. The - regularization consists in using the BC
p( ) is the pressure at the point
(7.3) instead of the BC (2.16).
In numerical calculations with - regularization, it is reasonable to use the aforementioned
advantages of variables which express the particle velocity v q / w . In the local system x1 ' O1 ' x2 '
moving with the front, for points close to the front, the normal component of the velocity is:
1
p
. Evaluating the partial time derivative under constant x1 ' rather than at constant
v
D(w, p)
w
x1 '

16
global coordinate, the PDF (2.9) at points near the front takes the form (2.17). As noted in Sec. 2, it is
reasonable to transform PDF (2.17) to (2.18) by introducing the variable y w1 / , where the
exponent α accounts for asymptotic behavior of the opening at the front.
As
is small, the velocity v( ) is close to v* . This serves to employ the equation
(7.4)
v( ) v *
1 / 3 , p 4w , p0 4w( x* ) 0 ,
for iterations in v* at a time step. For the Nordgren problem,
x* ; then equations (2.18), (7.3) and (7.4) are easily transformed to (6.2), (6.7) and (6.8),
respectively.
8. Conclusions
The paper presents the following conclusions.
(i)
The speed function for fluid flow in a thin channel is given by the ratio of the flux through a
cross section to the channel width at the front. Its specification for hydraulic fracturing
follows from the equation of Poiseuille type connecting the flux with gradient of pressure.
The speed function, as the basis of the theory of propagating interfaces, facilitates
employing such methods as level set methods and fast marching methods. The speed
equation for hydraulic fracturing is a general condition at the liquid front not dependent on
a particular BC defined by a physical situation ahead of the front, liquid properties, leak-off
through the channel walls and/or presence of distributed sources. When there is no lag
between the liquid front and the crack tip, both the flux and opening are zero at the liquid
front; in this case, the SE is fulfilled in the limit when a point behind the front tends to the
front.
(ii)
Using the SE gives a key to proper choice of unknown functions when solving a hydraulic
fracture problem numerically. Specifically, the particle velocity, averaged over the opening,
is a good choice, because it is non-singular at the front and non-zero in entire flow region.
Using the velocity as an unknown may also serve to avoid the stiffness of the system of
differential equations obtained in a conventional way. Another proper unknown is the
opening to the degree 1/α, where α is non-negative exponent, characterizing the asymptotic
behavior of the opening near the front.
(iii) Existence of the SE also discloses the crucial feature of the problem: for zero or small lag,
at the points of the front we actually have prescribed two rather than one BC; the BV
problem appears ill-posed when solved numerically. It requires a proper regularization to
have accurate and stable numerical results.
(iv)
Self-similar formulation of the Nordgren problem, reducing PDE to ODE, unambiguously
demonstrates that the problem is ill-posed when considered as a BV problem. Numerical
experiments confirm this theoretical conclusion: by no means could we obtain more than
two correct digits when solving the problem without regularization, and the results always
strongly deteriorate near the liquid front.
(v)
The solution of the self-similar problem, obtained when solving it as a Cauchy (initial
value) problem by the Runge-Kutta method, provides benchmarks with at least five correct
digits at entire liquid including its front. These numerical results may serve for testing
methods of hydraulic fracture simulation.
(vi)
Studying of the Nordgren problem, accounting for the SE, suggests a means to overcoming
the difficulty by - regularization. It consists in employing the SE together with a
prescribed BC on the front to formulate a new BC at a small distance behind the front
rather than on the front itself. Application of - regularization and comparing the results
with the benchmarks shows that it is robust and provides highly accurate and stable results.
(vii) Numerical experiments also disclose that using coarse spatial meshes is beneficial and
serves as a specific regularization when high accuracy is not needed. Still, although

17
involving analytical work, - regularization is superior in the possibility to guarantee
accurate and stable results and to evaluate the accuracy of calculations; it is competitive in
time expense.
(viii) The suggested - regularization may serve as a means to obtain accurate and stable results
for simulation of hydrofracture with zero or small lag. When tracing 2-D hydrofracture
propagation, - regularization is obtained by writing equations in the local coordinate
system moving with the liquid front.
Acknowledgment. The author appreciates the support of the EU Marie Curie IAPP program (Grant #
251475).
References
Adachi J.I. & Detournay E. 2002 Self-similar solution of plane-strain fracture driven by a power-law
fluid. Int. J. Numer. Anal. Meth. Geomech., 26, 579-604.
Adachi J., Siebrits E. et al. 2007 Computer simulation of hydraulic fractures. Int. J. Rock Mech.
Mining Sci., 44, 739-757.
Batchelor, G.K. 1967 An Introduction to Fluid Dynamics. Cambridge University Press.
Carter R.D. 1957 Derivation of the general equation for estimating the extent of the fractured area,
Appendix to “Optimum fluid characteristics for fracture extension” by G.C. Howard and C.R.
Fast, Drill. and Prod. Prac. API, 261-270.
Crowe C.T., Elger D.F., Williams B.C. & Roberson J.A. 2009 Engineering Fluid Mechanics, 9th ed.,
John Willey & Sons, Inc.
Descroches J., Detournay E. et al. 1994 The crack tip region in hydraulic fracturing. Proc. Roy Soc.
London, Ser. A, 447, 39-48.
Detournay E., Adachi J. & Garagash D. 2002. Asymptotic and intermediate asymptotic behaviour near
the tip of a fluid-driven fracture propagating in a permeable elastic medium. Structural Integrity
and Fracture, Dyskin A., Hu X. and Sahouryeh E. (eds). The Netherlands: Swets & Zeitlinger, 918.
Epperson J.F. 2002 An Introduction to Numerical Methods and Analysis. 556p. New York: John Wiley
& Sons, Inc.
Garagash D.I. 2006 Propagation of a plane-strain hydraulic fracture with a fluid lag: Early time
solution. Int. J. Solids Struct. 43, 5811-5835.
Garagash D.I. & Detournay E. 2000 The tip region of a fluid-driven fracture in an elastic medium.
ASME J. Appl. Mech. 67, 183-192.
Garagash D.I., Detournay E. & Adachi J.I. 2011 Multiscale tip asymptotics in hydraulic fracture with
leak-off. J. Fluid Mech. 669, 260-297.
Geertsma J. & de Klerk F. 1969 A rapid method of predicting width and extent of hydraulically
induced fractures, J. Pet. Tech. 21, 1571-1581.
Hadamard J. 1902 Sur les problemes aux derivees partielles et leur signification physique. Princeton
University Bulletin, 49-52.
Howard G.C. & Fast C.R. 1970 Hydraulic Fracturing. Dallas: Monograph Series Soc. Petrol. Eng.
Hu J. & Garagash D.I., 2010 Plane strain propagation of a fluid-driven crack in a permeable rock with
fracture toughness. ASCE J. Eng. Mech., 136, 1152-1166.
Jamamoto K., Shimamoto T. & Sukemura S. 2004 Multi fracture propagation model for a threedimensional hydraulic fracture simulator. Int. J. Geomech. ASCE, 1, 46-57.
Khristianovich S.A.& Zheltov V.P. 1955 Formation of vertical fractures by means of highly viscous
liquid. Proc. 4-th World Petroleum Congress, Rome, 579-586.
Kovalyshen Y. & Detournay E. 2009 A re-examination of the classical PKN model of hydraulic
fracture. Transport in Porous Media, 81, 317-339.
Kovalyshen Y. 2010 Fluid-Driven Fracture in Poroelastic Medium. PhD Thesis. 215 p. Minnesota
University.

18
Lavrent‟ev M.M. & Savel‟ev L.Ja. 1999 Theory of Operators and Ill-posed Problems. (Novosibirsk:
Institute of Mathematics im. S. L. Sobolev, ), ISBN 5-86134-077-3 [in Russian].
Lenoach B. 1995 The crack tip solution for hydraulic fracturing in a permeable solid. J. Mech. Phys.
Solids, 43, 1025-1043.
Linkov A.M. 2011a Speed equation and its application for solving ill-posed problems of hydraulic
fracturing. Doklady Physics, 56, 436-438.
Linkov A.M. 2011b On numerical simulation of hydraulic fracturing. Proc. XXXVIII Summer SchoolConference, Advanced Problems in Mechanics, Repino, St. Petersburg, 291-296.
Mitchell S.L., Kuske R. & Pierce A.P. 2007 An asymptotic framework for analysis of hydraulic
fracture: the impermeable fracture case. ASME J. Appl. Mech., 74, 365-372.
Nolte K.G. 1988 Fracture design based on pressure analysis. Soc. Pet. Eng. J., Paper SPE 10911, 1-20.
Nordgren R.P. 1972 Propagation of a vertical hydraulic fracture, Soc. Pet. Eng. J., August, 306-314.
Perkins K. & Kern L. F. 1961 Widths of hydraulic fractures, J. Pet. Tech., 13, 937-949.
Pierce A.P. & Siebrits E. 2005 A dual multigrid preconditioner for efficient solution of hydraulically
driven fracture problem. Int. J. Numer. Meth. Eng., 65, 1797-1823.
Savitski A. & Detournay E. 2002 Propagation of a fluid driven penny-shaped fracture in an
impermeable rock: asymptotic solutions. Int. J. Solids Struct., 39, 6311-6337.
Sethian J.A. 1999 Level Set Methods and Fast Marching Methods. Cambridge: Cambridge University
Press, 1999.
Spence D.A. & Sharp P.W. 1985 Self-similar solutions for elastohydrodynamic cavity flow. Proc. Roy
Soc. London, Ser. A, 400, 289-313.
Tychonoff A.N. 1963 Solution of incorrectly formulated problems and the regularization method.
Soviet Mathematics, 4, 1035-1038 [Translation from Russian: А.Н. Тихонов, Доклады АН
СССР, 1963, 151, 501-504].

19

Well
Hydraulic fracture

h
x

w(x,t)
O
x*(t)
y

Fig. 1 Scheme of the problem on hydraulic fracture propagation

dψ 3 /dξ
0.6

0.4

0.2

0
0

0.2

0.4

Fig. 2 Graphs of the normalized velocity d

0.6

3

/d

0.8

1 ξ

in self-similar formulation of

Nordgren‟s problem
benchmark solution
solution of ill-posed boundary value problem

