4.1 Problem formulation

In order to better understand the influence of turbulent flow inside a hydraulic fracture, the aim of this section is to analyse the tip region of a hydraulic fracture in situations where the fluid flow can become turbulent. The tip region of a hydraulic fracture is known to obey a multiscale behaviour even without turbulence (Garagash, Detournay & Adachi Reference Garagash, Detournay and Adachi2011; Dontsov & Peirce Reference Dontsov and Peirce2015), while the addition of turbulence will enable one to better understand the transition from a laminar to a turbulent solution.

As follows from Garagash et al. (Reference Garagash, Detournay and Adachi2011), the problem of the tip region is equivalent to analysing a semi-infinite hydraulic fracture that propagates steadily in plane strain elastic conditions. With reference to figure 2, the volume balance for a steady semi-infinite hydraulic fracture that propagates with the velocity $V$ can be written as (see, e.g., Peirce & Detournay Reference Peirce and Detournay2008; Garagash et al. Reference Garagash, Detournay and Adachi2011)

(4.1) $$\begin{eqnarray}\displaystyle U=V+2C^{\prime }V^{1/2}\frac{x^{1/2}}{w}, & & \displaystyle\end{eqnarray}$$

where $U$ is the mean fluid velocity inside the fracture, $C^{\prime }=2C_{L}$ is the scaled Carter leak-off coefficient, while $w$ is the fracture aperture. The fluid velocity and fluid pressure are related through the flow (3.5) as (the flow is in the negative $x$ direction)

(4.2) $$\begin{eqnarray}\displaystyle U=\frac{w^{2}}{{\it\mu}^{\prime }\tilde{f}_{D}(Re,r)}\frac{\text{d}p}{\text{d}x}, & & \displaystyle\end{eqnarray}$$

where ${\it\mu}^{\prime }=12{\it\mu}$ is the scaled fluid viscosity and $\tilde{f}_{D}(Re,r)$ is given in (3.3). An elasticity equation that automatically satisfies the linear elastic fracture mechanics solution at the fracture tip can be written as (Dontsov & Peirce Reference Dontsov and Peirce2015)

(4.3a,b ) $$\begin{eqnarray}\displaystyle w=\frac{K^{\prime }}{E^{\prime }}x^{1/2}-\frac{4}{{\rm\pi}E^{\prime }}\int _{0}^{\infty }F(x,s)\frac{\text{d}p}{\text{d}s}\,\text{d}s,\quad F(x,s)=(s-x)\ln \left|\frac{x^{1/2}+s^{1/2}}{x^{1/2}-s^{1/2}}\right|-2x^{1/2}s^{1/2}, & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

where $K^{\prime }=\sqrt{32/{\rm\pi}}K_{Ic}$ , $K_{Ic}$ denotes the fracture toughness and $E^{\prime }=E/(1-{\it\nu}^{2})$ is the plane strain Young’s modulus.

Figure 2.

Schematics of a semi-infinite plane strain brittle hydraulic fracture with leak-off and turbulent flow.

By following Dontsov & Peirce (Reference Dontsov and Peirce2015) and introducing the scaling

(4.4a-d ) $$\begin{eqnarray}\displaystyle \tilde{w}=\frac{E^{\prime }w}{K^{\prime }x^{1/2}},\quad {\it\chi}=\frac{2C^{\prime }E^{\prime }}{V^{1/2}K^{\prime }},\quad \tilde{x}=(x/l)^{1/2},\quad \tilde{s}=(s/l)^{1/2},\quad l=\left(\frac{K^{\prime 3}}{{\it\mu}^{\prime }E^{\prime 2}V}\right)^{2}, & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

where $l$ is the length scale and the dimensionless parameter ${\it\chi}$ characterizes the leak-off, the governing equations (4.1)–(4.3) can be reduced to

(4.5a,b ) $$\begin{eqnarray}\displaystyle \tilde{w}(\tilde{x})=1+\frac{8}{{\rm\pi}}\int _{0}^{\infty }G\left(\frac{\tilde{s}}{\tilde{x}}\right)\tilde{f}_{D}\left[\frac{1}{\tilde{w}(\tilde{s})^{2}}+\frac{{\it\chi}}{\tilde{w}(\tilde{s})^{3}}\right]\text{d}\tilde{s},\quad G(t)=\frac{1-t^{2}}{t}\ln \left|\frac{1+t}{1-t}\right|+2. & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

