## Abstract

This paper develops a closed-form approximate solution for a penny-shaped hydraulic fracture whose behaviour is determined by an interplay of three competing physical processes that are associated with fluid viscosity, fracture toughness and fluid leak-off. The primary assumption that permits one to construct the solution is that the fracture behaviour is mainly determined by the three-process multiscale tip asymptotics and the global fluid volume balance. First, the developed approximation is compared with the existing solutions for all limiting regimes of propagation. Then, a solution map, which indicates applicability regions of the limiting solutions, is constructed. It is also shown that the constructed approximation accurately captures the scaling that is associated with the transition from any one limiting solution to another. The developed approximation is tested against a reference numerical solution, showing that accuracy of the fracture width and radius predictions lie within a fraction of a per cent for a wide range of parameters. As a result, the constructed approximation provides a rapid solution for a penny-shaped hydraulic fracture, which can be used for quick fracture design calculations or as a reference solution to evaluate accuracy of various hydraulic fracture simulators.

## 1. Introduction

Hydraulic fractures are the fluid-filled cracks that propagate under the influence of a fluid pressure acting along the crack’s surface. The most common and well-known application of hydraulic fracturing is the stimulation of oil and gas wells to enhance production of hydrocarbons [1]. Other industrial applications include waste remediation process [2], waste disposal [3] and preconditioning in rock mining [4]. Hydraulic fractures also occur in nature in the process of magma ascent through the lithosphere owing to a buoyancy force [5–8,9] or as fluid-filled cracks in glacier beds [10].

Different hydraulic fracture geometries have been considered over the years. Research efforts shifted from the development of simple models, such as the Khristianovich–Zheltov–Geertsma–De Klerk (KGD) model for a plane-strain crack [11], towards more complex models that consider a planar hydraulic fracture [12–14], multiple hydraulic fractures [15–17] or fracture networks [18]. Detailed reviews of various hydraulic fracturing models can be found in [19–21]. In addition, there are techniques in which cracks are not modelled explicitly, such as the phase field [22,23], distinct element method [24,25] and peridynamics [26]. The primary advantage of such methodologies is the ability to tackle complex crack geometries easier than the conventional methods. At the same time, predictions of such approaches should be thoroughly tested against reference solutions to ensure that the new techniques are able to capture all the features of the conventional approaches.

Even for the simplest geometries, hydraulic fractures are known to obey a complex multiscale behaviour, see e.g. a thorough review paper [27]. This multiscale nature arises both in time, where multiple timescales determine fracture evolution, and space, where the solution may experience variations at different length scales in the tip region. As pointed out in [27], the time and length scales are related; that is, a particular timescale in the fracture evolution corresponds to dominance of one length scale in the tip region. Recognizing the significance and multiscale nature of the tip region, many studies were devoted specifically to quantify hydraulic fracture behaviour in the tip region [28–35]. On the other hand, time evolution and regimes of propagation were studied for a plane-strain KGD fracture in [36–40], and for a penny-shaped hydraulic fracture in [40–44]. A recent review paper [27] provides a more detailed summary of the findings and shows the complexity of the behaviour that occurs even for simple fracture geometries.

In view of the multiscale behaviour of hydraulic fractures, the primary goal of this paper is to quantify such a behaviour for the case of a penny-shaped hydraulic fracture, where the latter is driven by a Newtonian fluid and propagates in a permeable medium assuming no fluid lag. Most of the previous studies that addressed the problem of a penny-shaped fracture focused on the limiting regimes of propagation and asymptotic closed-form solutions (or accurate approximations) for the problem [40–43]. One exception is the study [44], in which a numerical solution for the full problem is obtained. The numerical solution is, however, relatively difficult to obtain owing to the temporal and spatial multiscale nature of the solution. In contrast, this study develops a closed-form approximate solution that provides results virtually instantly and accurately captures the complex behaviour of a radial hydraulic fracture at all length scales and timescales. In particular, the developed solution is able to describe all existing limiting solutions and all possible transitions between them, so that it covers all the trajectories in the parametric space for the problem under consideration. This development is made possible by using a closed-form approximation for the tip asymptotic solution obtained in [34], which is used to approximate fracture width profile. Once the fracture geometry is known, the global fluid volume balance is used to determine behaviour of the solution.

