Turbulent flow effects in hydraulic fracture propagation in permeable rock

This chapter considers a model for a radial hydraulic fracture propagation in a permeable, linear elastic rock formation driven by a point source fluid injection. The linear elastic fracture mechanics theory controls the quasi-static propagation. The hydraulic fracturing fluid is slickwater – pure water solution with polymeric additives which allow reducing the fluid flow friction in the wellbore and fracture in reservoir field applications. We focus on the possible transformation of the fluid flow regime inside the fracture channel from laminar to turbulent with distance from the fracture front. We assume that the turbulent friction of slickwater is described by the maximum drag reduction asymptote, while Carter’s law governs the leak-off into the permeable rock. The solution is obtained numerically using the algorithm based on the Gauss-Chebyshev quadrature and Barycentric Lagrange interpolation techniques. We compute solution examples for typical field cases and demonstrate a significant impact of the turbulent flow regime during the initial fewminutes of propagation, namely, shorter radius andwider maximum aperture than the laminar model provides. Moreover, we observe higher fluid pressure values at the wellbore within tens ofminutes of the start of the injection. This leads to a larger hydraulic pumping power requirement than the laminar model predicts. We also find that the fluid leak-off into the permeable rock enhances the turbulent flow effect in the fracture when compared to the impermeable rock case. In order to analyze the parametric dependence of the general solution, we convert the governing equations into the dimensionless form. We perform an extensive exploration of the normalized solution in space of two non-dimensional parameters, leak-off and characteristic Reynolds numbers, and normalized time. Specifically, we determine the applicability domains of the limiting propagation regimes to frame the general solution, investigate the alterations of the crack characteristics depending on the governing parameters, and identify zones where the turbulent flow is important.