The primary difference between the formulation in Dontsov & Peirce (Reference Dontsov and Peirce2015) and (4.5) is the presence of the multiplier $\tilde{f}_{D}(Re,r)$ , which captures the effect of turbulent fluid flow. It should be noted that the Reynolds number and roughness parameter (3.4) can be written using (4.4) as

(4.6a-d ) $$\begin{eqnarray}\displaystyle Re=\tilde{Re}(\tilde{w}+{\it\chi})\tilde{s},\quad r=\frac{\tilde{r}}{\tilde{s}\tilde{w}},\quad \tilde{Re}=\frac{24{\it\rho}K^{\prime 4}}{{\it\mu}^{\prime 2}E^{\prime 3}},\quad \tilde{r}=\frac{{\it\epsilon}{\it\mu}^{\prime }E^{\prime 3}V}{2K^{\prime 4}}. & & \displaystyle\end{eqnarray}$$

The phenomenological model for turbulent flow (3.3) applies for all Reynolds numbers and $r\lesssim 0.05$ . In this situation, large variations of the Reynolds number in (4.6) are admissible, while large variations of the roughness parameter are not permitted. At the same time, the roughness does not noticeably influence the laminar flow in the adopted model (3.3) (see figure 1), and therefore large values of $r$ are permitted as soon as the flow is not turbulent. To provide a quick estimate for the case of no leak-off (in which case $U=V$ ), one can calculate the value of $r$ from (3.4) as

(4.7) $$\begin{eqnarray}\displaystyle r=\frac{{\it\epsilon}{\it\rho}V}{{\it\mu}Re}. & & \displaystyle\end{eqnarray}$$

By substituting $Re=2000$ (just before the transition to turbulent flow), $V=1~\text{m}~\text{s}^{-1}$ , ${\it\epsilon}=10^{-4}$  m, and the properties of water ${\it\rho}=10^{3}~\text{kg}~\text{m}^{-3}$ , ${\it\mu}=10^{-3}$  Pa s into (4.7), one obtains $r=0.05$ . This shows that care must be taken when applying this approach to make sure that the model is not used beyond its applicability region.

4.2 Estimation of the laminar zone size

In order to estimate the influence of turbulence, it is instructive to know at what distance away from the fracture tip the turbulent effects become important. To this end, (4.6) can be rewritten as

(4.8) $$\begin{eqnarray}\displaystyle \frac{Re}{\tilde{Re}}=(\tilde{w}(\tilde{x},{\it\chi})+{\it\chi})\tilde{x}. & & \displaystyle\end{eqnarray}$$

Here, $\tilde{w}(\tilde{x},{\it\chi})$ should be a solution of (4.5) for laminar flow, i.e. $\tilde{f}_{D}=1$ , since the turbulence becomes important only some distance away from the tip. To produce a rapid estimate, one can use a closed-form approximate solution of (4.5) for $\tilde{f}_{D}=1$ , which has been obtained in Dontsov & Peirce (Reference Dontsov and Peirce2015). Figure 3 shows the variation of $Re/\tilde{Re}$ versus $\tilde{x}$ that is calculated based on (4.8) and using the approximate width solution (referred to as the universal solution for laminar flow) from Dontsov & Peirce (Reference Dontsov and Peirce2015). The figure also shows the variation of the aforementioned approximate solution $\tilde{w}$ versus the normalized spatial coordinate $\tilde{x}$ . The values of the leak-off parameter used for the plot are ${\it\chi}=\{0,1,10,10^{2},10^{3}\}$ . The curves that correspond to larger values of ${\it\chi}$ are above those that represent smaller values of ${\it\chi}$ , which applies to both panels.

Figure 3.

Variation of $Re/\tilde{Re}$  (4.8) (a) and $\tilde{w}$ (b) versus $\tilde{x}$ for ${\it\chi}=\{0,1,10,10^{2},10^{3}\}$ (solid black lines). The fracture width solution shown in (b) is calculated using an approximation obtained by Dontsov & Peirce (Reference Dontsov and Peirce2015). The dashed coloured lines indicate the vertex solutions (4.9).