The importance of the obtained solution can be summarized as follows. It first shows that the tip region plays a crucial role in hydraulic fracture modelling. It also allows one to obtain a solution rapidly, which can be useful for quick estimations of fracture geometry for any values of fracture toughness, fluid viscosity and leak-off. Owing to a relatively simple implementation of the solution, it can be used as a reference solution to evaluate accuracy of other hydraulic fracturing simulators and, at the same time, as an initial condition to improve stability of the numerical schemes at early times. Finally, the developed approximation leads to construction of the solution map, which indicates validity regions of the limiting solutions and permits one to easily determine whether the solution for given problem parameters corresponds to one of the limiting cases.

This paper is organized as follows. Section 2 describes governing equations for a radial hydraulic fracture with leak-off. Section 3 outlines procedure for obtaining the approximate solution. Comparison of the developed approximation with the existing limiting solutions is presented in §4. Section 5 contains description of the solution structure, where applicability zones of the limiting solutions are indicated. Finally, §6 evaluates accuracy of the approximation by comparing its predictions to the reference numerical solution, which is followed by the summary.

## 2. Governing equations for a penny-shaped hydraulic fracture

This study considers propagation of an axisymmetric (‘radial’ or ‘penny-shaped’ are also used throughout the paper) hydraulic fracture in a permeable rock [27,44]. There are four primary material parameters that appear in the model, which can be introduced in the scaled form, for convenience, as
*μ* is the fluid viscosity, *E* is the Young modulus, *ν* is the Poisson’s ratio, *K*_{Ic} is the mode I fracture toughness of the rock and *C*_{L} is the Carter’s leak-off parameter.

Volume balance for an incompressible Newtonian fluid inside the crack can be written as
*w*(*r*,*t*) denotes the fracture width, *q* is the flux in the radial direction, the term that contains *C*′ accounts for leak-off via Carter’s model, *t*_{0}(*r*) is the time instant at which the fracture front was located at point *r*, *p* is the fluid pressure and *Q*_{0} is the fluid injection rate (assumed to be constant in time).

The elasticity equation, which characterizes elastic equilibrium of the rock, relates the fluid pressure inside the crack to the fracture width as [41,44,45]
*R* is the fracture radius, and the kernel is
*K*(⋅) and *E*(⋅) denote the complete elliptic integrals of the first and the second kind, respectively.

Fracture propagation is modelled by a classical result of the linear elastic fracture mechanics (LEFM), in which the fracture opening in the tip region follows the square root solution [46]
*q*(*R*,*t*)=0.

For future reference, it is useful to consider the global fluid balance, which can be obtained by integrating (2.2) with respect to time and space as
*q*(*R*,*t*)=0 and *w*(*R*,*t*)=0 were used to obtain the result.

## 3. Approximate solution for a penny-shaped hydraulic fracture

### 3.1. Outline of the methodology

