An approximate solution for a penny-shaped hydraulic fracture that accounts for fracture toughness, fluid viscosity and leak-off

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.


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][6][7][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][13][14], multiple hydraulic fractures [15][16][17] or fracture networks [18]. Detailed reviews of various hydraulic fracturing models can be found in [19][20][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][29][30][31][32][33][34][35]. On the other hand, time evolution and regimes of propagation were studied for a plane-strain KGD fracture in [36][37][38][39][40], and for a penny-shaped hydraulic fracture in [40][41][42][43][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][41][42][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. parameters that appear in the model, which can be introduced in the scaled form, for convenience, as where μ 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 where 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] where R is the fracture radius, and the kernel is The functions 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] w → K E (R − r) 1/2 , r → R, (2.5) which implies that the stress intensity factor is equal to the fracture toughness for a propagating fracture. The propagation condition (2.5) should also be complemented by a zero flux condition at the fracture tip, i.e. 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 where the relations q(R, t) = 0 and w(R, t) = 0 were used to obtain the result.
3. Approximate solution for a penny-shaped hydraulic fracture

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 where 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Ṙ(t), whereas λ is the parameter that is determined later. Here, the tip asymptotic solution, 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 ( 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 w a (s) ∝ sδ, where s = R − r is the distance from a point inside the fracture to the tip andδ is a slowly varying function, which is also part of the solution (i.e. known). This fact allows to reduce (3.1) to Motivated by the results of the limiting asymptotic solutions, for which R(t) ∝ t α and α is a number that is equal to either 1/4, 2/5, or 4/9 [27,44], 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 where ρ = r/R is the scaled spatial coordinate. The integrals in equation (3.3) can be evaluated, in which case (3.3) reduces to where B(a, b) denotes the beta function, and where 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Ṙ = αR/t. Also note that the expression forδ comes from the tip asymptotic solution. Then, by specifying λ and taking an initial guess for α = 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 α = d log(R)/d log(t). After that, equation (3.4) is solved again. The iteration procedure is performed until the desired level of convergence is reached, which is typically achieved quickly owing to relatively modest dependence on α. 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 To calculate fluid pressure, one can substitute (3.2) into (2.3) and evaluate the integral, so that where the function F (ρ, λ,δ) can be evaluated numerically, and the kernel 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.

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 where g δ is a known function (see appendix A), andṘ = αR/t is used for velocity of the crack tip. By introducing the following scales and the dimensionless quantitieŝ the global fluid balance equation (3.4) can be rewritten as the parameter λ(K,Ĉ, α) is specified later, andδ = 1 2 (1 + (K,Ĉ)) (see appendix A). The first equation in (3.7), the first relation in (3.9) and (3.10) can be combined to yield the system of equations To obtain the solution, the system of equations (3.12) is solved numerically forK(t) andĈ(t) using Newton's method. Note thatQ 0 is a parameter, and the initial value of α is taken as 4/9, which is then updated iteratively using α = d log(R)/d log(t) (so that the system (3.12) is solved a few times until α converges). Once the time historiesK(t) andĈ(t) are obtained, equation (3.7) is used to findŝ(t), and the relations in (3.9) are used to calculateR(t) andŵ a (t). The values ofδ and λ (see §4.5) are also expressed in terms ofK(t) andĈ(t). The pressure and the fracture opening solutions can be obtained from (3.6) and (3.2), respectively, asp whereas the efficiency can be expressed using (3.5), (3.10) and (3.11) as Finally, the solution is given in terms of the fracture radiusR(t), fracture widthŵ(ρ,t), fluid pressurê p(ρ,t) and efficiency η(t). Equations (3.8) and (3.9) can be used to convert the solution to the unscaled form if needed.

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) [27,44] negligible. At the same time, the efficiency η = 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 and η = 0 corresponds toĈ → ∞, as can be seen from (3.11) and (3.14).
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 (orK andĈ), to zero makes the scaling (3.8) and (3.9) deficient. It is still possible, however, to use a combination of variables, such asR/ŵ 2 a , because it may not depend onĈ and, consequently, on 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).

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 Because the M vertex solution corresponds to the case of no toughness and no leak-off, it is possible to setK =Ĉ = 0 in (4.2), in which case α = 4/9. Note that B(0, 0, 4/9) = 0.1979 and g δ (0, 0) = 1/β 3 m = 0.0321 (see appendix A) are constants, where β m = 2 1/3 3 5/6 . By introducing scaling that is associated with the M vertex solution [13,41,44] where α = 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 where the fact thatδ = 1 2 (1 + (0, 0)) = 2 3 is used. Comparison of the fracture width at the wellbore Ω 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 which can be written in the dimensional form using (4.3) as where the function F is defined in (3.6).