The dashed red, blue and green lines correspond to the vertex solutions (asymptotic solutions of (4.5) for laminar flow $\tilde{f}_{D}=1$ , see, e.g., Garagash et al. (Reference Garagash, Detournay and Adachi2011)), which are given in the scaled form (4.4) by

(4.9a-c ) $$\begin{eqnarray}\displaystyle \tilde{w}_{k}=1,\quad \tilde{w}_{\tilde{m}}={\it\beta}_{\tilde{m}}{\it\chi}^{1/4}\tilde{x}^{1/4},\quad \tilde{w}_{m}={\it\beta}_{m}\tilde{x}^{1/3}, & & \displaystyle\end{eqnarray}$$

where the constants are ${\it\beta}_{\tilde{m}}=4/(15^{1/4}(\sqrt{2}-1)^{1/4})$ and ${\it\beta}_{m}=2^{1/3}3^{5/6}$ , $\tilde{w}_{k}$ stands for the toughness solution, $\tilde{w}_{m}$ represents the viscous solution and $\tilde{w}_{\tilde{m}}$ corresponds to the leak-off solution. As can be seen from figure 3, the estimation of the normalized Reynolds number from the leak-off solution (green dashed line in a) coincides with the results of the corresponding toughness and viscous solutions. This occurs because $1/\tilde{w}^{2}\ll {\it\chi}/\tilde{w}^{3}$ (or ${\it\chi}\gg \tilde{w}$ ) in the leak-off regime (see (4.5) and discussion in Dontsov & Peirce (Reference Dontsov and Peirce2015)), which leads to $Re/\tilde{Re}=(\tilde{w}+{\it\chi})\tilde{x}\approx {\it\chi}\tilde{x}$ . This shows that the estimate of the normalized Reynolds number is independent of the fracture opening solution for the leak-off regime. As a result, one can use the ‘ $mk$ ’ solution (which is the no-leak-off laminar solution of (4.5), i.e. ${\it\chi}=0$ and $\tilde{f}_{D}=1$ ) in (4.8) for the purpose of obtaining a crude estimate, i.e. $Re/\tilde{Re}=(\tilde{w}_{mk}+{\it\chi})\tilde{x}$ . The latter estimate for $Re/\tilde{Re}$ captures the small-leak-off case accurately, while the large-leak-off limit is captured automatically as discussed above. It should be noted that such an approximation (replacing $\tilde{w}$ by $\tilde{w}_{mk}$ in (4.8)) causes only a 10 % error in $Re/\tilde{Re}$ values for the range of parameters shown in figure 3. To simplify the calculations, the ‘ $mk$ ’ solution can be approximated as $\tilde{w}_{mk}\approx (1+{\it\beta}_{m}^{3}\tilde{x})^{1/3}$ , which introduces an additional error that is of the order of 1 %, see Dontsov & Peirce (Reference Dontsov and Peirce2015).

The procedure for estimating the size of the laminar zone can be formulated as follows. First, the fraction $Re/\tilde{Re}$ and leak-off parameter ${\it\chi}$ should be calculated based on the values of the problem parameters, and taking $Re=2000$ (turbulent effects become important for $Re\gtrsim 2000$ ). Then, figure 3 can be used to calculate $\tilde{x}$ , which can be converted to a physical length using (4.4) as $d=l\tilde{x}^{2}$ . For instance, by taking $K_{Ic}=0.16~\text{MPa}~\text{m}^{1/2}$ , $E^{\prime }=6.8$  GPa, ${\it\rho}=1000~\text{kg}~\text{m}^{-3}$ , $C^{\prime }=0$ , ${\it\mu}=1.8\times 10^{-3}~\text{Pa}~\text{s}$ and $V=2~\text{m}~\text{s}^{-1}$ , which correspond to the parameters used in Tsai & Rice (Reference Tsai and Rice2010) for modelling glacial cracks, one obtains $d\approx 2$  m ( $d\approx 5$ m if one accounts for the bi-material effect as in Tsai & Rice (Reference Tsai and Rice2010)). In other words, the size of the laminar zone is 2 m; however, this does not imply that the solution is turbulent immediately further away from the tip since there is always a transition zone. It should be noted that the size of the laminar zone strongly depends on the velocity, and for $V=0.5~\text{m}~\text{s}^{-1}$ (other parameters are kept the same) becomes $d\approx 30$  m. To address petroleum applications, consider a slick-water treatment with $K_{Ic}=1~\text{MPa}~\text{m}^{1/2}$ , $E^{\prime }=20$  GPa, ${\it\rho}=1000~\text{kg}~\text{m}^{-3}$ , $C^{\prime }=0$ , ${\it\mu}=10^{-3}$  Pa s and $V=1~\text{m}~\text{s}^{-1}$ . In this situation, $d\approx 7$  m. As expected, the size of the laminar zone also strongly depends on the viscosity, and for ${\it\mu}=10^{-1}$  Pa s one obtains $d\approx 730$  m. Therefore, no turbulence is expected near the tip. In order to provide similar estimates for different parameters, one may use figure 3 (and the procedure described in this section) to estimate the size of the laminar zone.

