A radial hydraulic fracture driven by a Herschel–Bulkley fluid

We analyse the influence of fluid yield stress on propagation of a radial (penny-shaped) hydraulic fracture in a permeable reservoir. In particular, the Herschel-Bulkley rheological model is adopted that includes yield stress and non-linearity of the shear stress. The rock is assumed to be linear elastic, and the fracture is driven by the point source fluid injection with a constant volumetric rate. The fracture propagation condition follows the theory of linear elastic fracture mechanics, and Carter’s leak-off law is selected to govern the fluid exchange process between the fracture and formation. The numerical solution for the problem is found using the algorithm based on Gauss-Chebyshev quadrature and Barycentric Lagrange interpolation techniques. We also construct an approximate solution with the help of the global fluid balance equation and the near-tip region asymptote. The latter approximation is computationally efficient, and we estimate its accuracy by comparing the primary crack characteristics such as opening, pressure, and radius with that provided by the full numerical solution. We present examples corresponding to typical field cases and demonstrate that the addition of yield stress can lead to shorter radius and wider opening compared to the corresponding case with a simpler power-law fluid rheology. Further, we quantify the limiting propagation regimes (or vertex solutions) characterised by dominance of a particular physical phenomenon. Relative to the power-law results, there are two new vertices that are associated with domination of yield stress: storage-yield-stress and leak-off-yield-stress. To understand the influence of various problem parameters, we utilise the constructed approximate solution to investigate the dimensionless parametric space for the problem, in which the applicability domains of the limiting solutions are quantified. This enables one to quickly determine whether the yield stress provides a strong influence for a given problem parameters.


Introduction
Hydraulic fractures are fluid-filled tensile cracks propagating in a solid material, such as Earth's crust, and their dynamics is driven by the high-pressure fluid injection. Hydraulic fractures exist in nature in the form of magmafilled dykes (Spence and Turcotte, 1985;Lister, 1990;Rivalta et al., 2015;Dontsov, 2016b) and fluid-filled cracks in glacier beds (Tsai and Rice, 2010;van der Veen, 2007). However, most often they are human-made, and are utilised in the oil and gas fields to enhance production of hydrocarbons (Economides et al., 1989(Economides et al., , 2002. Horizontal wells with multistage hydraulic fracturing have quickly become popular in low permeability formations because, in this case, the connection between the reservoir and wellbore has large area (Vishkai and Gates, 2019).
With respect to petroleum applications, hydraulic fracturing fluid is injected through a wellbore into the reservoir rock resulting in the breakdown. A number of parameters can be adjusted to optimise the treatment, such as volumetric injection rate and hydraulic fracturing fluid properties. Besides breaking the rock, the hydraulic fracturing fluid also carries proppant inside the crack channel, which prevents complete crack closure after shut-in. The fracturing fluid rheology and pumping schedule are typically engineered to achieve the following goals (Economides et al., 1989;Barbati et al., 2016): creation of the sufficient crack aperture for proppant placement, minimisation of proppant particle settling (Osiptsov, 2017), reducing the risk of the bridging , reduction of fluid leak-off rate into the permeable reservoir, decreasing the energy consumption required to pump fluid , and setting the stage for favourable conditions for flowback . rheological behaviour. The main aims of the present work are the following: (i) to implement a solver for calculating an accurate numerical solution for the radial fracture, (ii) to construct an approximate computationally efficient solution for rapid estimations, (iii) to derive limiting propagation regimes occurring in the model, (iv) to explore the problem parametric space and (v) to analyse variations of the fracture characteristics depending on the yield stress and leak-off intensity for different values of the flow index.
The paper is organised in as follows. Section 2 outlines the problem formulation and the governing equations. In Section 3, we describe the methodology for calculating numerical solution, and, after that, we provide insights into getting the simplified approximate solution. Section 4 revisits the known limiting propagation regimes for a radial crack in a permeable rock and introduces the new members associated with dominance of yield stress. In that section, we also establish the conceptual representation of the problem parameter space and discuss the solution trajectories. Section 5 presents analysis of the results, including normalisation of the governing equations, estimation of the admissible ranges for the dimensionless parameters, extensive exploration of the parameter space, as well as discusses examples for several typical field cases. Finally, we summarise the findings of the paper and discuss potential applications of the developed model in Section 6.

Model formulation 2.1. Problem definition
The present paper considers a radial (penny-shaped) hydraulic fracture model whose sketch is shown in Figure 1. The crack grows along the plane perpendicular to the far-field confining stress due to fluid injection through a point source. The fracture is axisymmetric, i.e. there is a symmetry relative to the axis passing through the source and perpendicular to the fracture plane. For the model construction, we introduce the in-plane polar coordinate system ( , ) with the origin at the source. Due to symmetry, the spatiotemporal characteristics do not depend on the polar angle .
The volumetric rate of fluid injection is constant and is denoted by 0 , so that the injected volume is equal to: inj ( ) = 0 . The ambient rock is taken as linear elastic with Young's modulus and Poisson's ratio . We assume that the size of the zone, where the dissipation processes associated with the rock failure happen, is small compared to other length scales realised in the model, e.g., linked with the viscous fluid flow and leak-off. Consequently, the linear elastic fracture mechanics (LEFM) theory can be applied to model quasi-static fracture propagation in a solid medium with toughness . The fracture surface is exerted by the fluid pressure ( , ) from the internal side, while the far-field confining stress is . The radial fracture model is fully characterised by the opening profile ( , ), radius ( ), net pressure profile ( , ) = ( , ) − and efficiency parameter ( ) = crack ( )∕ inj ( ), where crack ( ) is the fracture volume.
Hydraulic fracturing fluid has Herschel-Bulkley rheology given by the constitutive relation: where is the shear stress, 0 is the yield stress,̇ is the shear rate, is the consistency index, and is the flow index. When the yield stress is non-zero ( 0 > 0), and the shear stress is less than the yield stress ( < 0 ), the unyielded (or plug) zone is formed where the fluid behaves like a solid. Since the shear stress is linear across the aperture and equals zero at the centerline, the plug zone is formed in the middle of the crack channel, and its width 2 ( , ) is apriori unknown function of position and time ( Figure 1). The Herschel-Bulkley rheological model has three limiting cases: (i) Newtonian fluid when 0 = 0 and = 1 ( corresponds to the dynamic viscosity), (ii) power-law fluid when the yield stress is zero 0 = 0, and, finally, (iii) it reduces to Bingham model when 0 ≠ 0 and = 1. The fluid flow inside the fracture channel is controlled by lubrication theory (Batchelor, 1967), its velocity and the flow rate are denoted by ( , ) and ( , ) = ( , ) ( , ), respectively. The fracture and fluid fronts are assumed to coincide in the model, which means that the lag filled by the vapour (Garagash and Detournay, 2000) or pore fluid (Detournay and Garagash, 2003;Kanin et al., 2020d) adjacent to the fracture front is negligibly small.
The host rock is porous (with porosity ) and permeable (with permeability ). The fluid exchange process between the fracture and the formation is taken in the form of the fracturing fluid filtrate leak-off, and its properties are assumed to be similar to water (viscosity and compressibility ). The leak-off is modelled by Carter's law (Carter, 1957) which presents the one-dimensional and pressure-independent process. Following this law, the leak-off rate ( , ) is proportional to the inverse square root of the "exposure time", i.e. the interval between the current time and the time instant at which the fracture tip was at a given point on the fracture plane. The proportionality coefficient is called Carter's leak-off coefficient . When the filter-cake does not form, and the leaked fluid properties are identical to the pore fluid characteristics (Kovalyshen, 2010;Kanin et al., 2020d), the Carter's coefficient is calculated from: where is the far-field pore pressure, and = ∕( ) is the diffusivity coefficient. The reader can find the expressions for Carter's leak-off coefficient in other cases, e.g., filter-cake presence, in (Economides et al., 1989).