To construct the approximate solution, it is assumed that the fracture evolution is primarily determined by the near-tip behaviour and the global fluid balance (2.6). In this situation, the approximation for the fracture width can be taken in the form
*w*_{a} is the tip asymptotic solution that, in addition to distance to the tip *R*−*r*, depends on the material parameters (2.1) and time through *R*(*t*) and *w*_{a}, is the solution of (2.2), (2.3) and (2.5) in the tip region defined as (*R*−*r*)/*R*≪1, which, as shown in [13], is equivalent to the problem of a semi-infinite hydraulic fracture that propagates steadily under plane-strain elastic conditions [32,34]. As a result, because *w*(*R*−*s*,*t*)=*w*_{a}(*s*) for *s*≪*R*, the approximation (3.1) precisely solves the governing equations (2.2), (2.3) and (2.5) in the tip region and approximates the solution inside the fracture away from the tip.

In order to construct the approximation, the closed-form approximation for *w*_{a}, obtained in [34], is employed. As shown in [34], the three-process tip asymptotic solution, which captures the effects of fracture toughness, fluid viscosity and leak-off, has the maximum error of 0.14% for all regimes of propagation and all possible transitions between them. One of the properties of the tip solution, *w*_{a}, is that *s*=*R*−*r* is the distance from a point inside the fracture to the tip and *R*(*t*)∝*t*^{α} and *α* is a number that is equal to either 1/4, 2/5, or 4/9 [44,27], it is further assumed that *R*(*t*)∝*t*^{α}, where *α* is a slowly varying function. With this approximation, the function *t*_{0}(*r*) is determined from the relation *r*/*R*=(*t*_{0}/*t*)^{α}, which together with (3.2) reduces the global fluid balance (2.6) to
*ρ*=*r*/*R* is the scaled spatial coordinate. The integrals in equation (3.3) can be evaluated, in which case (3.3) reduces to
*B*(*a*,*b*) denotes the beta function, and
*B*(*x*;*a*,*b*) is the incomplete beta function.

Procedure for obtaining the approximate solution can be outlined as follows. First, material parameters (2.1), *Q*_{0} and the function *w*_{a}(*R*) should be specified. Note that *w*_{a}(*R*) also depends on the material parameters (2.1) and *α*=4/9 (corresponds to zero toughness and zero leak-off solution), equation (3.4) can be solved for *R*(*t*) (e.g. by using Newton’s method). Once *R*(*t*) is obtained, the value of *α* is updated according to *α*. Once *R*(*t*) is calculated, the width can be calculated using (3.2). Efficiency, which is defined as the ratio between the current fracture volume and the total amount of injected fluid, can be calculated as
*M*(*ρ*,*s*) is given in (2.4). The accuracy of such pressure calculation is examined later in the paper.

It should be noted here that λ is taken as a parameter at this point, and its value can be varied to achieve a better approximation. At the same time, in order to capture the elliptic solution for a uniformly pressurized crack, one should set λ=1/2. As is shown later, the values of λ close to 1/2 provide the most accurate results.

### 3.2. Solution in scaled variables

To solve (3.4) numerically, it is useful to introduce scaling, in which the tip asymptotic solution *w*_{a} is expressed. In particular, as shown in [14,17], the solution *w*_{a} is given implicitly by the equation
*g*_{δ} is a known function (see appendix A), and *α* is taken as 4/9, which is then updated iteratively using *α* converges).

Once the time histories

## 4. Comparison with vertex solutions