For simple fracture geometries, such as a plane strain or a radial fracture, the value of $d$ can be compared with the fracture length to quantify the global importance of the laminar tip. At the same time, for numerical schemes the value of $d$ needs to be compared with the size of the tip element, which would allow one to select a suitable asymptotic solution for calculations. For instance, if the element size is $1$  m, while $d=7$  m, then the laminar asymptotic solution should be used for the tip element to advance the fracture front even if turbulent flow develops further inside the fracture.

4.3 ‘Turbulent’ asymptotic solutions

This section focuses on calculating approximate solutions in situations where the effect of turbulent flow dominates the response. Two possibilities are considered: (i) a viscosity-dominated turbulent solution, labelled $\tilde{w}_{t}$ , and (ii) a leak-off-dominated turbulent solution, labelled $\tilde{w}_{\tilde{t}}$ . The corresponding governing equations can be obtained from (4.5) by retaining only the viscous or leak-off terms, and can be written respectively as

(4.10a,b ) $$\begin{eqnarray}\displaystyle \tilde{w}_{t}(\tilde{x})=\frac{8}{{\rm\pi}}\int _{0}^{\infty }G\left(\frac{\tilde{s}}{\tilde{x}}\right)\frac{\tilde{f}_{D}}{\tilde{w}_{t}^{2}}\,\text{d}\tilde{s},\quad \tilde{w}_{\tilde{t}}(\tilde{x})=\frac{8}{{\rm\pi}}\int _{0}^{\infty }G\left(\frac{\tilde{s}}{\tilde{x}}\right)\frac{{\it\chi}\tilde{f}_{D}}{\tilde{w}_{\tilde{t}}^{3}}\,\text{d}\tilde{s}. & & \displaystyle\end{eqnarray}$$

It should be noted that there is no toughness-dominated turbulent solution since for the latter the flow does not affect the solution. The friction factor (3.3), reduced for the turbulent regime ( $Re\gtrsim 3000$ ), and the friction factor (3.2) can be written as

(4.11a,b ) $$\begin{eqnarray}\displaystyle \tilde{f}_{D}\approx \frac{Re}{8.511^{2}}(f_{0}(Re,r))^{-2},\quad \tilde{f}_{D,2}=1.88\times 10^{-3}r^{1/3}Re. & & \displaystyle\end{eqnarray}$$

Further, the Reynolds number (4.6) for the viscous turbulent solution can be reduced to $Re=\tilde{Re}\tilde{w}\tilde{s}$ , while its counterpart for the leak-off turbulent solution is $Re={\it\chi}\tilde{Re}\tilde{s}$ , and the roughness parameter is the same for both solutions, $r=\tilde{r}/(\tilde{s}\tilde{w})$ .

By seeking the solution in the form $\tilde{w}_{t}={\it\beta}_{t}\tilde{x}^{{\it\gamma}}$ , the substitution of (4.11) into (4.10) leads to the following expressions for the viscosity-dominated and leak-off-dominated turbulent solutions for both friction factors in (4.11):

