Fluid velocity based simulation of hydraulic fracture: a penny shaped model—part I: the numerical algorithm

In the first part of this paper, a universal fluid velocity based algorithm for simulating hydraulic fracture with leak-off, previously demonstrated for the PKN and KGD models, is extended to obtain solutions for a penny-shaped crack. The numerical scheme is capable of dealing with both the viscosity and toughness dominated regimes, with the fracture being driven by a power-law fluid. The computational approach utilizes two dependent variables; the fracture aperture and the reduced fluid velocity. The latter allows for the application of a local condition of the Stefan type (the speed equation) to trace the fracture front. The obtained numerical solutions are carefully tested using various methods, and are shown to achieve a high level of accuracy. Electronic supplementary material The online version of this article (10.1007/s11012-018-0899-y) contains supplementary material, which is available to authorized users.


Introduction
Hydraulic fracture (HF) is the phenomenon of a fluid driven crack propagating in a solid material. It can be encountered in various natural processes, such as subglacial drainage of water or during the extension of magmatic intrusions in the earth's crust. Simultaneously the underlying physical mechanism is very important in numerous man-made activities. Hydrofracturing can appear as an unwanted and detrimental factor in underground CO 2 or waste repositories [1]. On the other hand, intentionally induced hydraulic fractures constitute the essence of fracking technology -a method used when stimulating unconventional hydrocarbon reservoirs [2] or for geothermal energy exploitation [3]. All of these applications create demand for a proper understanding and prediction of the process of hydraulic fracture.
As a result of the multiphysical nature of the underlying physical phenomenon and complex interactions between the component physical fields, the mathematical modeling of hydraulic fractures represents a significant challenge. The main difficulties arise due to: (1) strong non-linearities resulting from interaction between the solid and fluid phases, (2) singularities in the physical fields, (3) moving boundaries, (4) degeneration of the governing equations at the crack tip, (5) leak-off to the rock formation, (6) pronounced multiscaling effects, vii) complex geometry.
The first theoretical models of hydraulic fracture were created in 1950s (see for example [4] and [5]).