The problem of a penny-shaped fracture with fluid loss and no lag has four limiting regimes of propagation (or vertex solutions) [44,27]: the storage viscosity (*M*), storage toughness (*K*), leak-off viscosity (*η*≤1 as
*η*=1 corresponds to the situation in which all fluid is stored in the fracture, whereas *η*=0 signifies the case in which all the injected fluid leaks into the formation. It is important to note that *η*=1 corresponds to *η*=0 corresponds to

To evaluate accuracy of the approximate solution, and to specify the values for λ, the remainder of this section is devoted to a comparison between the limiting regimes (4.1) of the approximate solution and the existing reference solutions for these cases. In addition, setting some of the parameters, such as *K*′ and *C*′ (or *C*′. This fact will be used for the remaining part of this section to write the limiting solutions in terms of scaling (3.8) and (3.9).

### 4.1. *M* vertex limit of the solution

Before proceeding with the *M* vertex limit, it is useful to rewrite equations (3.9) and (3.12) in the form
*M* vertex solution corresponds to the case of no toughness and no leak-off, it is possible to set *α*=4/9. Note that *β*_{m}=2^{1/3}3^{5/6}.

By introducing scaling that is associated with the *M* vertex solution [13],41,44]
*α*=4/9 (because *L*_{m}∝*t*^{4/9}) is used to obtain the result. The solutions for the fracture width and pressure variations can be written using (3.13) as
*Ω*_{m}(0) and fracture radius *γ*_{m} that are calculated using (4.4) and (4.5) with those obtained in [41] for the *M* vertex case implies that λ_{m}=0.487 provides the best approximation, for which the error for both *Ω*_{m}(0) and *γ*_{m} is below 0.5%.

The *M* vertex limit of the approximate solution can be summarized as

### 4.2. M ∼ vertex limit of the solution

In order to obtain the *η*=0. In this situation, equation (3.14) implies that
*g*_{δ} at large

By introducing scaling that is associated with the

The

### 4.3. *K* vertex limit of the solution

The *K* vertex solution corresponds to the case when *η*=1 (or *α*=2/5.

By introducing scaling that is associated with the *K* vertex solution [13,41,44]
*K* vertex solution is an ellipse [41], then λ_{k}=1/2 leads to the exact solution.

Given that *K* vertex solution in [41,44]. The above results can be written in the unscaled form as
*K* vertex limit is an ellipse and that the approximation (3.1) captures elliptic width profile precisely, the obtained limited solution is the exact solution for the *K* vertex limit and should not be tested against a reference solution.

### 4.4. K ∼ vertex limit of the solution

The *η*=0. As a result, the efficiency relation (3.14) reduces to

By introducing scaling that is associated with the *K* vertex limit, *K* vertex case, the fracture width for the

Because *K* vertex solution, the

### 4.5. Parameter λ interpolation

This section aims to describe an interpolation scheme to determine values of the parameter λ, knowing its values for the limiting cases. Comparison between the developed approximate solution and the vertex solutions indicates that λ_{m}=0.487, _{k}=0.5 and *g*_{δ} depends on *M* to *K* (*η*=0) is less important for the interpolation purposes, because λ_{k}≈λ_{m} and any interpolation would work. The transitions from *M* to *K* to

## 5. Structure of the solution

Solution in the original scaling is given by the quantities *η*, which depend on the dimensionless time *ϕ* is defined as
*γ*(*τ*), fracture width *Ω*(*ρ*,*τ*), fluid pressure *Π*(*ρ*,*τ*) and efficiency *η*(*τ*).

To visualize the structure of the solution, figure 1 plots the approximate width solution evaluated at the wellbore *Ω*(0,*τ*) versus *τ* for a range of values *ϕ*. Blue, red, green and magenta regions indicate, respectively, validity zones of the *M* vertex solution (4.6)–(4.7), *K* vertex solution (4.20)–(4.21), *i* corresponds to either *M*, *K*,

Figure 1 is a map of the solutions and allows one to determine the regime of propagation of a radial fracture for a given problem parameters. For instance, for the parameters *E*′=10 GPa, *K*′=1 MPa⋅m^{1/2}, *Q*_{0}=10^{−2} m^{3} s^{−1}, *t*=10^{3} s, *μ*′=1 Pa⋅s and *C*′=10^{−6} m s^{−1/2}, the dimensionless time and leak-off parameters, respectively, are *τ*=10^{−5} and *ϕ*=1, in which case, according to the figure 1, the *M* vertex solution still can be used to estimate the fracture parameters. At the same time, if the viscosity is changed to *μ*′=0.01 Pa⋅s (other parameters are kept the same), then *τ*=1 and *ϕ*=10^{−6}, so that none of the vertex solutions can be used. Instead, the developed approximation can be used. In addition, figure 1 can be used to indicate possible time evolution paths. The fracture starts at the *M* vertex and ultimately reaches the *K* vertex, through the

*Transition along the MK edge.* To quantify the boundaries between the *M* and *K* vertices, it is necessary to consider the *MK* edge transition. This *MK* edge solution corresponds to the situation of no leak-off, or *η*=1. As a result, *MK* edge solution. It is important to note that the solution depends on the parameter *τ*. As a result, the dimensionless time that quantifies the *MK* transition is
*MK* transition, shown in figure 1, are calculated numerically as
*MK* transition (near the *M* vertex), whereas *K* vertex).

*Transition along the* *edge.* To quantify the boundaries between the *M* and *τϕ*^{9/14}. In this case, the dimensionless time that quantifies the *M* vertex), whereas