(4.12) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \tilde{w}_{t}={\it\beta}_{t}\tilde{x},\quad {\it\beta}_{t}=0.4165\tilde{Re}^{1/2}\left[f_{0}\left(\tilde{Re}{\it\beta}_{t,0}\tilde{x}^{2},\frac{\tilde{r}}{{\it\beta}_{t,0}\tilde{x}^{2}}\right)\right]^{-1},\quad {\it\beta}_{t,0}\approx 0.08\tilde{Re}^{1/2},\\ \displaystyle \tilde{w}_{t,2}={\it\beta}_{t,2}\tilde{x}^{5/7},\quad {\it\beta}_{t,2}=0.221\tilde{Re}^{3/7}\tilde{r}^{1/7},\\ \displaystyle \!\!\!\tilde{w}_{\tilde{t}}={\it\beta}_{\tilde{t}}\tilde{x}^{1/2},\quad {\it\beta}_{\tilde{t}}=0.737\tilde{Re}^{1/4}{\it\chi}^{1/2}\left[f_{0}\!\left(\!{\it\chi}\tilde{Re}\tilde{x},\frac{\tilde{r}}{{\it\beta}_{\tilde{t},0}\tilde{x}^{5/3}}\right)\right]^{-1/2},\quad {\it\beta}_{\tilde{t},0}\approx 0.33\tilde{Re}^{1/4}{\it\chi}^{1/4},\\ \displaystyle \tilde{w}_{\tilde{t},2}={\it\beta}_{\tilde{t},2}\tilde{x}^{5/13},\quad {\it\beta}_{\tilde{t},2}=0.504\tilde{Re}^{3/13}{\it\chi}^{6/13}\tilde{r}^{1/13},\end{array}\right\} & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

where the additional subscript ‘ $2$ ’ refers to the case where the second approximation for the friction factor in (4.11) is used (solutions with no additional subscript utilize the first approximation), the viscosity-dominated turbulent solutions are labelled with the subscript ‘ $t$ ’ and the leak-off-dominated turbulent solutions are augmented with the subscript ‘ $\tilde{t}$ ’. Here, the logarithmic term is treated as a slowly varying function, and the numeric factors in ${\it\beta}_{t,0}$ and ${\it\beta}_{\tilde{t},0}$ are calculated from estimating $|0.9\log (7/Re)|$ for $Re=2000$ .

It should be noted that the unscaled solutions, which can be calculated from (4.12) using the scaling (4.4), vary with distance as $w_{t}\propto x$ , $w_{t,2}\propto x^{6/7}$ , $w_{\tilde{t}}\propto x^{3/4}$ and $w_{\tilde{t},2}\propto x^{9/13}$ . It should also be noted that the solution $\tilde{w}_{t,2}$  (4.12) is consistent with the one used in Tsai & Rice (Reference Tsai and Rice2010), where the latter was calculated based on the solution for a power-law fluid (Desroches et al. Reference Desroches, Detournay, Lenoach, Papanastasiou, Pearson, Thiercelin and Cheng1994). In order to obtain identical agreement one needs to remove the bi-material effect from the solution by Tsai & Rice (Reference Tsai and Rice2010), and to account for the numerical factors appearing in the expressions for $\tilde{Re}$ and $\tilde{r}$ in (4.6).

4.4 Numerical solution

To illustrate the applicability of the laminar zone size estimate covered in § 4.2 and the ‘turbulent’ asymptotic solutions described in § 4.3, figure 4 shows the numerical solutions of (4.5) for various parameters. In particular, the solutions are calculated for the normalized Reynolds numbers $\tilde{Re}=\{0,10,10^{3},10^{5},10^{7}\}$ , for leak-off parameter ${\it\chi}=0$ and ${\it\chi}=10$ , and for normalized roughness $\tilde{r}=0$ and $\tilde{r}=60/\tilde{R}e$ . The latter choice for the roughness is based on the estimate (4.7), for which $r=0.03$ is chosen, and ensures that the adopted model is not used beyond its applicability region. It should be noted that the aforementioned estimate is valid for the no-leak-off case and becomes an approximation for the case ${\it\chi}=10$ . The estimate of the laminar zone size that is based on (4.8), the vertex asymptotic solutions (4.9) and the ‘turbulent’ asymptotic solutions (4.12) are all shown as well.