Governing equations
In this section, we discuss the system of governing equations. It is formulated for the unknown crack radius ( ), opening ( , ) and net fluid pressure ( , ) profiles. The fracture characteristics depend on time , distance from the point source , the set of material parameters: as well as the injection rate 0 and the yield stress 0 . In equation (2), ′ is the plane strain elastic modulus, ′ and ′ are toughness and viscosity parameters, and ′ is the leak-off parameter.

Crack elasticity
The elasticity equation expresses the net fluid pressure ( , ) in terms of the crack aperture ( , ) and radius ( ) (Arin and Erdogan, 1971;Cleary and Wong, 1985;Savitski and Detournay, 2002): The integral kernel ( , ) in equation (3) has the following form: where K( ) and E( ) are the complete elliptic integrals of the first and the second kind, respectively.

Fluid flow
Based on lubrication theory (Batchelor, 1967), we write out the width-averaged mass conservation equation for the fluid flow inside the crack channel: The Reynolds equation is obtained by substitution of the fluid velocity (7) into the continuity equation (5).

Fracture propagation
The LEFM theory states that for quasi-static propagation of a hydraulic fracture, the stress intensity factor matches rock toughness: = .
The alternative form to this condition can be expressed in terms of the near-tip asymptotic behaviour of the fracture opening ( , ) (Irvin, 1957): If the crack propagation velocity = ∕ =̇ is zero, the tip asymptote has the following form: where ′ is an unknown scaled stress intensity factor, whose value is less than the critical value ′ .

Boundary conditions
At the fracture inlet, the volumetric flow rate should be equal to the specified value 0 : In turn, the crack tip is characterised by zero opening and no-flux condition:

Global fluid volume balance
By integrating the continuity equation (5) with respect to time and distance and taking into account the boundary conditions (9), (10), we derive the global fluid balance equation: For implementation of the numerical solution, we utilise the elasticity equation with an extended interval: where the opening profile ( , ) continues symmetrically to negative , and the integral kernel has the following form: It can be easily shown that this formulation of the elasticity equation is identical to (3).

Numerical solution
By following the approach of Liu et al. (2019), we present gradient of the fracture opening profile in the form: where ( ) is the required unknown function, and the weight function ( ) includes the tip behaviour of The function ( ) is odd, which ensures that the fracture opening profile ( ) is symmetric. Further, we perform discretisation of the computational domain ∈ [−1, 1] by introducing two systems of nodes corresponding to the chosen weight function (Viesca and Garagash, 2018) . Now, we move on to the discretisation of the governing equations. For brevity, the matrix notation is utilised. Using the derivations of Liu et al. (2019), we write out the discretised form of the elasticity equation (14): where "×" denotes matrix multiplication; G is the elasticity quadrature matrix, and its representation composed of the following matrices: Ψ ( ) − Ψ (0) ′ , Ψ ( ) = cos [( + 1) ]∕( + 1), = arccos , ′ = 2∕ ⋅ sin ′ ∕ sin ′ ( + 1)∕ (in this section, symbol is a summation index, not formation permeability); Let us now focus on the discretisation of Reynolds equation (16), (17). Similarly to Liu et al. (2019), we integrate equation (16) with respect to from each node to the tip, i.e. = 1. As a result, we obtain the following discretised form: where the matrices S and R are: In equation (22), the last term on the right-hand-side, i.e., the integral with respect to time, is computed using Simpson's rule.
Finally, we write out the propagation condition in the matrix form (see details in ) and differentiate it with respect to time: Here, we should make a comment regarding the application of the propagation condition (23). The LEFM asymptote (19) is always valid near the fracture tip. In certain situations, e.g., vanishingly small toughness, its spatial applicability domain can be tiny. Even when the vanishing LEFM region can not be accurately resolved by the method, Chebyshev's nodes high density near the tip accurately captures the dominant asymptotic behaviour.
Further, we combine equations (21), (22), (23) into the system of ordinary differential equations (ODEs) which can be written as: where the vector consists of the unknown parameters, and we use the fact that =̇ .
The value of should be odd since the opposite choice leads to the presence of the infinite components in the elasticity matrix G (20) (see in the denominator). The function ( ) is odd, and the the vector has in the following form The total number of the independent unknown parameters is ( − 1)∕2 + 1, i.e., radius and independent components of the vector , meaning that the system of ODEs should be composed of the first ( − 1)∕2 − 1 discretised Reynolds equations (21), i.e., corresponding to the complementary nodes 1 , … , ( −1)∕2−1 , the global fluid balance equation (22), and the propagation condition (23). We choose = 101, and the storage-viscosity dominated regime ( -vertex) for the power-law fluid is taken as an initial condition for all considered values of the flow index (the fracture properties in this regime are provided by Section 4.1). Despite the fact that for < 0.5 the early time solution corresponds to the storage-toughness regime (vertex, see details in Section 4.1), the algorithm works fine with the selected initial condition, even though requiring a certain time-span to adjust the solution to the actual trajectory. The specified time interval is discretised uniformly on a logarithmic mesh, and the solver is applied within each time segment. Python programming language is used for implementation of the numerical algorithm, and the system of resulting ODEs is solved via "solve_ivp" function of SciPy library (Virtanen et al., 2020).

Rapid approximate solution
This section outlines a rapid approximate solution for the problem. This approach is based on the idea that the near-tip region behaviour predominantly determines the finite fracture characteristics. As a result, the crack opening profile is presented in the following form (Dontsov, 2016a): where ( ) is the opening asymptote near the fracture tip stemming from the solution for the semi-infinite geometry, and it is a function of the distance from the front , material parameters (2), yield stress 0 and time (through ( ) anḋ ( )). In equation (25), we also utilise parameters and̄ : the first one originates from comparisons with the accurate solution for limiting cases (will be defined later), and the second one is a slowly varying parameter from the tip asymptotic solution, i.e. ( ) ∝ ̄ . Since the fracture radius is a power-law function of time in the known limiting propagation regimes (see e.g. a review paper (Detournay, 2016)), we assume that ( ) ∝ , where is a slowly varying function of time. Consequently, the inverse radius function is expressed in the form: 0 ( ) = 1∕ . By substituting the opening profile (25) and 0 ( ) into the global fluid balance equation (18), we obtain: The integrals in the above equation can be estimated using special functions, and the result becomes: where 0 ( ; , ) = ( , ) − ( ; , ), ( , ) is the beta function, and ( ; , ) is the incomplete beta function. As can be seen from the above result, one of the building blocks in the approximate solution is the tip asymptote ( ). We utilise an approximate near-tip region solution developed in Bessmertnykh and Dontsov (2019), which is computed using a combination of the yield stress dominated tip solution ( ) and the one corresponding to the power-law fluid ( pl ): where the function (̃ , ) is defined as follows (Dontsov and Kresse, 2018): The parameter̄ in the global fluid balance equation (27) can be computed from the following relation: where the function 1 is defined in the equation (29), and it enters approximation for the power-law solution: pl ( ) ∝ ( 1 +1)∕2 , where is the distance from the tip. Further, we solve the non-linear algebraic system of equations composed of the global fluid balance (27) and the tip asymptote (28). Let us discuss the overall approach in more details. Initially, the desired time interval is discretised uniformly on a logarithmic scale, and the target parameters  = , , are computed iteratively step-by-step for each time instant. Let us consider the example of time step . From the previous result at −1 , we know the values of  −1 = −1 , −1 , ( ) −1 , and they are used as an initial guess for the current computation. We calculate the solution for two time instants: and * , where the latter is a fraction of the current time moment, e.g., it can be * = 0.9 ⋅ . Further, the value of initially taken as −1 is updated by using the equation: = log ( )∕ log ( ) = log ( ) − log ( * ) ∕ log ( ) − log ( * ) . After that, the system is solved again with the new initial guess presented by the parameters from this iteration. The process continues until the convergence in the value of is reached. It is also necessary to mention that for the first time step ( = 1), we take the initial guess 0 = (2 + 2)∕(3 + 6), 0 and ( ) 0 corresponding to the storage-viscosity dominated limiting propagation regime (Section 4.1).
Once the solution  is computed, the fracture opening profile ( , ) is evaluated using equation (25), while the pressure ( , ) profile is calculated based on the transformed elasticity (13) equation (Dontsov, 2016a): where the function  is evaluated numerically. Finally, we also compute the fracture efficiency as: where the function  is introduced in equation (27). For brevity of presentation, the exact form of the function is given in Appendix A. The calculation procedure for also requires the results of Section 4.1, in which the limiting propagation regimes are outlined for the present model.