*Transition along the* *edge.* To quantify the boundaries between the *K* and *τϕ*^{5/6}. In this case, the dimensionless time that quantifies the *K* vertex), whereas

*Transition along the* *edge.* To quantify the boundaries between the *η*=0 or *τϕ*^{−1/2}. In this case, the dimensionless time that quantifies the

## 6. Comparison with numerical solution

To check the accuracy of the developed approximate solution, it is necessary to compare predictions of the approximation to a reference solution. In this study, the reference solution is calculated numerically by solving the original equations governing the problem of a radial hydraulic fracture (2.2) and (2.3). One of the key features of the scheme is the ability to capture multiscale tip behaviour by using the tip asymptotic solution as a propagation condition. Brief description of the numerical scheme is given in appendix B.

Because §4 examined accuracy of the approximation at the limiting cases of vertex solutions, the aim of this section is to check the accuracy in the transition regions. For this reason, the choice of the parameters for comparison is *ϕ*={10^{−18},10^{−15},10^{−12},10^{−9},10^{−6},10^{−3},1, 10^{3},10^{6}} and 10^{−7}<*τ*<10^{7}, which cover the transition region according to figure 1. It is important to note that the numerical solution is more difficult to obtain near the *K* and *K* and

Figure 2 shows comparison between the developed approximate solution (dashed grey lines) and the numerical solution (solid black lines) in terms of time histories of the fracture radius (panel (*a*)), efficiency (panel (*b*)), width at the wellbore (panel (*c*)) and pressure at *ρ*=0.5 (panel (*d*)), all for different values of *ϕ* (arrows indicate the direction in which the curves shift as *ϕ* increases). The fracture width and pressure are evaluated as
*M* vertex solution (4.6), which can be written in terms of the scaling (5.1) and (5.3) as
*M* vertex), red (*K* vertex), green (*M* vertex (i.e. the blue lines), and then transition to either *K*, *ϕ*=0 case) reach the

As can be seen from figure 2, there is no visual difference between the numerical solution and the approximation, which indicates that the developed approximate solution is accurate. Despite this fact, it is still important to quantify the level of accuracy of the approximation. To help doing this, the error is defined as
*A* is either radius *γ*(*τ*), efficiency *η*(*τ*), wellbore width *Ω*^{0}=*Ω*(0,*τ*) or pressure *Π*^{0}=*Π*(0.5,*τ*).

Figure 3*a* shows the error that is calculated using (6.2) for different values of the parameter *ϕ* for the finest mesh (for numerical solution) considered. The magnitude of the error for *γ*, *Ω*^{0} and *Π*^{0} does not vary notably versus *ϕ*, which indicates that the interpolation procedure for λ (4.29) is sufficiently accurate. The accuracy of all quantities, but the pressure, is under a per cent. The pressure is less accurate because it is calculated using the elasticity integral (2.3), which is very sensitive to the fracture width profile. To better understand the effect of mesh size that is used to calculate numerical solution, figure 3*b* shows the maximum error versus the number of spatial discretization points of the numerical solution *N*. Here the maximum is calculated for different values of *ϕ*, which can be obtained from figure 3*a* and its analogues for different meshes. Results in figure 3*b* demonstrate that error for the fluid pressure *Π*^{0} and fracture length *γ* reached the plateau values, but the efficiency *η* and the wellbore width *Ω*^{0} did not reach their respective plateau values for *N*=200. Despite this fact, it is still clear that the error of the developed approximation is under a per cent (for the considered error measure), which is sufficient for most practical cases.