Figure 4.

Numerical solutions of (4.5) calculated using (3.3) for $\tilde{Re}=\{0,10,10^{3},10^{5},10^{7}\}$ (solid black lines). The lowest curve in each panel corresponds to the purely laminar solution ( $\tilde{Re}=0$ ), while the highest curves represent the $\tilde{Re}=10^{7}$ cases. The (a,c) solutions without leak-off, i.e. ${\it\chi}=0$ ; (b,d) solutions for ${\it\chi}=10$ . (a,b) Solutions in smooth fractures, for which the roughness is zero, i.e. $\tilde{r}=0$ ; (c,d) solutions for which the roughness is calculated as $\tilde{r}=60/\tilde{R}e$ . The estimate of the laminar zone size, based on (4.8), is highlighted by the black filled circular markers. The vertex asymptotic solutions (4.9) are indicated by dashed red, blue and green lines, while the ‘turbulent’ asymptotic solutions (4.12) are indicated by dashed and solid magenta and grey lines. The square markers indicate the locations at which the Reynolds numbers associated with the flow in the fracture are $Re=\{10^{4},10^{5},10^{6}\}$ .

As can be seen from figure 4, larger values of $\tilde{Re}$ correspond to situations in which turbulent effects occur closer to the tip. The solution may ‘depart’ from either one of the vertex solutions (4.9) (dashed red, blue and green lines), or from a transition zone between the aforementioned solutions. The estimate of the laminar zone size (4.8) is able to estimate the validity of the laminar solution accurately for almost all cases. One exception is when the turbulent flow occurs while hydraulic fracture is propagating in the toughness regime. For the latter case, the estimate (4.8) is on the safe side, as it underestimates the validity region of the laminar solution. It should be noted that this occurs because the fluid flow does not affect the solution in the toughness regime, so that one can use the toughness solution even if the flow is turbulent. The ‘turbulent’ asymptotic solutions (4.12) calculated using (3.3), $\tilde{w}_{t}$ , and $\tilde{w}_{\tilde{t}}$ , provide a reasonably accurate estimate as soon as they apply. It should be noted that the accuracy is smaller than for the vertex solutions, which is due to additional approximations associated with the logarithmic term. On the other hand, the asymptotic solutions based on (3.2), $\tilde{w}_{t,2}$ and $\tilde{w}_{\tilde{t},2}$ , do not apply for smooth fractures, while they provide a good estimate for rough fractures. This is consistent with the results shown in figure 1, where both models give similar values for the friction factor in the range $10^{-3}<{\it\epsilon}/(2w)<5\times 10^{-2}$ . It should be noted that the roughness parameter becomes smaller for wider fractures, in which case the solutions $\tilde{w}_{t,2}$ and $\tilde{w}_{\tilde{t},2}$ inevitably divert from the numerical solutions. It should also be noted that all leak-off asymptotes require ${\it\chi}\gg 1$ , while relatively modest ${\it\chi}=10$ is used for plotting figure 4(b,d). For this reason, there is an additional discrepancy between the numerical solution and the asymptotes.

Since elasticity produces a non-local effect and there is a transition zone between laminar and fully developed turbulent flows, the turbulent asymptotic solutions do no apply immediately after the flow becomes non-laminar. There is a transition zone, which, according to figure 4, occurs approximately until $Re\gtrsim 10^{4}$ (the values are larger for solutions ‘departing’ from the toughness asymptote). Due to the presence of the transition zone, the estimate of the laminar zone size, described in § 4.2, cannot be used directly to quantify the distance from the fracture tip after which the turbulent asymptotic solution can be applied. In order to have a better estimate, one should multiply the calculated laminar zone size $\tilde{x}$ by a factor of 3–5 (this corresponds to the distance between the circular marker and the first square marker), in which case the final answer becomes at least an order of magnitude larger due to the scaling (4.4) since $x\propto \tilde{x}^{2}$ . This again shows that care must be taken when using turbulent flow models and the associated asymptotic solutions for hydraulic fracture modelling to make sure that the solutions are not used beyond their applicability region.