Limiting propagation regimes
Two different physical mechanisms govern the propagation regime of a finite hydraulic fracture (see a review paper Detournay (2016) and references therein). The first one is related to the distribution of the total dissipated energy between the creation of new fracture surface and viscous fluid flow including the movement of the solid plug inside the fracture channel. The latter energy component, namely movement of the solid plug, is included since the fracturing fluid has non-zero yield stress. In this case, the whole fracture can be filled with the un-yielded solid material. The second mechanism is the partitioning of the injected fluid volume between the fracture and the host permeable rock (due to leak-off). When the fracture grows, the partitioning of the dissipated energy and the injected volume change over time, leading to the emergence of the various limiting propagation regimes with one dissipation (out of two) and one storage (out of two) mechanisms at different time moments. The leak-off parameter ′ influences distribution of the injected fluid volume, while the viscosity ′ , toughness ′ , and the yield stress 0 have an effect on the partitioning of the dissipated energy.
Let's now consider scalings associated with the limiting propagation regimes. First of all, we present the main crack characteristics as: where = ∕ ( );  1 ,  2 ,  3 are dimensionless evolution parameters depending on , material parameters (2), 0 and 0 ; ( ) is the length scale (the same order as the crack radius) and ( ) is a small dimensionless parameter with the meaning of a characteristic strain in the rock.
Further, we substitute the expressions from equation (32) into the system of governing equations written in terms of the normalised distance from the source , i.e. the elasticity (13), Reynolds (16), (17), global fluid balance (18) equations and the propagation condition (19): • Elasticity: • Reynolds: • Global fluid balance: • Propagation condition: Here, we utilise the function 0 ( , ) = 0 ( )∕ and introduce five dimensionless numbers:  and  parameters quantify the fluid volume stored in the fracture and the volume leaked into the ambient permeable rock, respectively. In turn, the numbers ( ,  ) and  are related to the energy dissipation in the fluid flow inside the crack channel in overcoming fluid viscosity and solid plug yield strength and in the brittle rock failure, correspondingly. To derive various scalings, i.e., ( ) and ( ) in (32), we set one out of two fluid storage parameters ( ,  ) be equal to one. Similarly, among the three parameters responsible for energy dissipation ( ,  ,  ), one should be made equal to one. The remaining three dimensionless groups are the evolution parameters  1 ,  2 ,  3 mentioned earlier.
When the evolution parameters are approximately zero in a given scaling, the corresponding limiting propagation regime is realised. From the scaling analysis, we can determine only the dimensional multipliers for the radial crack characteristics, , , , in a given 'i-th' scaling. In order to approximately quantify the opening and pressure profiles, we rely on the approximations of Section 3.2: where the subscript indicates the particular limiting regime, and the function  is defined in equation (30). In order to express the prefactors values * , * through , we substitute ( ) and ( ) ( ) into the global fluid balance equation (27) (accounting the conditions inherent to the analysed regime) and the appropriate near-tip region asymptote for the crack opening. Further, the parameter can be found with the help of the accurate numerical solution (Section 3.1). Before discussing different scalings and limiting propagation regimes, we introduce the relative errors in the estimation of the radius and opening characteristics provided by the simplified approach (34) compared to the fully numerical solution: where Δ( ) is the radius error, Δ( 0 ) is the error of width at the inlet, Δ( ) is the maximum relative width error. The subscript denotes the considered limiting regime, and the bar symbol means the fully numerical solution.
We begin with the storage-viscosity scaling and the −vertex solution (subscript " "). We set  =  = 1 and obtain the following formulas for the length scale and the small parameter: Further, we substitute and into the remaining dimensionless groups: where we reassign  ,  ,  as the dimensionless leak-off  , yield stress  and toughness  . After that, we obtain formulas for the prefactors * , * : where = (2 + 2)∕(3 + 6),  = ( ,̄ ) (see equation (27)), the value of is provided by equation (29), and = 2∕(2 + ) is taken from the storage-viscosity tip asymptote (Dontsov and Kresse, 2018). The parameter is obtained by fitting the approximate solution above (34) and (40) to the full numerical solution (Section 3.1) in terms of and (0, ). We computed for different values of flow indexes ∈ [0.25, 1], and the obtained curve is shown in Figure 2(a). The relative errors defined by equations (35)-(37) are depicted in Figure 2 Figures 3(a) and 3(d) present comparison of the approximate normalised profiles of ( , ), ( , ) (dashed lines) with the corresponding properties evaluated by the fully numerical solution (solid lines) for = 1 (light-blue colour) and = 0.3 (blue colour). We use the normalisation of the width by ( ) ( ), pressure by ( ) ′ as per (32) and of the radial distance from the source by ( ) scale; these scales account for the limiting solution dependence on time and largely on the flow behaviour index . The characteristics for the Newtonian fluid are very close to each other, while the accurate and approximate profiles for = 0.3 intersect for both opening and pressure; moreover, there is a noticeable difference for the pressure values near the fracture front.
Let us consider the case of the leak-off-viscosity scaling and the corresponding̃ −vertex solution (subscript̃ ). This scaling corresponds to  =  = 1. Hence, The evolution parameters have the meaning of the dimensionless storage ̃ , yield stress ̃ and toughness ̃ : For thẽ −vertex, we derive the following prefactors * ̃ , * ̃ : where ̃ = 1∕4, ̃ is given by equation (29), and we also utilisē ̃ = (4 + )∕(4 + 4 ) known from the leak-offviscosity near-tip region asymptote (Dontsov and Kresse, 2018). The parameter ̃ featured in the ̃ profile is found from the comparison of ̃ (0, ) (equation (34) combined with (43)) with the numerical solution (Section 3.1). This  (e)), and̃ (panels (c) and (f)) calculated by the fully numerical (solid lines) and approximate (dashed lines) approaches. We apply ( ) ( ), ( ) ′ and ( ) dimensional prefactors for the normalisation of ( , ), ( , ) profiles, and distance from the source, and here, the subscript denotes the analysed regime. In the case of and̃ vertex solutions, we look at = 1 or = 0.3 and show the computed properties in plots (a), (d) and (b), (e) by the light-blue or light-green and blue or green colours, correspondingly. The normalisation coefficients take into account the value for the flow behaviour index. We use olive and maroon colours for the yield stress dominated regimes and̃ in the charts (c), (f). approach yields to the values of ̃ ( ) presented in Figure 4(a). Moreover, we also estimate the relative error in the crack aperture values ̃ ( , ) provided by equation (25) with ̃ ( ) (Figure 4(b)), and we utilise equations (36), (37) for Δ( 0̃ ) and Δ( ̃ ).
In the panels (b), (e) of Figure 3, we demonstrate the spatial variations of the width and pressure calculated by the accurate solution (solid lines) and provided by the approximate approach (dashed lines), i.e., ̃ ( , ), ̃ ( , ), for two cases of the flow behaviour index = 1 (light-green colour) and = 0.3 (green colour). The dimensional prefactors ̃ ( ) ̃ ( ), ̃ ( ) ′ , ̃ ( ) are applied for the normalisation of the aperture, pressure and distance from the source, and they include the value of the flow behaviour index. We can mention that for both analysed values of the flow index, the accurate and approximate opening and pressure profiles are very close to each other.
Next, we discuss the storage-toughness scaling that corresponds to the −vertex solution (subscript " "). Here, we require  =  = 1 leading to: and the evolution parameters can be interpreted as the dimensionless leak-off  , yield stress  and viscosity  : 1∕ (5 ) .
The LEFM asymptote governs the tip behaviour in the storage-toughness limiting propagation regime, i.e.̄ = 1∕2. Moreover, in this case, the crack opening profile is elliptical (Savitski and Detournay, 2002). Consequently, the simplified form of the radius crack opening profile (25) allows describing it accurately by taking = 1∕2. Using the global fluid balance equation (27) and the "square-root" asymptote (Irvin, 1957), we find out the following coefficient in the radius and opening profiles: Next, we construct the leak-off-toughness scaling and focus on thẽ −vertex solution (subscript "̃ "). We put  =  = 1, and get the expressions for the length scale ̃ and the small parameter ̃ : The evolution parameters are identified as the dimensionless storage ̃ , yield stress ̃ and viscosity ̃ : Similarly to the −limiting propagation regime, the crack opening behaviour near the tip is described by the LEFM asymptote (̄ ̃ = 1∕2), while the whole profile is elliptical ( ̃ = 1∕2). The derivation of the prefactors for the radius and opening leads to: * = Finally, next we reach the cases that correspond to dominance of yield stress. The storage-yield-stress scaling (subscript " ") corresponds to the case:  =  = 1, which results to: In turn, the evolution parameters are the dimensionless leak-off, viscosity and toughness: Further, we should determine the prefactors * , * from the global fluid balance (27) and the yield stress dominated near-tip region asymptote (Bessmertnykh and Dontsov, 2019) where  = ( ,̄ ) and̄ = 1. Using the radius and wellbore opening characteristics computed via the accurate numerical solution (Section 3.1), we find out that = 1.029 provides the most accurate approximation. This value together with equations (25), (52) result in Δ( ) = 3% (equation (35)) while the (0, ) is captured precisely. At the same time, the absolute difference between the fully numerical and approximate opening profiles ( , ) normalised by the opening at the wellbore (0, ) can reach Δ( ) = 5% (equation (37)), which demonstrates that the accuracy of the approximation (34) reduces for such values of . The last case is the leak-off-yield-stress scaling and thẽ −vertex solution associated with it. We set  =  = 1, which results in the following parameters: Here, the evolution parameters are identified as the dimensionless storage ̃ , viscosity ̃ and toughness ̃ : Then, we retrieve the coefficients for the radius and opening (equation (34)) based on the global fluid balance (27) and the tip asymptote (Bessmertnykh and Dontsov, 2019) Similarly to thẽ −vertex solution, we estimate ̃ via the tuning of the simplified opening profile (25) to the accurate numerical solution (Section 3.1). We compute ̃ = 1.077, and it provides the best approximation which reproduces precisely the fracture opening at the wellbore ̃ (0, ) and gives the error in the spatial variation ̃ ( , ) (in relation to ̃ (0, )) up to Δ( ̃ ) = 6% (equation (37)). Figures 3(c) and 3(f) depict the aperture and pressure profiles obtained by the accurate method (solid lines) and that computed by the approximate one (dashed lines) for the yield stress dominated regimes (olive colour) and̃ (maroon colour). For both vertices, we observe the intersection between the fully numerical and simplified solutions and a noticeable difference between them, which is caused by the loss of accuracy of the width approximation utilised to construct the simplified solution (34). One possible mitigation of this issue is to introduce a more sophisticated expression for fracture width that captures such behaviour for the and̃ limits. The existing approximation was initially proposed for smaller values of̄ and that correspond to Newtonian fluids. However, for the purpose of this study, i.e., for exploring the problem parameter space, such accuracy is acceptable, especially since the boundaries of the limiting solutions vary on a logarithmic scale.
The fracture characteristics observed in the limiting propagation regimes are summarised in Appendix B. In addition to the dimensional form, we also present the solutions in the normalised form obtained by the application of the -scaling introduced in Section 5.1.