4.2.M vertex limit of the solution
In order to obtain theM limit, one should consider the limiting case of no fracture toughness, i.e.K = 0, and zero efficiency, i.e. η = 0. In this situation, equation (3.14) implies that 3 2 ). (4.8) Because the zero efficiency case corresponds to large values ofĈ, one can use the asymptotic behaviour of g δ at largeĈ, result, equations (3.9) and (3.12) can be rewritten in the form The combination of (4.8) and (4.9) finally leads tô 3/8t 1/16 , (4.10) in which case it is clear that α = d log(Rm)/d log(t) = 1/4, and hence the beta function can be evaluated as B( 1 2 , 3 2 ) = π/2. By introducing scaling that is associated with theM vertex solution [13,44] relations (4.10) can be reduced to and solutions for the fracture width and pressure variations can be written using (3.13) as whereδ = 5/8 is used, which corresponds to the caseK = 0 andĈ → ∞. Comparison of the fracture width at the wellbore Ωm(0) and the fracture radius γm in (4.13) with the reference solution for theM vertex [13,44] indicates that λm = 0.397 provides the most accurate approximation, in which the error of Ωm(0) is below 0.1%, whereas the value of γm is captured precisely. TheM vertex limit of the approximate solution can be summarized from (4.12) and (4.13) as Ωm(ρ) = 1.0574(1 + ρ) 0.397 (1 − ρ) 5/8 , Πm(ρ) = 3.0931F (ρ, 0.397, 5 8 ), γm = 0.4502, (4.14) and can be written in the dimensional form using (4.11) as and where the function F is defined in (3.6).
Given that F (ρ, 1 2 , 1 2 ) = π/ 8 √ 2 and B(1, 0, α) = 3 √ 2 −1 , the solution can be summarized as which coincides with the K vertex solution in [41,44]. The above results can be written in the unscaled form as Note that because the fracture width for the 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
TheK vertex limit of the solution corresponds to the situation for whichK = 1 and η = 0. As a result, the efficiency relation (3.14) reduces to As for theM solution, α = d log(Rk)/d log(t) = 1/4, and the beta function can be evaluated as B( 1 2 , 3 2 ) = π/2. By introducing scaling that is associated with theK vertex solution [13,42,44] relations (4.23) can be reduced to 25) and the width and pressure spatial variations can be obtained using (3.13) as where, similarly to the K vertex limit,δ = 1 2 was used. As for the K vertex case, the fracture width for thẽ K vertex solution is an ellipse [42,44], and λk = 1/2 leads to the exact solution. which coincides with theK vertex solution in [42,44]. The unscaled form of (4.27) is and Similar to the K vertex solution, theK limit of the approximate solution coincides with theK limit of the exact solution, because the approximation (3.1) captures elliptic width profile precisely.

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, λm = 0.397, λ k = 0.5 and λk = 0.5. Owing to a small variation of λ, one may consider the following interpolation procedure: where η 0 (K,Ĉ, α) is the approximation for the efficiency (3.14) for which λ = 0.5 is used. The interpolation (4.29) captures the correct values of λ at all four vertices, and provides an approximation elsewhere. It is important to note that the choice ofK 4 in lieu of e.g.K is due to the fact that the function g δ depends onK 4 whenĈ → ∞ (transition betweenM andK regimes); see appendix A. Note that the transition from M to K (η = 0) is less important for the interpolation purposes, because λ k ≈ λ m and any interpolation would work. The transitions from M toM and from K toK are determined by the efficiency, hence the efficiency is used as a parameter for the interpolation in (4.29). Because even a constant approximation λ ≡ 0.5 is able to provide relatively accurate approximation owing to a small variation of λ, the accuracy of the interpolation (4.29) is expected to be sufficient. This statement will be verified later in §6, where the approximate solution is compared with the reference numerical solution.