Electronic supplementary material
The online version of this article (https://doi.org/10.1007/s11012-018-0899-y) contains supplementary material, which is available to authorized users. Subsequent research led to the formulation of the socalled classic 1D models: PKN [6,7], KGD (plane strain) [8,9] and penny-shaped/radial [10]. Up to the 1980s these very simplified models were used to design and optimize the treatments used in HF. The increasing number and size of fracking installations, alongside the simultaneous advance in computational techniques, resulted in the formulation of more sophisticated and realistic models of HFs. A comprehensive review of the topic can be found in [11].
Though superseded in most practical applications, the classic 1D models remain a significant avenue of research into the fundamentals of HF. They enable one to investigate some inherent features of the underlying physical process, the mathematical structure of the solution, and finally to construct and validate computational algorithms. Substantial advances have been achieved in this area throughout the last 30 years by way of a cyclical revision of these classic formulations. It was not until 1985 [12] that the importance of the solution tip asymptotics was first noticed, specifically for the KGD and penny shaped cracks. The explicit form of the tip asymptotic solution for the PKN model was given in 1990 by Kemp [13]. Moreover, in this publication the author remarked, for the first time, that when properly posed the Nordgren's model constitutes a Stefan-type problem and as such needs an additional boundary condition which equates the crack propagation speed with the velocity of the fluid front. However, this important idea was abandoned for the next 20 years until being rediscovered by Linkov [14] in 2011. The author proved that the general HF problem is ill-posed and proposed a regularization technique based on application of the aforementioned Stefan conditioncalled there the speed equation. The numerous investigations carried out since the beginning of the present century for the KGD [15][16][17][18] and penny-shaped models [19][20][21] have led to the importance of the problem's multiscale character being recognized. It is now well understood that the global response of the fluid driven fracture is critically dependent on the interaction between competing physical processes at various temporal and spatial scales. Depending on the intensity of various energy dissipation mechanisms, as well as the fracturing fluid and solid material properties, the hydraulic fracture evolves in the parametric space encompassed by the limiting regimes: (1) viscosity dominated, (2) toughness dominated, (3) storage dominated, (4) leak-off dominated.
Bearing in mind the whole complexity of the problem, it still remains an extremely challenging task to deliver credible solutions which reflect all of the desired features. The relative simplicity of the classic 1D models means that they are well suited to the task of creating benchmarks, used when developing and verifying more advanced solutions and algorithms. For the KGD and PKN models one can find in the literature a number of credible results, including recently developed simple and accurate approximate solutions, that can be utilized for the aforementioned purposes [22][23][24][25]. A more complete review of recent benchmarks will be provided in part II of this paper.
The aim of the first paper is to meet the demand for benchmark solutions to the radial HF model and: (1) to deliver a dedicated computational scheme capable of obtaining highly accurate numerical solutions, (2) introduce purely analytical solutions to the problem obtained for a predefined non-zero leak-off function, (3) introduce and verify an alternative measure of the computational error, to use when no analytical solution is available.
To this end the self-similar formulation of the penny-shaped model will be analyzed. The numerical computations will be performed according to a modified form of the universal algorithm introduced in [24,25]. It employs a mechanism of fracture front tracing, based on the speed equation approach [23], coupled with an extensive use of information on the crack tip asymptotics and regularization of the Tikhonov type (the technical details of both concepts can be found in [26,27]). The modular architecture of the computational scheme facilitates its adaptation to the problem of radial HF.
It is worth stating that the second part of this paper will introduce simple to use semi-analytical approximations of numerical benchmark solutions obtained for the case of an impermeable solid, and comparisons with other available data will be performed. Both parts are written in such a way that they can be read as individual, independent papers (for a unified version of the text, see arXiv:1612.03307).
The paper is organized as follows. The basic system of equations describing the problem is given in Sect. 2. Next, normalization to the dimensionless form is carried out. In Sect. 3, comprehensive information about the solution asymptotics is presented, which is heavily utilized in the subsequent numerical implementation. New computational variables, the reduced fluid velocity and modified fluid pressure derivative, are introduced. The advantages of both are outlined, and the problem is reformulated in terms of the new variables. In Sect. 4 the governing system of equations is reduced to the time independent selfsimilar form. This formulation is used in Sect. 5 to construct the computational algorithm. The accuracy and efficiency of computations are examined against newly introduced analytical benchmark examples. Alternative error measures are proposed for the cases where no closed-form analytical solution is available. Section 6 contains the final discussion and conclusions. Some additional information concerning the limiting cases of Newtonian and perfectly plastic fluids is collected in the appendices.

Problem formulation
Let us consider a 3D penny-shaped crack, defined in polar coordinates by the system fr; h; zg, with associated crack dimensions flðtÞ; wðtÞg as the fracture radius and aperture respectively, noting that both are functions of time. The crack is driven by a point source located at the origin, which has a known pumping rate: Q 0 ðtÞ. The fluid's rheological properties are described by the power-law model [28], as is common in the literature. The rational behind this choice is outlined in [25]. We have that, as the flow is axisymmetric, all variables will be independent of the angle h and it is sufficient to analyse the problem for only r ! 0.
The fluid mass balance equation is as follows: where q l ðr; tÞ is the fluid leak-off function, representing the volumetric fluid loss to the rock formation in the direction perpendicular to the crack surface per unit length of the fracture. We will assume it to be a predefined and smooth function which is bounded at the fracture tip. Meanwhile, q(r, t) is the fluid flow rate inside the crack, given by the Poiseuille law: with p(r, t) being the net fluid pressure on the fracture walls (i.e. p ¼ p f À r 0 , where p f is the total pressure and r 0 is the confining stress). The constant M 0 is the so-called modified fluid consistency index M 0 ¼ 2 nþ1 ð2n þ 1Þ n =n n M, where M stands for the consistency index (relating the shear stress and strain rate) and 0 n 1 is the fluid behaviour index. The non-local relationships between the fracture aperture and the pressure (elasticity equations) are as follows: pðr; tÞ ¼ E 0 lðtÞ A½wðr; tÞ; wðr; tÞ ¼ lðtÞ where E 0 ¼ Y=ð1 À m 2 Þ, with Y being the Young's modulus and m the Poisson ratio. The operator A and its inverse take the form: ð5Þ with the pertinent kernels being: Gðn; sÞ ¼ K, E are the complete elliptic integrals of the first and second kinds respectively, and F is the incomplete elliptic integral of the first kind, given in [29]. These equations are supplemented by the boundary condition at r ¼ 0, which defines the intensity of the fluid source, Q 0 : the tip boundary conditions: wðlðtÞ; tÞ ¼ 0; qðlðtÞ; tÞ ¼ 0; ð9Þ and appropriate initial conditions describing the starting crack opening and length: Additionally, it is assumed that the crack is in continuous mobile equilibrium, and as such the classical crack propagation criterion of linear elastic fracture mechanics is imposed [30]: where K Ic is the material toughness while K I is the stress intensity factor. The latter is computed according to the following formula [31]: Throughout this paper we accept the convention that when K Ic ¼ 0 the hydraulic fracture propagates in the viscosity dominated regime. Otherwise the crack evolves in the toughness dominated mode. Each of these two regimes is associated with qualitatively different tip asymptotics, which constitutes a singular perturbation problem as K Ic ! 0, and leads to serious computational difficulties in the small toughness range. Finally, noting (1) and (8), the global fluid balance equation is given by: Z lðtÞ 0 r wðr; tÞ À w 0 ðrÞ ½ dr þ Z t 0 Z lðtÞ 0 rq l ðr; sÞ dr ds The above set of equations and conditions represents the typically considered formulation for a pennyshaped hydraulic fracture [19]. In order to facilitate the analysis we shall utilize an additional dependent variable, v, which describes the average speed of fluid flow through the fracture crosssection [23]. It will be referenced to in the text as the fluid velocity (often referred to in the literature as the particle velocity, e.g. [24,25]), and is defined as: vðr; tÞ ¼ qðr; tÞ wðr; tÞ ; v n ðr; tÞ ¼ À We assume that the leak-off q l is such that the fluid velocity is finite at the crack tip, meaning that v has the following property: Additionally, given that the fracture apex coincides with the fluid front (there is no lag), and that the fluid leak-off at the fracture tip is bounded, the so-called speed equation [14] holds: This Stefan-type boundary condition constitutes an explicit method, as opposed to an implicit level-set method [32], and can be effectively used to construct a mechanism of fracture front tracing. The advantages of implementing such a condition have been shown in [24,25].

Problem normalization
In order to make the presentation clearer, we will assume that 0\n\1 in the main body of the text. Any modification to the governing equations and numerical scheme in the limiting cases n ¼ 0 and n ¼ 1 are detailed in ''Appendix 1''. We normalize the problem by introducing the following dimensionless variables: wherer 2 0; 1 ½ and l Ã is chosen for convenience. We note that such a normalization scheme has previously been used in [24,25], and that it is not attributed to any particular influx regime or asymptotic behaviour of the solution.
While this form of the elasticity operator has not previously been used in the case of a penny-shaped fracture, an analogous form of the elasticity equation for the KGD model has been utilized in [24,25], where its advantages in numerical computations have been demonstrated. Notably, the kernel function K exhibits better behaviour than the weakly singular kernel G (7), having no singularities for any combination ofr; y f g. Additionally, Eq. (19) can be easily transformed to obtain p 0 and then substituted into (26), meaning that the latter can be utilized without the additional step of deriving the pressure function needed for the classic form of the operator.
Next the boundary conditions (9), in view of (15), transform to a single condition: alongside the initial conditions (10): The source strength (8) is now defined as: It is worth restating that there are minor differences to both the asymptotics and fundamental equations in the limiting cases n ¼ 0 and n ¼ 1. These are explained in further detail in ''Appendix 1''.

Crack tip asymptotics, the speed equation and proper variables
A universal algorithm for numerically simulating hydraulic fractures has recently been introduced in [24,25] and tested against the PKN and KGD (plane strain) models for Newtonian and shear-thinning fluids. It proved to be extremely efficient and accurate. Its modular architecture enables one to adapt it to other HF models by simple replacement or adjustment of the basic blocks. In the following we will construct a computational scheme for the radial fracture based on the universal algorithm. To this end we need to introduce appropriate computational variables, and to define the basic asymptotic interrelations between them. For the sake of completeness detailed information on the solutions tip asymptotic behaviour, for different regimes of crack propagation, are presented below.

Crack tip asymptotics
In the viscosity dominated regime the crack tip asymptotics of the aperture and pressure derivative can be expressed as follows:w ðr;tÞ ¼w 0 ðtÞ 1 Àr 2 À Á a 0 þw 1 ðtÞ 1 Àr 2 À Á a 1 op or ðr;tÞ ¼p 0 ðtÞ 1 Àr 2 À Á a 0 À2 þp 1 ðtÞ 1 The crack tip asymptotics of the pressure function can be derived from the above, however this form is given due to its use in computations (this will be explained in further detail later). As a consequence the fluid velocity behaves as: vðr;tÞ ¼ṽ 0 ðtÞ þṽ 1t ð Þ 1 Àr Note that we requireṽ 0 ðtÞ [ 0 to ensure the fracture is moving forward. The values of constants a i , b i are given in Table 1. The general formulae for the limiting cases n ¼ 0 and n ¼ 1 remain the same as (35)-(37), with the respective powers a i , b i again being determined according to Table 1. Now, let us adopt the following notation for the crack propagation speed, based on the speed equation (20) and the tip asymptotics (37): Here LðwÞ [ 0 is a known functional and C is a positive constant. In the viscosity dominated regime we have that: Additionally, we can directly integrate (38) in order to obtain an expression for the fracture length: 3.1.2 Toughness dominated regime (K Ic [ 0) Near the fracture front the forms of the aperture and fluid velocity asymptotics remain the same as in the viscosity dominated regime (35), (37). Meanwhile the pressure derivative asymptotics yields: The values of a i , b i for this regime are provided in Table 1. The limiting cases n ¼ 0 and n ¼ 1 are discussed in ''Appendix 1'' (Eqs. 78, 87 respectively). We again use notation (38) for the crack propagation speed, however the values of the functional L and the C will in this case be: while the fracture length will be given by (40).

Reformulation in terms of computational variables
It is readily apparent that the choice of computational variables plays a decisive role in ensuring the accuracy and efficiency of the computational algorithm [23,24,26]. Let us introduce a new system of proper variables which are conducive to robust numerical computing.
• The reduced fluid velocity Uðr;tÞ: Uðr;tÞ ¼rṽðr;tÞ Àr 2ṽ 0 ðtÞ: It is a smooth, well behaved and non-singular variable that facilitates the numerical computations immensely. It is bounded at the crack tip and the fracture origin. The advantages of using an analogous variable in the PKN and KGD models have previously been demonstrated in [24,25].
• The modified fluid pressure derivative Xðr;tÞ: r n Xðr;tÞ ¼r n op or À X 0 ðtÞ; ð44Þ It reflects the singular tip behavior ofp 0r , having the same tip asymptotics as the pressure derivative, however it is bounded at the fracture origin. From (44) the pressure can be immediately reconstructed as: where the term C p follows from (25): This auxiliary variable will primarily be used in numerical computation of the elasticity operator.
The following interrelationship exists between the newly introduced variables: Since under this new scheme U is bounded at the fracture tip, the source strength (31) and the boundary condition (29) can now be expressed as: wð0;tÞUð0;tÞ ¼Q 0 ðtÞ 2p ;wð1;tÞ ¼ 0: By utilizing the boundary condition (49) 1 , the relationship between the new variables (48) can be represented in the form: Note that this is not only a more concise representation of (48) but it also does not depend on LðtÞ, which will be beneficial when analyzing the self-similar formulation. In this way the computational scheme will be based on: the crack opening,w, the reduced fluid velocity, U, and an auxiliary variable, the modified fluid pressure, X. By substituting the new variable U from (43) into the continuity Eq. (18), we obtain the modified governing equation: Additionally, the elasticity Eq. (26) can be now expressed as follows: where K is given in (27), while G n is defined by: It can be easily shown that this function is well behaved in the limits.

Self-similar formulation
In this section we will reduce the problem to its timeindependent self-similar version. For this formulation we will define the computational scheme used later on in the numerical analysis.
We begin by assuming that a solution to the problem can be expressed in the following manner:w where WðtÞ is a smooth continuous function of time. Such a separation of variables enables one to reduce the problem to a time-independent formulation when W is described by an exponential or a power-law type function. From here on the spatial components will be marked by a 'hat'-symbol, and will describe the selfsimilar quantities. It is worth noting that the separation of spatial and temporal components given in (54) ensures that the qualitative behaviour of the solution tip asymptotics remains the same as in the timedependent variant.

The self-similar representation of the problem
We wish to examine two variants of the time dependent function: In both cases the fluid leak-off function will be assumed to take the form: The self-similar reduced fluid velocity (43), modified fluid pressure derivative (44), (45) and pressure (46) are defined by: UðrÞ ¼rvðrÞ Àr 2v 0 ;rXðrÞ ¼r dp dr ÀX 0 ; ð57Þ with It is immediately apparent from (38) and (54) that the self-similar crack propagation speed is given by: Àŵ nþ1 dp dr Note again that the qualitative asymptotic behaviour of the aperture, pressure and fluid velocity asr ! 0 andr ! 1 remains the same as in the time dependent version of the problem (35), (36), (37), (41). In the self-similar formulation, the multipliers of respective terms are time-independent. The self-similar counterparts of the elasticity Eqs. (22) and (23)  As the integral and function G n ðrÞ both tend to zero faster than the square root term at the fracture tip, it immediately follows that, in the toughness dominated case (K Ic [ 0), the leading asymptotic term of the aperture (35) is given by: The self-similar fluid velocity takes the form: vðrÞ ¼ Àŵ nþ1 ðrÞ dpðrÞ dr The governing Eq. (51) becomes: with the value of q in each case, alongside the fracture length, provided in Table 2. Meanwhile the fluid balance condition (21) becomes: The self-similar stress intensity factor (25) is given by: Finally, the system's boundary conditions (49) transform to: In the general case with 0\n\1 these equations represent the full self-similar problem. Some modifications are necessary in the special cases when n ¼ 0 and n ¼ 1. These differences are outlined in ''Appendix 1''.

Numerical results
In this section we will construct an iterative computational scheme for numerically simulating hydraulic fracture. The approach is an extension of the universal algorithm introduced in [24,25]. The computations are divided between two basic blocks, the first of which utilizes the continuity equation and the latter using the elasticity operator. The previously introduced computational variables, alongside the known

Computational scheme
The algorithm is constructed using the approach framework introduced for the PKN and KGD models in [24,25]. The numerical scheme is realized as follows: 1. An initial approximation of the apertureŵ ¼ŵ jÀ1 is taken, such that it has the correct asymptotic behaviour and satisfies the boundary conditions. 2. The fluid balance Eq. (68) is utilized to obtain the asymptotic term(s)ŵ j 0;1 needed to compute the fluid velocityv j 0 using (61).

Having the above values the reduced fluid velocitŷ
U j is reconstructed by direct integration of (67). Tikhonov type regularization is employed at this stage. 4. Equations (57) and (66) is then used to obtain an approximation of the modified fluid pressure derivativeX, and the elasticity Eq. (64) serves to compute the next approximation of the fracture apertureŵ j . 5. The system is iterated until all variablesÛ,ŵ and v 0 converge to within prescribed tolerances.
We will demonstrate in this section that this scheme, combined with an appropriate meshing strategy, yields a highly accurate algorithm. A more detailed description of the algorithm's construction has been outlined in [24,25]. When iterated, the system of discretized equations converges to a final solution in under 20 iterations, with computation times on a standard laptop under 5 s when taking N ¼ 20, and under 30 s when N ¼ 300.
It is worth noting that, due to the degeneration of the Poiseuille equation when n ¼ 0, it can no longer be used to compute the fluid flow rate or the fluid velocity. However, thanks to the modular structure of the proposed algorithm, one can easily adapt it to this variant of the problem. In this case a special form of the elasticity Eq. (96) is utilized to obtain the aperture, with the fluid velocity being reconstructed using relations (97) and (98).

Accuracy of computations
In this subsection we will investigate the accuracy of computations delivered by the proposed numerical scheme. To this end a newly introduced set of analytical benchmark solutions with a non-zero fluid leak-off function will be used. Alternative measures for testing the numerical accuracy in the absence of exact solutions will then be proposed and analysed. Next, the problem of a penny-shaped hydraulic fracture propagating in an impermeable material will be considered. The accuracy of numerical solutions will be verified by the aforementioned alternative measures. Simple, semi-analytical approximations, which mimic the obtained numerical data to a prescribed level of accuracy will be provided in the second part of this paper, alongside a comparative analysis with other solutions available in the literature.

Analysis of computational errors against analytical benchmarks
The first method of testing the computational accuracy is by comparison with analytical benchmark solutions. Respective closed-form benchmarks with predefined, non-zero, leak-off functions are outlined in the supplementary material associated with this paper. They have been constructed for both the viscosity and toughness dominated regimes, for a class of shearthinning and Newtonian fluids. All of the analytical benchmarks used for comparison are designed to ensure physically realistic behaviour of the solution while maintaining the proper asymptotic behaviour. In all numerical simulations the power-law form of the time dependent function W 2 (55) 2 is used to ensure consistency with results in the literature (e.g. [19]). The accuracy of computations is depicted in Figs. 1  and 2, for varying number of nodal points N. A nonuniform spatial mesh was used, with meshing density increased near the ends of the interval (the same type of mesh was used for all n). The measures dw, dv, describing the average relative error of the crack opening and fluid velocity, are taken to be: whereŵ Ã andv Ã denote the exact solutions forŵ andv.
The results clearly show that the values of both error measures decrease monotonically with growing N. For a fixed number of nodal points N, dw is lower than dv, but within the same order of magnitude. One can observe a sensitivity of the results to the value of the fluid behaviour index n. Here, the level of error measures can vary up to one order for a constant N. This trend can be alleviated by adjusting the mesh density distribution to the value of n (i.e. to the varying asymptotics of solution), however such an investigation goes beyond the scope of this paper. In general, it takes fewer than N ¼ 300 nodal points to achieve the relative errors of the level 10 À7 .
In cases when the exact solution is not prescribed an alternative method of testing the solution accuracy is required. The method outlined here relies on the fact that the solution converges to the exact value at a known rate, with respect to the number of nodal points. The convergence rate has been established numerically to behave as 1=N 3 . As a result the following estimation holds: whereĝ 1 ðrÞ ¼ŵðrÞ andĝ 2 ðrÞ ¼vðrÞ. A i and B i are constants which can be found numerically through use of a least-squares approximation. Next, one can define the limiting value of (72) as: forĝ Ã 1 ðrÞ ¼ŵ Ã ðrÞ,ĝ Ã 2 ðrÞ ¼v Ã ðrÞ. Knowing this, the following alternative error measures can be proposed: Using this strategy, it is possible to identify the relative rate at which the solution converges: e w ðNÞ for the aperture and e v ðNÞ for the fluid velocity. The results are shown in Figs. 3 and 4. It is notable that both dw and e w , as well as dv and e v , provide estimates of a similar order for a fixed N. Thus, both of these error measures can be employed in the accuracy analysis, although only the latter (e w;v ) in the cases when no exact solution is available. As such, e w ðNÞ and e v ðNÞ will be used in the following investigations (it should also be noted that this approach has previously been shown to be successful for the PKN/KGD models of HF [24,25,33]).

Impermeable solid: reference solutions
With a suitable measure for testing the solution accuracy in place we move onto examining the solution variant most frequently studied in the literature, the case with a zero valued leak-off function and withQ 0 ¼ 1. Although there is no analytical solution to this problem, due to its relative simplicity, it is commonly used when testing numerical algorithms. For this reason it is very important that credible reference data is provided for this case, which can be easily employed to verify various computational schemes. Both the viscosity and toughness dominated regimes (for different values of the material toughness:K Ic ¼ 1; 10; 100 f g ) will be investigated. In the second part of the paper, accurate and simple approximations of the obtained numerical solutions will be provided.
The results for the crack opening and fluid velocity convergence rates are shown in Figs. 5, 6, 7 and 8.
As can be seen, over the analyzed range of N, the computations are very accurate and converge rapidly as the mesh density is increased. In the viscosity dominated regime it can be seen that there is a lower sensitivity of e w and e v to the value of n, however even in the toughness dominated mode the dependence of e w on the fluid behaviour index becomes less pronounced asK Ic grows. A general trend can be observed, in that the convergence rate is magnified as the self-similar material toughnessK Ic increases. This is due to the fact that, for large values ofK Ic , the solution tends to the limiting case of a uniformly pressurized immobile crack with a parabolic profile. To explain this tendency we present in Figs. 9, 10, 11 and 12 some additional data for a single value of the fluid behavior index (n ¼ 0:5).
It is immediately obvious that forK Ic [ 2 the fracture aperture is almost entirely described by the leading term of its crack tip asymptotics (forK Ic ¼ 2 the maximal deviation between them is approximately 1 percent). For the fluid velocity it can be seen that, while the effect is not as substantial as for the aperture, the crack propagation speedv 0 does become a better predictor of the parameter's behaviour for larger values of the material toughness. Meanwhile, the fluid pressure increases with growingK Ic , eventually becoming uniformly distributed overr. As a result of the decreasing pressure gradient the velocity of the fluid flow is reduced. In Fig. 12 it can be seen that the fluid flow rate rapidly converges to the limiting case with growingK Ic , however the rate of convergence is greater for larger values of n. Indeed, as can be seen in Fig. 13, for n ¼ 1 the curves forK Ic ¼ 1 andK Ic ¼ 100 are indistinguishable, which is not the case when n ¼ 0.
In fact, the behaviour of the solution asK Ic ! 1 can easily be shown to take the form:   where q is defined in Table 2. As a result the computations become far more efficient in this case and the resulting solution is calculated to a far higher level of accuracy.
Combining the results shown above in Figs. 1, 2, 3, 4, 5, 6, 7 and 8, it is clear that the computations presented here achieve a very high level of accuracy for both the aperture and fluid velocity regardless of the crack propagation regime. When using N ¼ 300 the accuracy of computations can almost always be assumed to be correct to a level of at least 10 À7 for the fracture aperture, and 2:5 Â 10 À7 for the fluid velocity. In this way the obtained data constitutes a very convenient and credible reference solution when It is worth mentioning that the efficiency of computations achieved by this algorithm means that this high level of accuracy does not come at the expense of simulation time. The final algorithm requires fewer than 20 iterations to produce a solution. Simulation times are also very short with this scheme.

Conclusions
In this paper, the problem of a penny-shaped hydraulic fracture driven by a power-law fluid has been analyzed. Following an approach similar to that in [24,25]   The following conclusions can be drawn from the conducted research.
• The universal algorithm for numerically simulating hydraulic fractures, introduced in [24], can be successfully adapted to the case of a penny-shaped fracture. It enables accurate and efficient modelling of HFs driven by the power-law fluids in both the viscosity and toughness dominated regimes. • The key elements of the algorithm, which contributed to its outstanding performance, are: i) choice of proper computational variables, including the reduced fluid velocity, ii) extensive utilization of the information on the solution asymptotics, combined with a fracture front tracing mechanism based on the Stefan-type condition (speed equation), iii) application of the modified form of the elasticity operator (26), which has a non-singular kernel, that can easily be coupled with the new dependent variable -the reduced fluid velocity. • The newly introduced analytical benchmark solutions, with a predefined non-zero fluid leak-off, can be adjusted to mimic the HF behaviour for a class of power-law fluids in both the viscosity and toughness dominated regimes. These solutions can be directly applied to investigate the actual error of computations when testing various computational schemes. While this work has allowed for the creation of highly accurate numerical benchmarks, they are not in a form which can be easily utilized. In the second part of this paper, approximate formulae for the case of an impermeable solid, constituting a set of accurate and easily accessible reference solutions when investigating other computational algorithms, will be delivered. Additionally, a brief comparison with alternative benchmarks available in the literature will be performed.

Compliance with ethical standards
Conflict of interest The authors declare that these publicly funded arrangements haven't created a conflict of interest.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix 1: Limiting cases: Newtonian and plastic fluids
Newtonian fluid: n ¼ 1 In the case of a Newtonian fluid the majority of the results remains the same as in the general case (setting n ¼ 1), but a few constants and functions will take alternate forms. These are detailed below.  Fig. 13 The self-similar fluid flow rate for different values of the fracture toughness when the fluid behaviour index is: a n ¼ 0 and b n ¼ 1