Representation of the problem solution space
The parametric space of solutions for the discussed radial fracture model can be conceptually presented in the form of a hexahedral pyramid, and the limiting propagation regimes are located at its vertices. The sketch of the pyramid ̃ ̃ ̃ is shown in Figure 5. Each edge of the pyramid links two vertex solutions, e.g., and , and it is possible to introduce a characteristic transition time between them by solving the following equation: ( ) = ( ) → = . Furthermore, the evolution along the edge is controlled by a single evolution parameter, which can be expressed as a , while  ( ) and  ( ) are evolution parameters related to this edge. By introducing the dimensionless time = ∕ , we can write out the relations: In the current problem, we have 9 transition time scales, and only 3 of them are independent, e.g., , ̃ , . Consequently, location of the solution inside the parametric space can be expressed as a function of the three dimensionless times = ∕ , ̃ = ∕ ̃ , = ∕ . Further, we introduce the parameters and as: and they can be interpreted as the dimensionless leak-off and yield stress numbers. By taking the set of parameters = , , , we can express ̃ and in terms of them in order to characterise the solution trajectories. Along each trajectory, the numbers and are constant while varies.
When the flow index is greater than one-half ( > 0.5), all solution trajectories start from the limiting propagation regime (early-time asymptote) since all the dimensionless groups (39) vanish when time tends to zero. However, their destination points (large-time asymptote) can be different: One can explain the endpoints of the solution trajectories by looking at the evolution parameters of the corresponding regimes (e.g., (45), (48), (51), (54)) which go to zero when time tends to infinity or they are identically zero due to values of , . For intermediate time intervals, the solution trajectory can be attracted to either of the vertices, and its behaviour strongly depends on the values of and . For example, when the yield stress is absent ( = 0), and̃ regimes can be considered in the general solution as intermediate asymptotes when the leak-off is small ( ≪ 1) and large ( ≫ 1), correspondingly. As another example, we can take zero leak-off case ( = 0) with the yield stress ( > 0), and -vertex attracts the solution trajectory when ≪ 1.
We want to emphasise that in the general case when the dimensionless numbers and are both non-zero, the large time asymptotic behaviour is always dominated by leak-off and yield stress, i.e.,̃ -vertex, regardless of the value of the flow index . The examples of the solution trajectories for > 0.5 and < 0.5 are shown in Figure 5 by dashed and dotted lines, respectively.
Finally, we should mention that the chosen set of parameters, i.e. = , , , is not applicable for the description of the solution trajectories behaviour for = 0.5 since they are not properly defined in this case. One possibility to mediate this issue is to introduce another set of variables, e.g., , where ′ is the dimensionless toughness, and ′ is the dimensionless yield stress. Nevertheless, since the problematic region consists of solely a single point, we proceed with the chosen set of parameters = , , .