In order to check accuracy of the pressure and width spatial variations, figure 4 plots pressure and width predictions stemming from the approximate solution (dashed grey lines) and those calculated numerically (solid black lines) for *ϕ*=10^{−15} (panel (*a*)), *ϕ*=10^{−6} (panel (*b*)) and *ϕ*=10^{3} (panel (*c*)). Both the pressure and the width are normalized by the *M* vertex solution (6.1) for convenience. The vertex solutions are shown by the blue (*M* vertex), red (*K* vertex), green (*τ*. For the *ϕ*=10^{−15} case (panel (*a*)), the solution transitions from the *M* vertex and almost reaches the *K* vertex. For the *ϕ*=10^{−6} case (panel (*b*)), the solution transitions from the *M* vertex and approaches the *ϕ*=10^{3} case (panel (*c*)), the solution transitions from the *M* vertex to the

## 7. Summary

This paper presents a closed-form approximate solution for a penny-shaped hydraulic fracture, whose behaviour is governed by a simultaneous interplay of fluid viscosity, fracture toughness and fluid leak-off into the formation. This development is made possible by combining the approximation for the fracture width profile, which uses the multiscale tip asymptotic solution in combination with the global fluid volume balance. First, the limiting regimes of the approximation are analyzed and the results are compared with the existing solutions. Then, a solution map is presented. This solution map indicates the areas of applicability of the limiting solutions and allows one to easily determine whether a problem with given parameters corresponds to one of the limiting cases or not. In addition, the developed solution captures the scaling that is associated with the transition from one limiting solution to another. In order to estimate accuracy of the approximate solution, a reference numerical solution is constructed. The latter numerical solution uses a moving mesh together with the implicit time stepping and a multiscale tip asymptotic solution to locate the moving fracture front. Results indicate that predictions of the fracture width and radius by the developed approximation are accurate to within a fraction of a per cent. The fluid pressure results are less accurate and lie within 2% error if measured half-way to the fracture tip. The importance of the developed approximate solution is twofold. First, it can be used to estimate fracture radius very quickly, in which case it can be used for a rapid fracture design calculations. Second, it can be used as a reference solution to test accuracy of more complex hydraulic fracture simulators (such as non-axisymmetric or non-planar), and can also be used as an initial condition for the aforementioned numerical codes.

## Data accessibility

Matlab code that produces the closed-form solution for a radial hydraulic fracture that is described in this paper is available from the Dryad Digital Repository at http://dx.doi.org/10.5061/dryad.gh469 [47].

## Competing interests

I have no competing interests.

## Funding

Start-up funds provided by the University of Houston are greatly acknowledged.

## Appendix A. Functions g δ ( K ^ , C ^ ) and Δ ( K ^ , C ^ )

This appendix provides expressions for the functions

With reference to the scaling (3.7) and the fact that *f* can be introduced as
*β*_{m}=2^{1/3}3^{5/6}. As mentioned in [34], the solution varies spatially as *δ* is given by
*δ*-corrected solution (3.7) can be written as

## Appendix B. Numerical scheme

To construct the numerical scheme, equation (2.2) is rewritten using *ρ*=*r*/*R*(*t*) and scaling (5.1)–(5.3) as

The spatial coordinate *ρ* is discretized as *j*=1…*N*, in which case *τ* is discretized uniformly on a logarithmic scale. Piecewise constant approximation for *Ω* is used, in which case *Ω*^{i} represents an array of values of *j*. In this situation, the elasticity equation (B2) is discretized as
*j*±1/2 are calculated as an average between the corresponding values of these quantities, and the source/leak-off term is
*δ*_{1j} is the Kronecker delta. The values of [*B**Ω*^{i}]_{N} and [** A**(

*Ω*^{i})

*Π*^{i}]

_{N}are obtained by integration of (B1) over the last element and using the no-flux condition at the tip.

In order to capture multiscale behaviour near the fracture tip, the following propagation condition is used
*Ω*_{a} is the scaled tip asymptotic solution. The latter equation implies that the numerical solution follows the asymptotic solution from the penultimate element to the tip, and allows one to determine the propagation velocity *V* . To successfully use this condition, pressure at the tip element *Ω*_{s}∝*d*^{(1+δ)/2} and *τ*−*τ*_{0}=*γ*(1−*ρ*)/*V* , in which case

The numerical scheme consists of solving (B3)–(B7) for *j*=1…*N*−1) and *γ*^{i}=*γ*^{i−1}+*V* *Δτ*.

- Received September 22, 2016.
- Accepted November 7, 2016.

- © 2016 The Authors.

Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.