Structure of the solution
Solution in the original scaling is given by the quantitiesR,ŵ,p and η, which depend on the dimensionless timet and the parameterQ. At this point, it is useful to rewrite all the variables in terms of the scaling that will be used in the reference numerical solution in §6. In particular, the scaling of the numerical solution is determined by the following quantities [44] ε = μ E t mk whereas the leak-off parameter φ is defined as In this setting, the relationship between the new scaling and the original scaling (3.8) and (3.9) can be summarized as where i corresponds to either M, K,M orK vertex solutions. Note that all quantities in formula (5.4) are written in terms of the scaling (5.3). 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 theK vertex. Along its way, the fracture may pass through the K vertex, through theM vertex, or be in the transition region between the two, as was described in [27,44]. Black lines in figure 1 show the boundaries of applicability of the vertex solutions and are determined next. Note that all the timescales that are associated with the transitions from one vertex solution to another (obtained in the remainder of this section using limits of the approximate solution) are consistent with the results in [44], which were obtained using a different method. This indicates that the developed approximate solution precisely captures the scaling that is associated with the transitions from one vertex solution to another.
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,Ĉ = 0 and equations (3.12) and (3.14) reduce tô Solution of equation ( where τ The associated boundaries of the MM transition, shown in figure 1, are calculated numerically as which determines the KK edge solution. The solution depends on the parametert/Q 2/3 , which can be rewritten using (5.2) and (5.3) as τ φ 5/6 . In this case, the dimensionless time that quantifies the KK transition is The associated boundaries of the KK transition, shown in figure 1, are calculated numerically as The associated boundaries of theMK transition, shown in figure 1, are calculated numerically as corresponds to the end of the transition (near theK vertex).

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 andK vertices, because the pressure is nearly constant, in which case the fluid flow within the fracture becomes irrelevant. At the same time, because the approximate width solution is exact at the K andK vertices (and also KK transition), the approximation is expected to be very accurate in these regions. In addition, the numerical scheme becomes less computationally efficient (i.e. requires finer mesh for convergence) near theM vertex solution. This occurs, because the efficiency becomes very small, in which case most of the injected fluid is leaks into the formation. The accuracy of the approximation at theM vertex, however, was evaluated earlier and there is no need to check it again with the numerical solution. 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 where the pressure is not evaluated at the origin owing to a singularity located at that point. All the quantities are normalized by the M vertex solution (4.6), which can be written in terms of the scaling (5.1) and (5.3) as The purpose of such normalization is to achieve better visual separation between the curves in figure 2.
The vertex solutions are shown by the blue (M vertex), red (K vertex), green (M vertex) and magenta (K vertex) lines. All the solutions originate near the M vertex (i.e. the blue lines), and then transition to either K,K,M or intermediate solution. Note that, eventually, all solutions (except the φ = 0 case) reach theK limiting case, as can be clearly seen from figure 1.
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 where · indicates average over all time points, subscripts 'apr' and 'num' correspond to the approximation and numerical solution respectively, whereas A is either radius γ (τ ), efficiency η(τ ), wellbore width Ω 0 = Ω(0, τ ) or pressure Π 0 = Π (0.5, τ ). 2) for the fracture radius γ (τ ), efficiency η(τ ), wellbore width Ω 0 = Ω(0, τ ) and pressure Π 0 = Π (0.5, τ ). (b) Variation of the maximum error (maximum is calculated over φ) versus number of mesh points N of the numerical solution for the fracture length γ , efficiency η, wellbore width Ω 0 and pressure Π 0 . Figure 3a 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 3b 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 3a and its analogues for different meshes. Results in figure 3b 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 (M vertex) and magenta (K vertex) lines for the first and the last value of τ . 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 theK vertex. For the φ = 10 3 case (panel (c)), the solution transitions from the M vertex to theM vertex. These observations agree with the results shown in figure 1. Fracture width variations are consistently accurate, whereas the fluid pressure estimations deviate notably in the tip region. This is the weakest point of the approximate solution, and, therefore, it should not be used to evaluate pressure for ρ 0.8. There might be a possibility to overcome this error by using (2.2) to estimate the pressure field. This, however, is a minor correction to the developed procedure and therefore is beyond the scope of this study.

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 leakoff 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]. In order to capture multiscale behaviour near the fracture tip, the following propagation condition is used where Ω 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 Π i N is treated as an unknown. Because the tip asymptotic solution satisfies (3.7), then where d = 3 2 γ ρ signifies the distance from the centre of the penultimate element to the tip. Equation (B 5) is solved forŝ using Newton's method, and the velocity of propagation is calculated. Because the solution in the tip element follows the asymptotic solution, then it is possible to determine its average opening as where it is used that Ω s ∝ d (1+δ)/2 and δ = (K,Ĉ). The factor 2/3 ensures that the calculation is performed only within the tip element (and does not continue the middle of the penultimate element). Leak-off in the tip element is calculated by taking τ − τ 0 = γ (1 − ρ)/V, in which case The numerical scheme consists of solving (B 3)-(B 7) for Ω i j (for j = 1 . . . N − 1) and Π i N , which is done iteratively for each time step. The fracture radius is then updated as γ i = γ i−1 + V τ .