Introduction
Hydraulic fractures (HFs) are tensile cracks that form and propagate in the solid media with the preexisting confining stress due to the high-pressure fluid injection. HFs manifest in nature as magmadriven dykes (Spence and Turcotte, 1985;Lister, 1990;Rivalta et al., 2015;Dontsov, 2016b) and waterdriven cracks in glacier beds (van der Veen, 2007;Tsai and Rice, 2010). They can also be humanmade, when hydraulic fracturing treatment of oil and gas wells is utilized to enhance production of hydrocarbons (Economides et al., 1989(Economides et al., , 2002. HFs ensure a larger hydraulically connected area between the rock and the well. In recent times, the horizontal drilling with the multistage hydraulic fracturing is frequently applied, especially in the development of formations with low permeability and porosity. When the water-based fracturing fluid is used in reservoir stimulation, the injection rate can be set relatively high. It is required to compensate for the undesirable effect of the relatively low viscosity of the HF fluid on the proppant settling and in order to create a crack aperture sufficiently large for the proppant placement (Barati and Liang, 2014). High injection rate of low-viscous fluid leads to the onset of the turbulent flow regime in the part of the fracture adjacent to the wellbore. The remaining part of the fracture, i.e., the region from the transition boundary to the fracture front, is occupied by laminar flow (Figure 1). Since Reynolds number for the plane channel flow depends on the aperture, the laminar flow regime always exists near the fracture front where the aperture value tends to zero.
Turbulent flow effects were first discussed in the context of hydraulic fracture propagation by Perkins et al. (1961) for the Perkins-Kern-Nordgren (PKN) fracture geometry. Further, Nilson (1981) analyzed a model for a plane strain gas-driven fracture, and this work discusses the flow regime transformation inside the crack channel from turbulent to laminar in the direction of the fracture front. The limiting propagation regimes associated with turbulent flow inside magma dykes were studied by Emerman et al. (1986), Lister and Kerr (1991). Tsai and Rice (2010) considered a plane strain HF with fully turbulent flow in application to natural water-driven fracture in glaciers. Different aspects of the turbulent flow regime impact on the HF growth with PKN geometry have been investigated in the following works (Anthonyrajah et al., 2013;Ames et al., 2015;Kano et al., 2015;Zia and Lecampion, 2017;Zolfaghari et al., 2017). Similar analyses have also been carried out for fractures with different geometries including a plane strain (Tsai and Rice, 2012;, radial Bunger, 2018, 2019;Lecampion and Zia, 2019), planar 3D (Dontsov and Peirce, 2017a) fracture. The aforementioned works consider finite crack models; however, the laminar-to-turbulent flow regime transformation has also been studied in the asymptotic near-tip region model. E.g., Dontsov (2016c) developed a model for a HF tip for the pure water case, while the slickwater case was examined by Lecampion and Zia (2019); Kanin et al. (2020a).
The fluid injection into the reservoir at a relatively high volumetric flow rate requires significant energy consumption for fluid pumping (Yang et al., 2019), especially during the simultaneous growth of multiple cracks, i.e., multistage HF treatment. The limitation of the operational costs can be achieved by adding the specific polymers to pure water resulting in a mixture called slickwater. The used additives slightly increase the fluid viscosity and considerably decrease the friction, up to 70% compared to that of the pure water (Nieuwstadt et al., 2016). The experiments of Virk (1971Virk ( , 1975 demonstrate how the flow friction factor declines depending on the used polymer additive type and its concentration. It was also determined in this work that the friction reduction has a so-called maximum drag reduction (MDR) or Virk's asymptote, which can be achieved at a relatively small concentration of the polymers. This chapter considers a radial hydraulic fracture propagation in a permeable reservoir. The model accounts for the flow regime transformation from laminar to turbulent when Reynolds number for the plane channel flow becomes greater than the critical value. The HF fluid is slickwater, and the MDR asymptote is used to approximate its frictional behavior in the turbulent flow regime. The fluid exchange between the crack and permeable formation is assumed in the form of the leak-off process governed by the pressure-independent, one-dimensional Carter's law (Carter, 1957). The main aim of the present examination is to analyze the combined effects of the laminar-turbulent flow of the HF fluid with drag reduction agents and the fluid leak-off into the ambient rock. For the problem formulation, we utilize the original model framework of Lecampion and Zia (2019), where the authors Figure 1: A radial hydraulic fracture model with the laminar-turbulent channel flow and leak-off. discuss a radial HF driven by the laminar-turbulent flow of slickwater in impermeable rock.
We organize the chapter in the following way. Section 2 outlines the problem formulation and governing equations, while Section 3 introduces the numerical scheme. The solution examples corresponding to the typical field cases are considered in Section 4. Further, we investigate the limiting propagation regimes in Section 5, namely, revisit the existing limiting solutions and derive the new ones. Finally, Section 7 describes the results of the problem parameter space investigation, which is accomplished by using the problem formulation in dimensionless variables introduced in Section 6.

Problem definition
We examine a radial (penny-shaped) hydraulic fracture propagation from a fluid point source along the plane perpendicular to the far-field confining stress direction. The model is axisymmetric, i.e., all fracture properties are the functions of the distance to the source r and time t. The sketch of the model is presented in Figure 1.
We consider a constant volumetric rate Q 0 of fluid injection such that the injected volume is V inj (t) = Q 0 t. The rock is assumed to be a homogeneous, linear-elastic solid characterized by Young's modulus E and Poisson's ratio ν. We assume that the size of the region near the fracture tip, where the rock failure occurs, is small compared to the other lengthscales realized in the model (e.g., linked with the viscous fluid flow inside the fracture, leak-off into the formation, etc.), which is generally true when the in situ confining stress σ o is much smaller than the rock tensile strength (Garagash, 2019). Therefore, we can apply the linear elastic fracture mechanics (LEFM) theory to model the quasistatic crack propagation in the solid media with toughness K Ic (Rice, 1968). The fracture surface is loaded by the fluid pressure p f (r, t) from the internal side. The considered radial crack model is fully described by evolution in time of the fracture radius R(t), opening w(r, t) and net pressure p(r, t) = p f (r, t) − σ o profiles, and efficiency parameter η(t) = V crack (t)/V inj (t). The latter allows evaluating the partition of the injected fluid volume V inj (t) between that of the crack V crack (t) and the fluid volume leaked into the rock V leak−off (t).
The HF fluid is considered to be slickwater with viscosity µ and density ρ which behaves like a Newtonian fluid in the laminar flow regime. When turbulent flow occurs, slickwater rheological response is modeled by the maximum drag reduction (MDR) asymptote which we will discuss further. The lubrication theory (Batchelor, 1967) controls the fluid flow inside the fracture. We neglect the fluid lag (Garagash and Detournay, 2000;Detournay and Garagash, 2003;Kanin et al., 2020c), assuming that its maximum value ∼ µV E 2 σ −3 o expressed in terms of the crack tip velocity V , is small compared to the fracture radius. In other words, the crack and the HF fluid fronts coincide in the model. Reynolds number for the plane channel flow, Re = ρvw/µ, is increasing with distance from the crack front where Re = 0. The laminar-to-turbulent flow regime transformation happens at distance λ from the tip (this value is a part of the solution), where Reynolds number equals to the critical value, Re = Re c . Consequently, the laminar flow domain of extent λ is adjacent to the front, i.e., observed within the interval r ∈ (R − λ, R) along which Re < Re c , while the non-laminar flow, i.e., transient and fully turbulent, happens inside the domain r ∈ (0, R − λ) characterized by Re > Re c . We implement the flow regime transformation into the model via the usage of the fluid friction factor f (Re) as detailed further in this Section.
The ambient reservoir is characterized by porosity ϕ and permeability k. For simplicity, we utilize Carter's leak-off law (Carter, 1957) to govern the fluid exchange process between the crack channel and formation. According to this law, the leak-off rate g(r, t) is proportional to the inverted square root of the 'exposure' time, i.e., the period between the current time moment and when the crack tip passes the considered point on the fracture plane. We consider the case in which the pore fluid has the same properties as that of the slickwater (viscosity µ, the compressibility c t ). Here, the proportionality (or Carter's) coefficient is the following (Collins, 1976): where p o is the far-field pore pressure, and c = k/(ϕc t µ) is the diffusivity coefficient. When it is required to take into account the filter-cake building or (and) the different properties of the HF and pore fluids in the pressure-independent approximation of the fluid exchange process, the reader can find the appropriate expressions for Carter's coefficient in (Economides et al., 1989).

Governing equations
In this section, the governing equations are formulated for the evolution of unknown crack radius R(t), opening w(r, t) and net pressure p(r, t) profiles. In doing so, we make use of the following set of material parameters in common with the laminar HF models: and the injection rate Q 0 . In equation (1), E ′ is the plane strain elastic modulus, K ′ and µ ′ are toughness and viscosity parameters, and C ′ is the leak-off parameter. Additional problem parameters related to turbulent flow will be introduced further in Section 2.2.2.

Crack elasticity
The elasticity equation links the net fluid pressure p(r, t) and aperture w(r, t) (Arin and Erdogan, 1971;Cleary and Wong, 1985;Savitski and Detournay, 2002): where integral kernel M (ρ, s) is defined as: and K(x) and E(x) are the complete elliptic integrals of the first and the second kind, respectively.

Fluid flow
Fluid flow inside the crack channel is described by width-averaged continuity equation for an incompressible fluid (Batchelor, 1967): where q(r, t) = w(r, t)v(r, t) is the flow rate, v(r, t) is the width-averaged fluid flow velocity, and g(r, t) is the fluid leak-off rate given by Carter's law. The latter is defined in terms of the time t 0 (r) when the fracture tip reaches distance r from the source and can be obtained as the inverse of the fracture radius function, i.e., R(t 0 (r)) = r.
We determine v(r, t) from the width-averaged momentum conservation equation in which the inertial terms are neglected (see Lecampion and Zia (2019) for details): where f is the Fanning friction factor. By comparing equation (5) with Poiseuille's law, v = −(w 2 /µ ′ ) (∂p/∂r), we obtain the friction factor for laminar flow (Re < Re c ) in plane channel: f lam = 12/Re. Using the definition of the normalized friction factorf = f /f lam , we rewrite equation (5): The combination of the continuity (4) and momentum conservation (6) yields Reynolds equation: In the turbulent flow regime (Re > Re c ), the slickwater frictional behavior is governed by the MDR asymptote that is a phenomenological relation proposed by Virk (1971Virk ( , 1975 on the basis of the experiments: where Re d = ρυd/µ is Reynolds number for pipe flow (d is a pipe diameter). In order to 'translate' (8) to the 'crack channel' flow geometry, the relation between apparent pipe Re d and channel Re is postulated by equating corresponding expressions for the laminar friction factors for a channel, f lam = 12/Re, and a pipe, f lam d = 16/Re d , yielding Re d = 4 3 Re (Jones Jr, 1976). Finally, the implicit equation (8) is approximated by explicit power-law relation suggested by Lecampion and Zia (2019): where f ′ 0 = f 0 (4/3) −n , and the numerical values for f 0 and n will be provided later. Behavior of the friction f and normalized frictionf = f /f lam in the entire range of Reynolds number values is approximated by piecewise functions: Figure 2: The friction factor dependence on the Reynolds number in ordinary (a) and Prandtl-Karman (b) coordinates for channel flow. The laminar branch is shown by blue color, while the turbulent branches for the slickwater flow described by the original Virk's MDR asymptote (8) (Virk, 1975) and its power-law approximation (9) (Lecampion and Zia, 2019) are depicted by solid green and red lines, respectively. Blasius correlation for turbulent flow of Newtonian fluid (pure water in our case) in a channel with smooth walls is plotted by solid yellow line.
with prefactor f ′′ 0 = f ′ 0 /12. We estimate the critical value of Reynolds number Re c for the slickwater case from the intersection of the laminar and turbulent branches, to ensure the continuity of f (Re): . To summarize the discussion on the slickwater frictional behavior, we plot Figure 2. Here, we show the friction factor f dependence on the Reynolds number Re for the channel flow during the laminar and turbulent flow regimes in the traditional (a) and Prandtl-Karman (b) notations. Function f (Re) defined by equation (10) is shown by the combination of blue (laminar part) and red (turbulent part) solid lines. In addition to Virk's asymptote (green line), we show Blasius correlation (yellow line) for pure water flow in a channel with smooth walls. Using the green and yellow curves, one can observe how the slickwater drag reduction agents drastically reduce the friction compared to that of water. The power-law form of the MDR asymptote (9) (red line) closely approximates the exact function (8) within the range Re ∈ (10 3 , 1.5 · 10 4 ) with the relative error of less than 5%.
Numerical values of the parameters in rheological relations (9), (10) for the slickwater (power-law approximation of MDR) and pure water (Blasius correlation) cases are taken after (Lecampion and Zia, 2019) and (Blasius, 1913), correspondingly: The critical Reynolds number for the channel flow of Newtonian fluid (pure water) is evaluated from the corresponding parameter for the pipe flow (Re d ) c = 2200: Re c = 2200 · 3/4 = 1650. One can notice that the laminar-to-turbulent transition for slickwater occurs at a smaller Reynolds number value than that for pure water.

Fracture propagation
According to the LEFM theory, the quasi-static fracture propagation is realized when the stress intensity factor matches the solid toughness: K I = K Ic . Alternatively, this condition can be expressed in terms of asymptotic behavior of the crack opening w(r, t) near the tip (Irvin, 1957): where elasticity E ′ and toughness K ′ parameters are defined in (1).

Boundary conditions
At the fracture inlet (r = 0), the volumetric flow rate is prescribed: while at the fracture tip (r = R(t)), the zero crack opening and the no-flow conditions are applied:

Global fluid volume balance
Global fluid balance follows from the integration of the continuity equation (4) with respect to time and distance from the source and application of the flux boundary conditions (13) and (14), and corresponds to the partition of the injected fluid volume (V inj ) between the crack (V crack ) and porous rock (V leak−off ).

Solution approach
The current section describes the numerical algorithm utilized in the solution of the radial hydraulic fracture propagation problem. The method is based on the Gauss-Chebyshev quadrature and Barycentric Lagrange interpolation techniques. Viesca and Garagash (2018) review these techniques and apply them to the solution of various fracture propagation problems with coupled physics. Liu et al. (2019) adapt the methodology to the radial fracture propagation. It is important to mention that the numerical scheme does not require an explicit implementation of the near tip region asymptote (if different from the classical LEFM one) compared to other approaches which use the 'tip logic', e.g., (Peirce and Detournay, 2008;Peirce, 2015;Dontsov, 2016a;Dontsov and Peirce, 2017b;Zia and Lecampion, 2020;Kanin et al., 2020b). In present work, we follow the method developed by Liu et al. (2019) and extend it to account for different HF fluid rheology and for the fluid leak-off.
Firstly, we introduce the dimensionless distance from the source as ξ = r/R(t), ξ ∈ [0, 1], and rewrite the system of governing equations through ξ variable using the transformation of the time and coordinate derivatives: The application of Gauss-Chebyshev quadrature implies the problem consideration on the spatial coordinate segment ξ ∈ [−1, 1]. In this regard, the spatiotemporal fracture characteristics and governing equations are extended symmetrically to negative ξ. Such modified form of the elasticity equation (2) is given by: where the corresponding form of the integral kernel is provided by Liu et al. (2019): Spatial derivative of the opening profile is written as: where the singular weight function W(ξ) is chosen in accordance with the fracture opening behavior near the tip, i.e., ∂w/∂ξ ∼ 1/ 1 − |ξ|, ξ → ±1, and F (ξ, t) is an unknown regular function. Since w(ξ, t) and W(ξ) are even functions, function F (ξ, t) should be odd.
Next, the coordinate domain ξ ∈ [−1, 1] is discretized using two sets of nodes connected with the selected weight function (the details are provided by Viesca and Garagash (2018)): primary Hereafter we utilize bold symbols for vectors. Primary and complementary nodes are the roots of the Chebyshev's polynomials of the first and the second kind. We further make use of the values of F (ξ, t) at the primary nodes: Let us turn to the discretization of the governing equations. We focus on the Reynolds and global fluid balance equations since the remaining equations (elasticity, propagation criterion, and boundary conditions) are identical to those in (Liu et al., 2019) and given by their equations (3), (10), (8). Integrating Reynolds equation (7) in space from a given position along the crack ξ to the tip ξ = 1, applying discretization and substituting discretized elasticity (equation (41) Here, we utilize the sign '×' for the matrix multiplication. Components of the matrices S, D, G 2 are provided by Liu et al. (2019) in their equations (34), (39), (43), while matrix R is given by: The rate form of the global fluid balance equation (15) can be discretized as follows: where the matrix S H is defined by Liu et al. (2019) in their equation (38), while R H is given by: In the right-hand-side of the global fluid balance equation we have the integral with respect to time, and it can be evaluated via application of the Simpson's rule.
The discretized Reynolds equations (with the elasticity equation and boundary conditions already taken into account in it) together with the global fluid balance and propagation condition can be combined into a system of ordinary differential equations (ODEs): where X is the solution vector. Once the above system of ordinary differential equations is integrated, the crack radius evolution is known, R = R(t), while the solution for opening, net pressure profiles, and efficiency is evaluated from F (ξ, t) using the following representation: where matrix S A is defined in (Liu et al., 2019) by their equation (37).
Due to the form of the elasticity equation (Liu et al., 2019), the number of primary nodes N should be odd. Since function F (ξ, t) is odd, the target vector F has the form: The total number of the independent unknown parameters is (N − 1)/2 + 1 which constrains the system of governing ODEs to be composed of the Reynolds equations at the complementary nodes z 1 , . . . , z (N −1)/2−1 , the global fluid balance equation and the propagation condition. We choose N = 101, and the initial condition is represented by the storage-viscosity-turbulent limiting propagation regime originally developed by Lecampion and Zia (2019), and further discussed in Section 5. The time computational domain is discretized on a logarithmic scale, and we apply ODE solver for each time step. The numerical algorithm is implemented in Python programming language, the system of ODEs is solved via solve_ivp function of SciPy library.

Solution examples for typical field applications
In this section, we demonstrate the results of the simulations for a radial hydraulic fracture growth in which the model parameters are close to the typical field cases. We would like to achieve the following objectives. First of all, we investigate how the laminar-to-turbulent flow regime transformation inside the fracture channel changes the problem solution compared to the fully laminar flow case for both impermeable and permeable reservoirs. Moreover, we perform calculations not only for slickwater but also for pure water driven cracks to examine the influence of the frictional behavior on the fracture parameters.
The injection volumetric flow rate is assumed to be a fairly high Q 0 = 0.1 m 3 /s (e.g, see representative ranges in the reviews (Barati and Liang, 2014;Barbati et al., 2016)). For simplicity, we consider a situation when the polymer molecules do not change the solvent (pure water) viscosity, i.e., they modify the frictional behavior during turbulent flow only (the MDR asymptote instead of the Blasius approximation). The following rock and fluid properties are chosen: Using the set of parameters (18), one can estimate the Carter's leak-off number corresponding to the permeable rock case, C ′ = 2.04 · 10 −4 m/ √ s, while C ′ = 0 for impermeable rock. Using these figures, one can notice that the laminar-to-turbulent flow regime transformation inside the crack channel affects the fracture parameters during an initial period of the propagation, leading to a shorter radius, larger opening at the wellbore, and pressure values compared to the fully-laminar model. Turbulent effects eventually become negligible, and the turbulent-laminar and fully-laminar solutions become nearly identical. It is important to highlight that the threshold time is different for radius, opening, and pressure for a particular choice of the governing parameters of the model.
Let us define a time moment t * when a time-dependent fracture characteristic A(t) (radius, opening, pressure) corresponding to the turbulent-laminar solution can be closely approximated by A lam (t) from the fully-laminar solution as follows: |A(t * ) − A lam (t * )|/|A lam (t * )| = 5%. We summarize the values t * for the crack parameters depicted in Figure 3 for all considered cases, namely, slickwater  It is evident from Table 1 that the noticeable difference in the fracture radius values R(t) and R lam (t) appears on the timescale less than one second, suggesting that in practice, the deviation of the fracture size in the turbulent-laminar solution from that assuming fully-laminar flow can be neglected. The flow regime transformation impacts opening at the wellbore during approximately 5 and 15 seconds for slickwater fracturing in an impermeable and permeable rock, correspondingly. These quantities for pure water fracturing are around 1 and 12 minutes, respectively. The interpretation of the results for w(0, t) leads to the following conclusions: (i) the leak-off process prolongs the duration of the turbulent flow regime influence on the crack aperture near the wellbore, and (ii) the turbulence effects are more pronounced for the pure water case since they continue during much longer period of time. However, even when the crack geometries (radius and aperture) are approximately the same in the turbulent-laminar and fully-laminar models, the pressure values near the wellbore (r = 0.1 m) can still differ significantly. This observation means that larger amount of energy is required to create a hydraulic fracture driven by turbulent-laminar flow compared to the value predicted by the fully-laminar model. From Figure 3 and Table 1, we identify for the slickwater case that p ≈ p lam (less 5% difference) after 5 minutes from the fracture initiation in an impermeable formation, while in a permeable rock, this period reaches 2.5 hours. In turn, when the HF fluid is pure water, the alignment of the pressure values occurs after 0.5 and 18 hours (beyond typical HF treatment duration in the field), correspondingly.
Figure 3(d) shows the extent λ(t) of the region r ∈ (R − λ, R) with laminar flow as a fraction of the crack radius R(t). It is an increasing function on time evolving from λ/R ≈ 0 at early time and approaching λ/R ≈ 1 at large time. In other words, the turbulent flow regime is realized along the entire fracture at the beginning of the propagation, while the laminar regime is spatially dominant at large time. For example, we observe in the permeable rock case that near the end of a typical fluid injection (t ∼ 10 4 s) only a small part of the crack ∼ 0.1R near the wellbore is occupied by turbulent flow, whereas the remainder of the crack ∼ 0.9R supports laminar flow. Turbulent flow spatial extent is yet smaller for the case of impermeable rock. By comparing slickwater and water injection cases in Figure 3(d), we notice that the spatial domain with turbulent flow is larger for slickwater HF, i.e., λ/R at all times is smaller, yet the water HF corresponds to larger deviations of R(t), w(0, t), p(0.1, t) from the laminar solution.
We now discuss how the spatial distributions of the fracture aperture and net pressure along the fracture evolve with time. These distributions are shown in Figure 5 for three time moments t = {1, 10, 100} s. As it has been already mentioned, a radial crack driven by turbulent-laminar flow has a shorter radius compared to the fully-laminar model, at early time, the opening and pressure profiles for t = 1 s in Figure 5 confirm this fact. In the subsequent time moments, i.e., t = 10 and 100 s, the difference between R(t) and R lam (t) values is imperceptible. The turbulent-laminar solutions become very close to their laminar 'analogs' starting from a certain distance from the fluid source. In turn, in the remainder of the crack, i.e., along the zone adjacent to the wellbore, a considerable difference, decreasing over time, is observed where the turbulent-laminar crack has a wider opening and larger net pressure. In general, the HF in a permeable rock (Figures 5(b) and (d)) has a smaller volume due to the leak-off compared to the crack propagation in the impermeable reservoir, which is manifested in a shorter radius and aperture in the leak-off cases. On the other hand, the net pressure profile has greater values in the cases with leak-off.
Using Figures 5(c) and (d), one can observe the presence of the pressure singularities at the wellbore and the fracture front. For example, the pressure behavior near the wellbore (r ≪ R) is governed by p ∼ | log r| for the fully-laminar crack, while p ∼ r −3/10 and p ∼ r −3/4 are applicable for slickwater and pure water fracturing, respectively, in the turbulent-laminar model. The near-tip region of a hydraulic fracture (R − r ≪ R) has multiscale nature (see, e.g., (Garagash et al., 2011)) such that the dominant tip singularity depends on both fracture length and propagation speed. It can take form of either toughness asymptote p ∼ −| log r| or storage-viscosity asymptote (Garagash et al., 2011). The latter is given by p ∼ −(R − r) −1/3 (Desroches et al., 1994)

Limiting propagation regimes
Two different mechanisms control the propagation regime of a finite hydraulic fracture. The first one regulates the total dissipated energy distribution between the creation of new fracture surfaces at the tip and viscous fluid flow inside the fracture channel. The second mechanism is related to the partitioning of the injected fluid volume between the fracture and host permeable rock (due to leak-off). During the crack growth, the allocation of the dissipated energy and injected volume can change over time resulting in the realization of the limiting propagation regimes characterized by one dissipation (out of two) and one fluid balance (out of two) mechanisms. The partitioning of the dissipated energy is influenced by the viscosity µ ′ and toughness K ′ parameters, while the leak-off parameter C ′ affects the distribution of the injected volume.
The toughness dominated regimes K andK do not depend on the flow regime type. Therefore, the remaining two limiting regimes emerge in the viscosity-dominated case, when the entire fracture is occupied by turbulent flow: • T -storage-viscosity-turbulent -C ′ = K ′ = 0, Re > Re c ; •T -leak-off-viscosity-turbulent -C ′ → +∞, K ′ = 0,, Re > Re c .
Although the laminar flow always exists near the hydraulic fracture front, the above two regimes correspond to the limit when the laminar domain becomes negligibly small. Savitski and Detournay (2002) derived solutions for M and K vertices, whileK andM solutions are given by Bunger et al. (2005), and by (Madyarova, 2003;Peirce and Detournay, 2008), respectively. Moreover, Dontsov (2016a) present approximate solutions for all laminar limiting regimes. The semianalytical and approximate solutions for the storage-viscosity-turbulent regime, T -vertex, are found by Lecampion and Zia (2019).
Further, we summarize scalings in different limiting regimes of the solution following (Detournay, 2016) for the laminar regimes and (Lecampion and Zia, 2019) for the turbulent ones. Firstly, we introduce the dimensionless radius γ, opening Ω, net pressure Π and fluid velocity V as follows: where ξ = r/R(t) is the normalized distance from the source, P = {P 1 , P 2 , P 3 } are dimensionless evolution parameters depending on time, material parameters (1), fluid density ρ and injection rate Q 0 . L(t) is a lengthscale of the same order as the crack radius, and ϵ(t) is a small dimensionless parameter.
Next, we substitute formulas (19) into the governing equations and obtain their normalized form: • Elasticity: • Reynolds: • Propagation condition: • Boundary conditions: , Ω(1, t) = 0, Ω(1, t)V(1, t) = 0; • Global fluid balance: where τ 0 (t) = t 0 (r)/t is the dimensionless inverse radius function, ∇Π = ∂Π/∂ξ is the pressure gradient, and R = 12ρϵL 2 /(µ ′ t) is the characteristic Reynolds number. We also introduce five dimensionless groups: Numbers G v and G c quantify the fluid volume stored inside the fracture and leaked into the rock, respectively, in reference to the total injected volume. Numbers G k , G lam m , and G turb m with the meaning of the non-dimensional rock toughness, fluid viscosity, and equivalent viscosity of turbulent flow, correspondingly, relate to the partition of the dissipated energy between solid and fluid.
To derive scalings linked to the limiting regimes (i.e., lengthscale L(t) and factor ϵ(t) in equation (19)), we assign one out of two storage numbers (G v , G c ) and one out of three dissipation numbers (G lam m , G turb m , G k ) equal to one. The remaining three groups constitute the evolution parameters P. When all three of these parameters tend to zero in a given scaling, one can observe the emergence of the corresponding vertex solution.
We consider the storage-viscosity-turbulent scaling G v = G turb m = 1. The following formulas for the lengthscale and small parameter follow: while three non-dimensional evolution parameters P (i.e., toughness G k , leak-off G c , and laminarviscosity G lam m ) are: The solution in the storage-viscosity-turbulent regime (T -vertex) is then obtained from the normalized equations (20)-(24) in the above scaling when the evolution parameters (27) are set to zero. The solution is computed numerically using the method defined in Section 3, and is given below for the time evolution of the fracture radius, opening at the wellbore, and pressure at the half-radius: The storage-viscosity-laminar scaling G v = G lam m = 1 and corresponding scales L m and ϵ m can be alternatively obtained by setting n = f ′′ 0 = 1 in expressions (26) for L t , ϵ t . Fracture characteristics in the M -vertex solution (G k = G c = G turb m = 0) have the following form (Savitski and Detournay, 2002): where, as before, the numerical prefactors are evaluated from the full numerical solution of the problem.
Next, we look at the leak-off-viscosity-turbulent scaling G c = G turb m = 1 and derive the formulas for the lengthscale and small parameter provided below: Fracture characteristics in theT -vertex solution for the leak-off-viscosity-turbulent regime (G v = G k = G lam m = 0) can be obtained in the following form: where the numerical prefactors (for opening and pressure) are computed from the full numerical solution of the governing equations (20)-(24) in theT -scaling when the evolution parameters equal zero; the coefficient for the fracture radius is found analytically (exact).
The substitution n = f ′′ 0 = 1 in equations (30) for Lt, ϵt leads to the leak-off-viscosity-laminar scaling G c = G lam m = 1 and corresponding scales Lm and ϵm. The fracture properties in theM -vertex solution (G k = G v = G turb m = 0) are given below: where the numerical coefficients for the opening and pressure are calculated from the full numerical solution of the problem.
The toughness dominated vertex solutions (K andK) do not depend on the realized flow regime inside the fracture channel since the corresponding energy dissipation in the fluid flow is assumed negligibly small. For the completeness, we write out expressions for the radius, opening at the wellbore and pressure (uniform along the whole crack) in these regimes (Savitski and Detournay, 2002;Bunger et al., 2005): Normalized solution (19) in one of the vertex-scalings of Section 5 depends on three non-dimensional evolution parameters which quantify the departure from the corresponding vertex solution. Either of these three parameters can be regarded as non-dimensional 'time' in the problem. Alternative, so-called 'mixed' scalings have also been used (Madyarova, 2003;Bunger et al., 2005;Adachi and Detournay, 2008;Dontsov, 2016a) to compute solution evolution in time between the limiting regimes (vertices).
One of such scalings, the mk-scaling (Detournay, 2016), is defined in terms of the timescale t mk quantifying the solution transition period between M and K regimes. Specifically, t mk is the time moment when the K−vertex and M −vertex lengthscales are the same, i.e., L k (t) = L m (t) at t = t mk . The corresponding length L mk = L m (t mk ) = L k (t mk ) and non-dimensional small parameter ϵ mk = ϵ m (t mk ) = ϵ k (t mk ) are used to define the mk-scaling: where the characteristic scales evaluated as: We note that either vertex scalings (Section 5) or the alternative 'mixed' mk-scaling can be utilized to obtain general solution and its evolution in time. In the following, we apply the mk-scaling, since it allows for more direct interpretation of the normalized solution (since the scales L mk and ϵ mk are constants compared to the time-dependent vertex scales L(t) and ϵ(t)).

Problem parameter space analyses
Let us investigate the parametric space of the model for a radial crack driven by turbulent-laminar flow of slickwater in both impermeable and permeable reservoirs. We perform the analyses using the problem formulation in the dimensionless form (Section 6) where the parameter space is three dimensional with the axes: time τ , leak-off number ϕ, and characteristics Reynolds number R.
Firstly, we determine the applicability domains of the limiting propagation regimes (Section 5) and present them as the regime maps. Such analysis is useful to frame the general solution inside the parameter space resulting in better understanding of the propagation conditions as a function of time.
Similar to Dontsov (2016a), we utilize the following criterion for the determination of the validity zone of the considered limiting solution i: where the subscript i can be M,M , K,K, T,T . In other words, the measure of the relative difference between the general numerical solution and given limiting case i is taken to be less than 1%, (38), for the limiting solution to be considered a valid approximation, and the fracture to be said to propagate in the corresponding limiting regime.

Zero leak-off case (impermeable rock)
Lecampion and Zia (2019) has already considered a radial crack propagation in an impermeable formation; however, the authors did not include the regime map for the corresponding reduced parameteric space (τ, ϕ = 0, R), and we would like to fill this gap. Figure  Further, we would like to clarify how exactly we determine the applicability boundaries shown in Figure 6(a). We estimate the borders numerically via criterion (38) at sample points at small and large R and then fit the evaluated points by the appropriate analytical functions derived from the transition timescales between the limiting regimes. As an example, let us consider the validity zone of the T -vertex solution shown by orange color in Figure 6(a). It is bounded by two power-law functions of time which are straight lines in the log-log scale. The first boundary located in the region R < 10 4 relates to the T M -transition, i.e., between T and M vertices. We expect that it corresponds to t ∼ t mt , where the time-scale t mt is the solution of L t (t) = L m (t). This boundary expressed in the normalized time is τ ∼ t mt /t mk = R 9/4 , where the proportionality coefficient is found from the fitting procedure at small R. The functional dependence for the second curve is obtained from the T K transition analysis and has the following form: τ ∼ t tk /t mk = R (5n−5)/(2n−4) , where the numerical prefactor is found from a similar fitting procedure at large R. Both curves are extended until the intersection point, bounding the domain depicted in Figure 6(a). The full usage of criterion (38) in the entire parametric space of the solution yields a very similar region but with a smoother boundary (without sharp corner point evident in the domain approximation shown in Figure 6(a)).
We also display two time bounds in Figure 6(a). The first one τ 0 (R, ϕ = 0) (dashed black line) indicates the time moment when the spatial extent of the laminar flow domain is a small fraction λ = λ 0 = R/25 of the fracture radius. As a result, the crack can be approximately considered as fully-turbulent for τ < τ 0 (R, ϕ = 0) (and λ < λ 0 ). One can observe that τ 0 (R, ϕ = 0) is an increasing function of the characteristic Reynolds number R, i.e., the time interval, within which the length of the laminar flow region is small, and turbulent effects dominate, grows with an increase of R. For R < 10 3 , the dashed line is very close to the border of the T -vertex validity domain. The second bound τ ∞ (R, ϕ = 0) (dotted black line) corresponds to the time moment following which the crack geometry (radius and aperture) in the turbulent-laminar case is approximately the same as the fully-laminar model provides. In other words, one can interpret the fracture as fully-laminar for τ > τ ∞ (R, ϕ = 0). This bound is computed via equation (38) where i should be understood as a fully-laminar solution. The function τ ∞ (R, ϕ = 0) is an increasing one (similarly to τ 0 (R, ϕ = 0)), and for R < 10 2 , the dotted line coincides with the boundary of the M -vertex domain.
Next, we discuss the time bound τ 0 (R, ϕ = 0) in the context of the tip element concept. Many numerical models for hydraulic fracture growth are based on the 'tip logic' (Peirce and Detournay, 2008;Peirce, 2015;Dontsov, 2016aDontsov, , 2017Zia and Lecampion, 2020). It means that specialized neartip region model is applied within such algorithms to determine the fracture front position at each time step (propagation criterion) and to describe crack characteristics (aperture, pressure) near the tip. The fracture tip model should include all physical phenomena that can be realized during the fracture evolution and resolve their influence on the propagation. It is assumed that the near-tip region model is valid along small distance from the fracture front, and in the case of 1D models (plane strain or penny-shaped) the typical tip element length equals λ 0 . When we analyze the turbulent-tolaminar flow transformation within the crack channel, two different options for the near-tip region model exist. Firstly, when λ > λ 0 or for τ > τ 0 (R, ϕ = 0) laminar flow occurs along the entire tip element, one can apply the model of Garagash et al. (2011) or its approximate and computationally efficient version (Dontsov and Peirce, 2015). On the other hand, when λ < λ 0 or for τ < τ 0 (R, ϕ = 0), it is recommended to utilize the laminar-turbulent tip model developed by Lecampion and Zia (2019).  Figure 8, one can notice that the validity regions of the turbulent limiting propagation regimes, T andT , expand with an increase of R and gradually reduce the applicability zones of the laminar regimes, M andM , until their complete disappearance (see panels (c) and (d) in Figure 8). Similar behavior is observed for the toughness-dominated regimes, K andK, from the regime map in Figure 8(d); however, they do exist for all R and large enough τ .
In Section 7.1, we have already discussed the physical meaning of the temporal bounds τ 0 and τ ∞ , and here, we analyze their behavior with the alteration of the governing parameters R and ϕ (see Figure 8). For R = const, both timescales are independent of ϕ when it is small. However, for large values of the leak-off number ϕ, we determine that τ 0 , τ ∞ ∼ ϕ 2 . When the characteristic Reynolds number grows, the bounds shift towards larger time since the turbulent flow regime prevails inside the fracture channel during longer time period for greater values of R. We should also mention that it is recommended to apply the turbulent-laminar fracture tip model (  The dashed black lines illustrate the time τ 0 (R, ϕ) at which the flow regime transformation inside the crack channel occurs at small distance λ 0 = R/25 from the tip, i.e., the length of the laminar flow spatial domain is small λ < λ 0 for τ < τ 0 (R, ϕ). The dotted black lines shows the time τ ∞ (R, ϕ) starting from which the fracture radius and aperture is approximated by the fully-laminar solution.
In panel (b), the grey dash-dotted lines emphasize the considered solution trajectories discussed in the current section.
agation criterion in more complex numerical models for HF growth (such as Planar3D) to simulate the crack evolution when the laminar flow domain is small, corresponding to the parametric zone bounded by the dashed black line in Figure 8, i.e., when τ < τ 0 (R, ϕ). Figure 9 demonstrates the crack characteristics in the turbulent-laminar (R = 10 2 ) and fully-laminar (dashed grey lines) cases for different values of the leak-off number ϕ = {10 −20 , 10 −5 , 10 5 , 10 10 }.
Since the characteristic Reynolds number is constant, we focus on the impact of the leak-off on the propagation of a radial crack driven by turbulent-laminar flow. The storage-viscosity-laminar (M ) limiting solution is utilized to normalize the time-dependent properties such as the radius γ(τ ), opening at the wellbore Ω(0, τ ) and pressure at the half-radius Π(1/2, τ ) in panels (a) -(c) in Figure 9. It is evident from Figure 9(a) that larger values of leak-off number ϕ lead to shorter time duration over which the fracture radius differs significantly from the fully-laminar case (i.e., when turbulent effects on the fracture run-out are significant). The situation is opposite for the opening at the wellbore and pressure at the half-radius (see Figures 9(b) and (c)) for which increase of the leak-off extends the influence of the turbulent flow effects. Finally, the fracture efficiency is roughly independent of the laminar-to-turbulent flow regime transformation in all considered cases, as exposed in Figure  9(d).

Conclusions
In this work, we have analyzed how the turbulent-laminar flow inside the fracture channel impacts the propagation of a radial hydraulic fracture in a permeable formation. We account for the flow regime transformation from the laminar to turbulent moving away from the fracture tip. Both flow regimes always coexist in the model -the turbulent at the wellbore and the laminar one near the tip. The fluid exchange process between the fracture and ambient rock is modeled in the form of Carter's leak-off law. The hydraulic fracturing fluid is slickwater, i.e., the water-based fluid with polymeric additives increasing the solvent viscosity and significantly modifying the turbulent flow frictional behavior. The latter is governed by the phenomenological maximum drag reduction (Virk's) asymptote. We carry out several simulations of the radial crack growth for the values of the model parameters corresponding to typical field cases. The problem solution is found numerically via the Gauss-Chebyshev quadrature and barycentric interpolation techniques. The obtained results demonstrate that turbulence changes the crack characteristics near the wellbore zone during the initial propagation period. E.g., in the slickwater fracturing, the turbulent-laminar solution differs significantly from the fully-laminar one during a couple of seconds, tens of seconds, and several minutes in terms of the radius, opening at the wellbore, and pressure, respectively. We reveal that the leak-off process prolongs the turbulence effects, and they continue much longer for pure water fracturing than for the slickwater driven crack. Although the fracture geometries (radii and aperture profiles) in the models with and without turbulent flow are practically the same after several minutes of the propagation, the pressure near the wellbore has larger values in the turbulent-laminar case during tens of minutes, meaning a greater amount of power is required for the fluid pumping than the laminar model predicts. The latter confirms that slickwater HF propagation, especially in permeable rock, remains in the turbulent-laminar regime on the timescale of typical treatment duration.
We have analyzed the limiting propagation regimes, characterized by the dominance of a single energy dissipation and a single fluid storage mechanism, and derive corresponding limiting solutions semi-analytically. Then, we convert the system of governing equations into the dimensionless form in which the solution depends on two numbers, normalized leak-off ϕ = C ′4 E ′11 µ ′3 Q 0 /K ′14 and characteristic Reynolds R = 12K ′4 ρ/(E ′3 µ ′2 ), respectively, in addition to the dimensionless time and distance from the source (spatiotemporal properties). Using the normalized variables, we explore the problem parameter space by determining the applicability domains of the limiting regimes Figure 9: The time-dependent characteristics of a radial hydraulic fracture driven by turbulentlaminar flow of slickwater in a permeable rock (ϕ > 0): (a) radius γ(τ ), (b) opening at the wellbore Ω(0, τ ), (c) pressure at the half-radius Π(1/2, τ ), and (d) efficiency η(τ ). The properties in panels (a) -(c) are normalized by storage-viscosity-laminar (M ) limiting solution. The solution profiles corresponding to R = 10 2 and ϕ = {10 −20 , 10 −5 , 10 5 , 10 10 } are shown. The analogous profiles from the fully-laminar solutions are depicted by the grey dashed lines. We plot the vertex solutions by the colored dotted lines; in panel (a), the leak-off dominated regimes (K,M ,T ) have the same color since the fracture radii in these regimes are governed by the same relation (e.g., see equation (31)). and looking at the variation of the crack characteristics with the change of the values of the governing parameters. We introduce the temporal bounds, τ 0 (R, ϕ) and τ ∞ (R, ϕ), which can be used to determine the duration of the approximately fully-turbulent fracture propagation τ < τ 0 and approximately fully-laminar one for τ > τ ∞ . The proposed model can be used as a benchmark solution for the numerical simulators of more realistic (complex) fracturing problems which include the laminarto-turbulent flow regime transformation, e.g., turbulent-laminar Planar3D model.