Admissible ranges of the dimensionless problem parameters.
Let us start by rewriting the problem formulation in terms of the dimensionless variables. For the normalisation of the governing equations, we rely on the -scaling consisting on the following parameters (Detournay, 2016): where is the transition timescale between and limiting regimes (Section 4.2), = ( ) is the lengthscale, and = ( ) is the small parameter (see equation (38)): . (58) The elasticity (13), Reynolds (16), (17), global fluid balance (18) equations and the propagation condition (19) written in terms of the dimensionless distance and using (57) and (58) can be written in the normalised form as: • Elasticity: Ω ; • Reynolds: where the dimensionless leak-off and yield stress are defined in equation (56), and the dimensionless inverse radius function is given by 0 ( , ) = 0 ( )∕ . For completeness, we also provide the definitions of and through the characteristics of the -scaling: • Propagation condition: We can use the numerical schemes presented in Sections 3.1 and 3.2 for the solution of the normalised system of governing equations assuming ′ = ′ = ′ = 0 = 1, ′ = 1∕4 , 0 = , and = .
We take three values of the flow behaviour index = {0.6, 0.8, 1} and compute the intervals for the governing parameters and by varying each of the aforementioned dimensional input parameters independently. We illustrate the evaluated domains in Figure 6(a). It can be noticed that the domain for = 0.6 (green one) includes the analogous zone for = 0.8 (red one) which, in turn, includes the region for = 1 (blue one). To understand the structure of the non-dimensional parametric domain for a fixed (how different dimensional input parameters affect the locations inside it) we look closely at the example of the Bingham fluid ( = 1) with the yield stress 0 = 15 Pa (yellow decagon in Figure 6(b); such decagons for different 0 ∈ [0, 15 Pa] fill completely the whole blue domain). We identify two sub-domains (octagons) corresponding to the minimum and maximum values of the plane-strain elastic modulus and frame them by coloured dash-dotted lines (Figure 6(b)); the octagons belonging to the intermediate values of ′ are located between them. Consequently, the resultant interval for the dimensionless leak-off number is controlled by the predefined range for ′ , while the interval for is primary governed by the yield stress parameter and then by ′ . Next, each octagon is limited by two hexagons corresponding to the limiting values of the rock toughness. In Figure  6(c), we show the octagon related to the minimum value of the plane-stain elastic modulus ( ′ = 10 GPa) by the orange colour and highlight the sub-domains for the minimum and maximum values of by the coloured dashed lines. One can also examine the interior of each hexagon and reveal the parallelograms corresponding to the extreme values of the consistency index and rock permeability which we emphasise by the coloured dotted boundaries and the coloured line fill in Figure 6

Analysis of the parametric space
This section considers analysis of the parametric space for a radial crack driven by a fluid with Herschel-Bulkley rheology. We utilise the dimensionless problem formulation (Section 5.1) in which the crack characteristics governed by the following dimensionless variables: time ( ), distance from the source ( ), leak-off ( ) and yield stress ( ).
In this section, we focus on the results related to the following values of the flow behaviour index: = 1 and = 0.3, and the examination includes the following two main components. Firstly, we identify the applicability domains of the limiting propagation regimes inside the parameter space and draw them as the regime maps. The construction of such kind of maps helps to understand the crack propagation conditions for the given problem parameters. The validity zones of the vertices are defined according to the following criterion (the same was used in Dontsov (2016a)): Equation (59) states that the combination of the relative differences of the opening at the wellbore and of the fracture radius between the numerical solution and the -th limiting regime should be less than 1% for the solution being considered at the vertex. It is important to highlight that for the purpose of this section we apply the simplified approach (Section 3.2) for calculations due to its computational efficiency. In Appendix C, we present validation of the corresponding approximate Table 1 The validity boundaries of the limiting regimes obtained by setting edge = const.
solution by comparing its predictions to that provided by the accurate numerical method (Section 3.1). We apply criterion (59) to estimate the applicability boundaries of all vertex solutions numerically. Moreover, to provide quick analytic estimates, we fit the obtained points by the appropriate analytical functional dependencies that are derived from the consideration of the transition timescales between the limiting regimes. As an example, let us consider the validity boundaries of the storage-viscosity ( ) and storage-yield-stress ( ) regimes framing the transition (another interpretation of the transition is the -edge of the solution space outlined in Section 4.2). The time normalised by the transition timescale has the form: = (2+ )∕(2 ) , providing the following relation for the discussed boundaries: ∝ −2 ∕(2+ ) . The prefactors are estimated numerically via the fitting procedure for and vertices separately. We summarise all possible relationships arising in the model as the applicability boundaries of the limiting solutions in Table 1.
The second element of the analysis is the examination of the variations of the main time-dependent crack characteristics such as radius ( ), opening at the wellbore Ω(0, ), pressure at the half-radius Π(1∕2, ), and efficiency ( ) with the dimensionless leak-off and yield stress parameters.

Bingham fluid ( = 1)
We start with a radial crack driven by the Bingham fracturing fluid ( = 1) in an impermeable formation ( = 0). The regime map is shown in Figure 7(a) in the coordinates ( , ). Generally, when the leak-off number = 0, only three storage vertices can be realised: viscosity ( ), toughness ( ), and yield stress ( ). For the flow behaviour index > 0.5, the problem solution evolves from to ( > 0) or to ( = 0) as it is confirmed by Figure 7(a) for = 1. Moreover, for > 0.5 and non-zero yield stress, the solution passes through the toughness dominated regime ( ) over an intermediate time range if yield stress ( ) is small enough. For example, the -vertex approximates the general solution along certain time intervals when ≲ 10 −6 for = 1, while the solution is dominated by dissipation in the fluid at all times, i.e., solution is given by to transition bypassing -vertex, when the yield stress is large enough, e.g., ≳ 1 for = 1 (Figure 7(a)). The time domain of the storage-viscosity regime gradually shrinks with increasing ≳ 10 −1 for = 1 (Figure 7(a), whereas the time domain of the storage-yield-stress regime expands.
Further, we look at the solution trajectories in the parameter space for = 1 and = 0 corresponding to the following values of the dimensionless yield stress number: = 10 −10 , 10 −5 , 1, 10 5 (the grey dash-dotted lines in Figure 7(a)). We calculate evolution of the time-dependent crack parameters (radius ( ), wellbore opening Ω(0, ), and pressure Π(1∕2, )) and plot them normalised by the -vertex solution in Figures 7(b) -(d). One can notice that the solutions are approximately independent of the yield stress and given by the -vertex solution during an initial time period which duration depends on the magnitude of . After that initial propagation stage, the presence of the non-zero yield stress leads to a radial crack with the smaller radius (slower fracture growth), larger wellbore opening and higher pressure compared to the = 0 case (see panels (b) -(d) in Figure 7 where the solution for = 10 −10 coincides with the = 0 case within the chosen time interval). Moreover, the relative differences between the computed properties for Newtonian ( = 0) and Bingham > 0 fluids grow with the increase of the magnitude of the dimensionless yield stress number . Various asymptotic regimes of the solution as discussed above with the help of the parametric map (Figure 7 Figure 8(a)), and plot time-dependent crack characteristics normalised by the storage-viscosity limiting solution in Figures 8(b) -(d). The results demonstrate that the non-zero yield stress affects radial fracture growth qualitatively similar to the Bingham fluid case: (i) the non-zero yield stress impacts the crack properties after a definite time range from the initiation which duration is governed by the yield stress number value; (ii) when the solution for the Herschel-Bulkley fluid deviates from that of for the Newtonian fluid, it is characterised by the reduced crack radius and increased maximum opening and pressure; (iii) the relative difference between the compared solutions, i.e., for > 0 and = 0, grows with increasing .

Permeable rock
Next, consider the radial crack problem in a permeable reservoir ( > 0). Here, we vary both governing parameters and to fully explore the problem parameter space. For each selected value of the flow behaviour index, i.e., = 1, 0.3, we carry out calculations for the following values of the dimensionless yield stress number: = 10 −10 , 10 −5 , 1, 10 5 and construct the regime maps in the coordinates ( , ). Then, we investigate variations of the major time-dependent crack properties for = 1 and different values of . Bingham fluid ( = 1) Figure 9 shows regime maps for a radial crack driven by fluid with = 1 and various values of yield stress. When the flow behaviour index > 0.5, the problem solution evolves from the storage-viscosity ( ) to the leak-off-toughness (̃ ) regime for the power-law fluid ( = 0) or to the leak-off-yield-stress (̃ ) vertex for the fluid with yield stress ( > 0). The applicability domains of the limiting propagation regimes for > 0 are shown by various colours in Figure 9, while the boundaries of the vertices for = 0 are depicted by coloured dashed lines. When the yield stress is absent, the storage-toughness ( ) and the leak-off-viscosity (̃ ) regimes are realised as the intermediate asymptotes, and for = 1, they emerge along the intervals ≲ 10 −15 and ≳ 10 3 , respectively (see the coloured dashed lines in Figure 9). For a non-zero yield stress, the yield stress dominated regimes ,̃ emerge earlier in the solution with an increase of (see the shifting of corresponding domains from Figure 9(b) to (d)). As a result, the validity zones of the toughness dominated regimes ,̃ shrink with growing , and eventually, they disappear completely (see Figure 9  the flow behaviour index within the interval ∈ (0.5, 1] leads to further separation of ,̃ , ,̃ regimes from the "center" of the map, i.e., the boundaries of and vertices framing the transition move to the left (i.e., smaller times) and to the right (i.e., larger times), respectively, while the boundaries of̃ and̃ regimes framing thẽ ̃ transition shift up-left and down-right, correspondingly. The value of starts to influence locations of the boundaries of and̃ regimes related to and̃ ̃ transitions only for large values of . Further, we take a closer look at variations of the time-dependent properties of the radial crack driven by Bingham fluid ( = 1) with respect to parameters and . (At this point, we should note that the following analyses are qualitatively applicable to the other values of the flow behaviour index inside the interval 0.5 < ≤ 1.) We consider the following set of the dimensionless leak-off values = 10 −20 , 10 −10 , 1, 10 10 and perform computations for both zero (Newtonian fluid) and non-zero (Bingham fluid) yield stress ( = 1). The obtained results are demonstrated in Figure 10. The solid black lines show the solution with = 1, while the case = 0 is presented by the dashed grey lines. First, we discuss the dependence of the solution on the leak-off number for a fixed value of the yield stress number ( = 0 or = 1). We observe that the solution is independent of (and thus given by the impermeable case, = 0) during an initial time period of fracture propagation which duration depends on and values. Outside this initial time span, the increase of leak-off leads to smaller radius, opening at the wellbore, and crack efficiency, whereas the pressure becomes larger (Figure 10). However, the pressure behaviour can differ from the description provided above for the non-zero yield stress case in which pressure for all values tends to reach the same yield stress dominated asymptote (see the cyan line in Figure 10(c)).
Next, we focus on the effect of the yield stress by comparing the solutions for Newtonian ( = 0) and Bingham ( = 1) fluids at a fixed value of the leak-off number (Figure 10). These solutions are insensitive to the fluid yield stress and thus coincide at the early stages of crack growth. Moreover, when the leak-off is very large ( = 10 10 ), two solutions match in the entire time domain for the radius and efficiency (Figures 10(a), (d)). At larger time the fluid yield stress ( = 1) leads to the reduced value of the radius and increased opening at the wellbore, pressure, and efficiency compared to the Newtonian fluid case ( = 0). The time domains where the fracture propagation is dominated by different limiting regimes can be identified in Figure 10 by comparing the computed fracture properties with the vertices (Appendix B) shown by the coloured dashed lines. We note that the limiting solutions for the radius are identical in all leak-off dominated regimes and thus shown in Figure 10(a) by a single (orange) colour. Similarly, the limiting solutions for pressure is all yield stress dominated regimes in Figure 10(c) are depicted by same colour (cyan) line.

Herschel-Bulkley fluid with = 0.3
We now turn to the discussion of the parametric dependence of the solution for a radial crack propagation in a permeable rock ( > 0) due to the injection of Herschel-Bulkley fluid with the flow behaviour index = 0.3. Figure  11 shows the computed propagation regime maps. When the flow behaviour index < 0.5, the problem solution starts at the storage-toughness ( ) regime and finishes at the leak-off-viscosity (̃ ) regime for the power-law fluid ( = 0) or at the leak-off-yield-stress (̃ ) regime for Herschel-Bulkley fluid ( > 0). The validity zones of the limiting regimes corresponding to > 0 cases are shown by various colours in Figure 9, where we also plot the boundaries of the corresponding zones for the power-law fluid case = 0 by dashed lines for comparison purposes. The general solution may approach the leak-off-toughness (̃ ), storage-viscosity ( ), leak-off-viscosity (̃ ), and/or storage-yield-stress ( ) regimes at an intermediate time depending on the values of and . The behaviour of the yield stress dominated regimes with an increase of is qualitatively similar to the already discussed case of the Bingham fluid, namely, and̃ domains expand towards smaller times. This shift leads to shrinking and, eventually, vanishing time domains of the viscosity dominated regimes as it can be noticed for = 0.3 from Figures 11(a), (b) as compared to the panels (c), (d), and to the shrinking of the toughness regimes time domains (see Figure 11(d)). In contrast to > 0.5 case, for < 0.5 the change of the dimensionless yield stress does not affect the boundaries framing the ̃ transition. We also should comment on the alterations of the regime map with the variation of the flow behaviour index inside the interval ∈ [0, 0.5) (not shown). When we increase the value of , the boundaries of and vertex solutions framing the transition move to the left (i.e., smaller times) and to the right (i.e., larger times) respectively, whereas the boundaries of̃ and̃ domains framing thẽ ̃ transition go up-left and down-right, correspondingly. Similarly to > 0.5, the flow behaviour index affect the boundaries of the yield stress dominated regimes framing and̃ ̃ transitions only for the large values of the dimensionless yield stress. Figure 12 illustrates the evolution of fracture characteristics with time for the Herschel-Bulkley fluid ( = 0.3, = 1) for various values of the leak-off number: = 10 −20 , 10 −10 , 1, 10 10 , 10 25 (see the corresponding dashdotted trajectories in the parametric map in Figure 11(c)). The solutions for the power-law fluid case ( = 0.3, = 0) are also shown for comparison by grey dashed lines. (Note that the solutions for the radius, maximum opening, and pressure at the half-radius in panels (a)-(c) in Figure 12 are normalised by the storage-viscosity limiting solution ( ).) One can notice that the fluid yield stress impacts the problem solution qualitatively similar to the Bingham fluid case: (i) for the fixed yield stress number, the increase of leak-off results in the fracture with smaller radius, maximum opening, and efficiency, but higher fluid pressure; (ii) when we fix the leak-off number and raise the yield stress number, we observe the reduction of the radius and increase of opening at the wellbore, pressure, and efficiency. The tendencies (i) and (ii) are applicable for the large enough time when the solutions for Herschel-Bulkley fluid ( > 0) and Newtonian

Quantitative estimations of the plug zone
In the current model for the hydraulic fracture driven by Herschel-Bulkley fluid, the plug zone is formed inside the regions where the shear stress is less than the yield stress 0 , and its width is equal to plug ( , ) = 2 = Figure 13: The charts present the isolines for Υ( , , ) = plug ∕ crack (the dashed black lines) corresponding to the quantities 10%, 50%, and 90%. Two cases of the flow behaviour index are analysed: = 1 (panels (a) and (c)) and = 0.3 (panels (b) and (d)). The top row reflects the impermeable reservoir case, while the bottom figures correspond to simulations with leak-off and for the yield stress number = 1. The regime maps demonstrated in Figures 7, 8, 9, 11 are applied as background in the diagrams.
The fully numerical solution (Section 3.1) is utilised for calculating Υ( , , ) time histories for various , since the approximate approach (Section 3.2) provides the pressure gradient profiles with a reduced accuracy. The results are shown in Figure 13 in the form of the isolines (the dashed black lines) Υ( , , ) = const where the constant can be varied within the segment [0, 1] and, in our case, takes values 0.1, 0.5, and 0.9. The plug fraction isolines are shown for two values of the flow behaviour index = 1 (the left column) and = 0.3 (the right column) for both impermeable (the top row) and permeable (the bottom row) formation cases. We demonstrate the results in the parameter space ( , ) when = 0 and in the space ( , ) for the fixed yield stress of = 1 for non-zero leak-off > 0. To better understand positions of the isolines relatively to the limiting regimes, we add their validity zones using the same colour palette as in the regime maps in Figures 7, 8, 9, 11. The selected isolines Υ( , , ) = const are predominantly located in the transition regions between the limiting propagation regimes. For a crack driven by Bingham fluid in an impermeable rock (panel (a) in Figure 13), and transitions contain the contour lines within the time intervals ≲ 10 −2 and ≳ 10 6 , respectively. The reversed situation is observed for the Herschel-Bulkley fluid with = 0.3 (panel (b) in Figure 13) in which the isolines are located inside the and transition zones for ≲ 10 −18 and ≳ 10 7 , correspondingly. The isoline Υ = 0.1 passes closely to the applicability domains of the storage limiting regimes, i.e., viscosity and toughness , in which Υ = 0. In turn, the level Υ = 0.9 lies near the validity region of the vertex associated with Υ = 1.
Generally, the function Υ( , , ) for the fixed values of the governing parameters and grows smoothly with time from 0 to 1 evolving from the viscosity ,̃ or toughness ,̃ limiting regimes to the yield-stress ,̃ vertices. In the model with the non-zero leak-off (panels (c) and (d) in Figure 13), we observe that the isolines belong to ( ) and̃ ̃ (̃ ̃ ) transitions when the leak-off number ≪ 1 and ≫ 1, respectively, and = 1 ( = 0.3). Moreover, one can notice that the contour lines are parallel to the boundaries of the applicability domains of the vertex solutions framing the transition zones enclosing the isolines.
Finally, we focus on the behaviour of the isolines near the toughness dominated limiting regimes when the leak-off number is fixed, and the yield stress number is small (for > 0.5) or large (for < 0.5). Using the panels (a) and (b) in Figure 13 for = 0 cases, it is possible to observe that the contour lines for Υ = 0.1, 0.5 penetrate into the validity zone of the storage-toughness regime. This phenomenon can be explained in the following way: the unyielded region starts to form when the crack geometry (radius and aperture) corresponds to the vertex solution in which the energy is spent primarily on the brittle rock failure and is independent of the fluid flow process inside the crack channel where the coexistence of both liquid and solid (plug) states can occur.

Simulations of the crack growth for typical field cases
This subsection outlines results of the radial crack propagation in terms of the dimensional variables (the system of the governing equations can be found in Section 2.2). We choose the values of the input parameters representative of typical field applications, and they consist of the geomechanical and filtration-storage properties of the porous rock, fluids (pore and hydraulic fracturing) characteristics, and the injection rate. The main aim of the analysis is to examine the impact of non-zero yield stress on the problem solution quantitatively for particular cases relevant to the field.
We consider fracture propagation during the first 6 ⋅ 10 3 seconds of injection with the volumetric rate of 0 = 0.01 m 3 ∕s. The set of the geomechanical parameters is: plane-strain elastic modulus ′ = 30 GPa, rock toughness = 1 MPa ⋅ √ m, far-field confining stress = 10 MPa. We take the fracturing fluid with the flow behaviour index = 0.7, consistency index = 1 Pa ⋅ s n providing ′ = 7.7 Pa ⋅ s n , and the yield stress value is 0 = 10 Pa (since we estimate the yield stress influence, the reference solution corresponds to the power-law fluid with the same and but 0 = 0). We analyse the crack propagation in both impermeable and permeable formations. In the latter case, we consider the pore fluid and formation characteristics to be the following: far-field pore pressure = 6 MPa, permeability = 10 mD, porosity = 20 %, viscosity = 5 cP (Newtonian fluid), compressibility = 10 −4 atm −1 . For simplicity, it is assumed that the filtrate of the fracturing fluid has the same properties as the pore fluid. Consequently, the leak-off parameter is equal to ′ = 9.1 ⋅ 10 −5 m∕ √ s in the permeable reservoir case. In terms of the dimensionless governing parameters ( , , ) introduced in Section 4.2, the analysed cases can be written as: ∈ (0, 1.04 ⋅ 10 −10 ), = 0 (impermeable rock) or = 4.7 ⋅ 10 14 (permeable rock), = 5.8 ⋅ 10 4 (Herschel-Bulkley fluid) or = 0 (power-law fluid). Figure 14 shows evolution in time of the fracture radius ( ), opening at the wellbore (0, ), net pressure at the half-radius ( ( )∕2, ), and efficiency ( ) calculated using fully numerical solution (see Section 3.1). The results for a radial crack driven by the Herschel-Bulkley fluid are depicted by the solid lines, for a crack driven by power-law fluid by the dashed lines. We utilise the blue colour for the zero leak-off case and the green one for when the rock formation is permeable.
We observe that the non-zero yield stress 0 = 10 Pa leads to the reduction of the crack radius, increase of the opening at the wellbore, pressure, and efficiency (there is no alteration in efficiency for the impermeable formation case) (see the solid lines versus the dashed ones of the same colour in Figure 14). Therefore, these findings are in the agreement with the analyses stated in Section 5.2. To characterise the discrepancies between two compared solutions, i.e., for Herschel-Bulkley ("hb") and power-law ("pl") fluids, quantitatively, we calculate the relative differences between the crack properties for the time moment = 6000 s using the following formula: = | hb − pl |∕ hb , where is an analysed parameter. The obtained metrics are summarised in Table 2. One can notice that the variations between the two compared solutions are more substantial during the fracture propagation in an impermeable rock. Here, for the radius and maximum aperture parameters is approximately 9 % and 20%, correspondingly. When leak-off is introduced into the model, the radius values are almost indistinguishable, while the opening at the wellbore and crack efficiency are higher for Herschel-Bulkley fluid by 11.4 % and 8.4 %, respectively. Finally, for the net pressure at the half-radius reaches approximately 27 % and 13 % for = 0 and > 0. All these clearly demonstrate the importance of yield stress for modelling practical cases of hydraulic fracturing.

Conclusions
In this paper, we constructed a numerical model for a radial hydraulic fracture driven by Herschel-Bulkley fluid in a permeable rock and carried out extensive analysis of the combined impact of the fluid yield stress and leak-off on the

Table 2
The relative differences ( ) between various radial fracture parameters (A) corresponding to the fracturing fluids with the Herschel-Bulkley ("hb") and power-law ("pl") rheologies that are calculated at = 6000 s. = | hb − pl |∕ hb (in %) Parameter (A) = 0 = 5.8 ⋅ 10 4 ( ) 9.4 0.8 (0, ) 19.9 11.4 ( ( )∕2, ) 27.2 12.6 ( ) 0 8.4 crack propagation. In most of the existing models, the fracturing fluid is taken as Newtonian or power-law; however, certain fracturing fluids, such as gels and foams, demonstrate shear-thinning rheology with the yield stress that can be described by the Herschel-Bulkley model. The system of governing equations is formulated in the dimensional and normalised forms, and, for the latter, all crack characteristics are governed by two dimensionless parameters: the leak-off number and the yield stress number . Two different approaches are used to obtain the problem solution: the accurate method based on Gauss-Chebyshev quadrature and Barycentric Lagrange interpolation techniques and the simplified algorithm that utilises the global fluid balance equation together with spatial continuation of the neartip asymptotic model. We examine the limiting propagation regimes characterised by certain energy dissipation and fluid balance mechanisms. Further, we investigate the problem parametric space using the dimensionless formulation by looking at the validity domains of the vertices and variations of different time-dependent crack characteristics versus the dimensionless problem parameters. The regime maps allow us to rapidly detect the dominance of various physical processes in a hydraulic fracture for certain values of the governing parameters without running numerical simulations. Using the obtained results, we conclude that the non-zero yield stress leads to a radial crack with a smaller radius, larger crack opening at the wellbore, pressure, and efficiency (for a permeable rock case) as compared to the power-law fluid case. Next, we simulate the radial fracture growth by taking the input parameters close to typical field applications and reveal that the yield stress can potentially result in notable deviations of the fracture parameters from the outcomes of zero yield stress model. The proposed radial crack model is applicable not only for computation of fracture parameters for the radial hydraulic fracture driven by fluid with the non-zero yield stress, but also for validating numerical calculations of more complex numerical simulators, such as a three-dimensional planar fracture (Planar3D) model, in which the Herschel-Bulkley rheology is also employed.

A. Calculation of for the approximate numerical solution (Section 3.2)
In Section 4.1, we outline the values of for all limiting propagation regimes realised in the model: = ( ) (see Figure 2(a)), ̃ = ̃ ( ) (see Figure 4(a)), = ̃ = 0.5, = 1.029, and ̃ = 1.077. As it was demonstrated in Section 5.2, the validity regions of the vertex solutions strongly depend on time and the governing parameters. Consequently, based on the values of for the vertices, it is necessary to construct an interpolation scheme allowing one to determine for all possible values of input parameters.
Let us define the auxiliary set of parameters: where the subscript denotes the considered limiting regime, and it can take values = ,̃ , ,̃ , ,̃ . The radius ( ) and wellbore opening (0, ) evolution in the vertex solutions are defined by the equations (34), and the function ( ,̄ ) given by equation (27). The parameter pl can be interpreted as the approximate efficiency for the crack filled by the power-law fluid, and the similar interpretation can be done for 0 which is more suitable for the yield-stress dominated regimes ,̃ .
Using the variables (61), we can construct the following interpolation function for : where is a fitting parameter which is the current study is set 3. We also should mention that for the zero leak-off case Ξ̃ = Ξ̃ = Ξ̃ = 0. Similarly, when we simulate the radial crack driven by the power-law fluid, we set Ξ = Ξ̃ = 0.

C. Validation of the approximate solution
This section presents the validation procedure for the constructed approximate numerical solution (Section 3.2) with the help of the accurate approach (Section 3.1). We check its accuracy in terms of the radius, opening at the wellbore, pressure at the half-radius, and efficiency time histories carrying out simulations for different values of the governing parameters , and the flow behaviour index . In Section 4.1, we have already looked at the accuracy of the approximate method in the limiting propagation regimes, and here, the main focus is aimed at the transitions between them.
For the impermeable formation case, we analyse the set = {1, 0.75, 0.3} and for each value of the flow behaviour index, the computations are performed for = 10 −10 , 10 −5 , 10 0 , 10 5 . Figure 15 shows the crack properties related to = 1, and they are normalised by the storage-viscosity limiting solution. We estimate the relative differences between the fracture characteristics in the fully numerical and approximate solutions according to the relation: = | numer − appr |∕| numer |, = { , Ω(0), Π(0.5)}, find the maximum and average quantities of across the time domain [ min , max ]. We summarise the calculations in Table 3 together with min and max values. Using the metrics from Table 3 in combination with Figure 15, we conclude that the approximate methodology provides reasonably accurate predictions of the radius evolution. The average error in the wellbore opening calculations is relatively small (the peak value among the considered cases is 2%); however, the maximum difference can reach 16% in the transition zone between and vertex solutions. The pressure profiles are computed with the average error lying within the interval (2, 11)%, and there is a tendency of the difference growth with an increase of the yield stress number. Similarly to the wellbore opening, the maximum error for Π(0.5) is reached between the viscosity and yield-stress regimes, and it is up to 16%. Note that max is not constant for all these cases, and it is linked with the ability of the fully numerical solution to perform calculations inside the yield-stress dominated regime where almost the whole fracture is occupied by the plug zone, and the fluid flow velocity approaches zero.
Let us move on to the validation of the approximate solution for crack propagation in a permeable rock. Both power-law ( = 0) and Herschel-Bulkley ( = 1) rheologies are analysed for = {1, 0.75, 0.3}, and we apply the following set of the leak-off number values: = 10 −20 , 10 −10 , 10 −5 , 10 0 , 10 5 , 10 10 . Figure 16 shows the comparison for = 1, = 1 (the radius, opening, and pressure profiles are normalised by the vertex solution), while the relative Figure 15: The comparison between the fully numerical (the solid black lines) and approximate (the dashed grey lines) numerical solutions corresponding to = 1 flow behaviour index, an impermeable rock ( = 0), and the set of the yield stress number values: = 10 −10 , 10 −5 , 10 0 , 10 5 . Three time-dependent crack parameters are considered: radius (a), opening at the wellbore (b), and pressure at the half-radius (c), and all properties are normalised by the storage-viscosity limiting solution.

Table 3
Table contains the values (in %) of the relative differences between the fully numerical and approximate numerical solutions in terms of radius , aperture at the wellbore Ω(0), and pressure at the half-radius Π(0.5). Two values of are written in each cell corresponding to the aforementioned parameters: maximum and average across the considered time span ∈ [ min , max ]. We analyse solutions for = {0.3, 0.75, 1} and the set of the yield stress number values = {10 −10 , 10 −5 , 10 0 , 10 5 }.
Further, we turn to the case = 1 (see metrics in Table 4 and Figure 16 for = 1), and here, the main trends Figure 16: The comparison between the fully numerical (the solid black lines) and approximate (the dashed grey lines) numerical solutions corresponding to = 1 flow behaviour index, the set of the leak-off number values = 10 −10 , 10 −5 , 10 0 , 10 5 , 10 10 , and = 1 yield stress number. Four time-dependent crack parameters are considered: radius (a), opening at the wellbore (b), pressure at the half-radius (c), and efficiency (d). The solution profiles in panels (a)-(c) are normalised by the storage-viscosity limiting solution.
are similar to the analysis for the power-law fluid. The average errors of the radius and efficiency computations are smaller than 1.5%, and the maximum differences reach 3% and 6%, correspondingly. The approximate approach gives the wellbore opening with the average accuracy less than 4%, while the peak differences are noticed betweeñ and regimes and comprise values up to 16%. The average precision of the pressure calculations is approximately 7%, and the difference between the fully numerical and approximate solutions for this parameter can achieve 13%. Table 4 Tables contain the values (in %) of the relative differences between the fully numerical and approximate solutions in terms of radius , aperture at the wellbore Ω(0), pressure at the half-radius Π(0.5), and efficiency . Two values of are written in each cell corresponding to the aforementioned parameters: maximum and average across the considered time span ∈ [ min , max ]. We analyse solutions for = {0.3, 0.75, 1}, and for each value of the flow behaviour index, both power-law ( = 0, table (a)) and Herschel-Bulkley ( = 10 0 , table (b)) rheologies are examined. In turn, the set of the leak-off number values = {10 −20 , 10 −10 , 10 −5 , 10 0 , 10 5 , 10 10 } is applied for calculations for each and .