2 Model description for flow in a narrow fracture

The HB model for a shear-thinning/thickening fluid with yield stress is

(2.1) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D70F}=(\unicode[STIX]{x1D707}_{0}\dot{\unicode[STIX]{x1D6FE}}^{n-1}+\unicode[STIX]{x1D70F}_{p}\dot{\unicode[STIX]{x1D6FE}}^{-1})\dot{\unicode[STIX]{x1D6FE}},\quad \unicode[STIX]{x1D70F}\geqslant \unicode[STIX]{x1D70F}_{p},\\ \dot{\unicode[STIX]{x1D6FE}}=0,\quad \unicode[STIX]{x1D70F}<\unicode[STIX]{x1D70F}_{p},\end{array}\right\}\end{eqnarray}$$

given in terms of the stress $\unicode[STIX]{x1D70F}$ and the strain rate $\dot{\unicode[STIX]{x1D6FE}}$ . A slightly more complicated description using tensor invariants is required for three-dimensional flows; this formulation is not reported here as we consider a one-dimensional problem below. The parameter $\unicode[STIX]{x1D707}_{0}$ , the consistency index, represents a viscosity-like parameter, while $\unicode[STIX]{x1D70F}_{p}$ is the yield stress of the fluid, and $n$ , the fluid behaviour index, controls the extent of shear thinning ( $n<1$ ) or shear thickening ( $n>1$ ); $n=1$ corresponds to the Bingham case. For flow through a narrow fracture (such as a Hele-Shaw cell) of width $L_{y}$ , as depicted in figure 1, the primary balance is between cross-gap quantities. The relevant relationship is that the velocity $u(x,y)$ in the $x$ -direction must satisfy

(2.2) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D70F}_{xy}=\left(\unicode[STIX]{x1D707}_{0}\left|{\displaystyle \frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}y}}\right|^{n-1}+\unicode[STIX]{x1D70F}_{p}\left|{\displaystyle \frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}y}}\right|^{-1}\right){\displaystyle \frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}y}},\quad \unicode[STIX]{x1D70F}_{xy}\geqslant \unicode[STIX]{x1D70F}_{p},\\ {\displaystyle \frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}y}}=0,\quad \unicode[STIX]{x1D70F}_{xy}<\unicode[STIX]{x1D70F}_{p},\end{array}\right\}\end{eqnarray}$$

where $\unicode[STIX]{x1D70F}_{xy}$ represents the cross-gap stress and $y$ is the cross-gap direction. Furthermore, the primary balance between the cross-gap stress and the along-fracture pressure gradient is given by $y(\unicode[STIX]{x2202}p/\unicode[STIX]{x2202}x)=\unicode[STIX]{x1D70F}_{xy}$ . We define the height $h(x,t)$ of fluid above a horizontal impermeable base. Then, under the relevant shallow water approximation, the pressure is to leading order hydrostatic, and the pressure gradient becomes $\unicode[STIX]{x2202}p/\unicode[STIX]{x2202}x=\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}\,g\,(\unicode[STIX]{x2202}h/\unicode[STIX]{x2202}x)$ , where $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$ is the density difference driving the flow, between fluid within and that outside the gravity current, and $g$ is the acceleration due to gravity. We assume that the wetting characteristics at the walls of the fracture are unimportant (or alternatively that the flow takes place in a prewetted fracture). Combination of these equations, using the condition that $\unicode[STIX]{x1D70F}_{xy}=\unicode[STIX]{x1D70F}_{p}$ as the fluid yields, leads to an expression for the location of the yield surface in the gap, $|y_{yield}|=\unicode[STIX]{x1D70F}_{p}/(\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}\,g|\unicode[STIX]{x2202}h/\unicode[STIX]{x2202}x|)$ , where $-L_{y}/2\leqslant y_{yield}\leqslant L_{y}/2$ is required.

To continue, we need to solve for the cross-gap flow structure in yielded and unyielded regions of the flow. The continuity of mass of the fluid layer may be written as

(2.3) $$\begin{eqnarray}{\displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}t}}=-{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}}(\bar{u}h),\end{eqnarray}$$

where $\bar{u}(x,t)$ is the gap-averaged velocity. Assuming a zero slip velocity, and upon introducing the expression for $\bar{u}(x,t)$ into (2.3), we can form an evolution equation for $h(x,t)$ alone, namely

(2.4) $$\begin{eqnarray}\displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}t}} & = & \displaystyle \left(\frac{L_{y}}{2}\right)^{(n+1)/n}\text{sgn}\left({\displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}x}}\right)\left(\frac{n}{2n+1}\right)\left(\frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}\,g}{\unicode[STIX]{x1D707}_{0}}\right)^{1/n}\nonumber\\ \displaystyle & & \displaystyle \times \,{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}}\left[h\left|{\displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}x}}\right|^{1/n}\left(1-\unicode[STIX]{x1D705}\left|{\displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}x}}\right|^{-1}\right)^{(n+1)/n}\left(1+\left(\frac{n}{n+1}\right)\unicode[STIX]{x1D705}\left|{\displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}x}}\right|^{-1}\right)\right],\end{eqnarray}$$

where $\unicode[STIX]{x1D705}=2\unicode[STIX]{x1D70F}_{p}/(\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}\,gL_{y})$ is a non-dimensional number representing the ratio between yield stress and gravity related stress, or the ratio between the Bingham and Ramberg numbers. Additionally, setting

(2.5) $$\begin{eqnarray}\unicode[STIX]{x1D6FA}=\left(\frac{L_{y}}{2}\right)^{(n+1)/n}\left(\frac{n}{2n+1}\right)\left(\frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}\,g}{\unicode[STIX]{x1D707}_{0}}\right)^{1/n},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FA}$ is a velocity scale, allows us to write this equation slightly more succinctly as

(2.6) $$\begin{eqnarray}{\displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}t}}=\text{sgn}\left({\displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}x}}\right)\unicode[STIX]{x1D6FA}{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}}\left[h\left|{\displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}x}}\right|^{1/n}\left(1-\unicode[STIX]{x1D705}\left|{\displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}x}}\right|^{-1}\right)^{(n+1)/n}\left(1+\left(\frac{n}{n+1}\right)\unicode[STIX]{x1D705}\left|{\displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}x}}\right|^{-1}\right)\right].\end{eqnarray}$$

To this equation we must add the further condition that when the flow is fully plugged (that is $y_{yield}=L_{y}/2$ ), then the velocity throughout the gap is zero ( $\bar{u}=0$ ), so therefore $\unicode[STIX]{x2202}h/\unicode[STIX]{x2202}t=0$ for such regions. This is the limiting version of the equation above in the limit $\unicode[STIX]{x2202}h/\unicode[STIX]{x2202}x=\unicode[STIX]{x1D705}$ , which is in turn the condition for the flow to be fully plugged. Consequently, the equation for the height of the current is continuous through such a plugging transition.

In the model, the slip contribution was neglected also because most experiments were conducted by roughening the Hele-Shaw cell with commercial transparent antislip tape, which is commonly adopted to make slippery surfaces safe. Independent rheometric measurements were conducted with plates roughened with sandpaper.

3 Self-similar solution

Various self-similar solutions exist to describe the spreading of gravity currents of constant or variable volume in cases similar to our problem where $\unicode[STIX]{x1D705}=0,n=1$ (e.g. King & Woods Reference King and Woods2003; Lyle et al. Reference Lyle, Huppert, Hallworth, Bickle and Chadwick2005) and where $\unicode[STIX]{x1D705}=0$ (Pascal & Pascal Reference Pascal and Pascal1993; Longo et al. Reference Longo, Di Federico, Chiapponi and Archetti2013b ; Longo, Di Federico & Chiapponi Reference Longo, Di Federico and Chiapponi2015b ; Ciriello et al. Reference Ciriello, Longo, Chiapponi and Di Federico2016). For the current study where both yield stress and shear-thinning/thickening effects are present, a generalized form of such solutions is not available since the presence of yield stress breaks self-similarity. However, one form of self-similar solution does exist if we allow a variable injection rate of fluid (a more general type of solution with self-similarity of the second kind might be possible, although it is not attempted here). To this end, suppose that the volume of fluid in the gravity current varies as

(3.1) $$\begin{eqnarray}L_{y}\int _{0}^{\infty }h(x,t)\,\text{d}x=Qt^{\unicode[STIX]{x1D6FC}},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FC},Q\geqslant 0$ . Then, the principal dimensions of the three parameters are $[\unicode[STIX]{x1D6FA}]=L/T,[Q/L_{y}]=L^{2}/T^{\unicode[STIX]{x1D6FC}},[\unicode[STIX]{x1D705}]=1.$ It is possible to rewrite our principal variables in dimensionless form as

(3.2a-c ) $$\begin{eqnarray}\widetilde{h}=h\left(\frac{Q}{L_{y}\unicode[STIX]{x1D6FA}^{\unicode[STIX]{x1D6FC}}}\right)^{1/(\unicode[STIX]{x1D6FC}-2)},\quad \widetilde{x}=x\left(\frac{Q}{L_{y}\unicode[STIX]{x1D6FA}^{\unicode[STIX]{x1D6FC}}}\right)^{1/(\unicode[STIX]{x1D6FC}-2)},\quad \widetilde{t}=t\left(\frac{Q}{L_{y}\unicode[STIX]{x1D6FA}^{2}}\right)^{1/(\unicode[STIX]{x1D6FC}-2)}.\end{eqnarray}$$

Introduction of these variables immediately reduces (2.6)–(3.1) to the dimensionless form

(3.3) $$\begin{eqnarray}{\displaystyle \frac{\unicode[STIX]{x2202}\widetilde{h}}{\unicode[STIX]{x2202}\widetilde{t}}}=-{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\widetilde{x}}}\left[\widetilde{h}\left|{\displaystyle \frac{\unicode[STIX]{x2202}\widetilde{h}}{\unicode[STIX]{x2202}\widetilde{x}}}\right|^{1/n}\left(1-\unicode[STIX]{x1D705}\left|{\displaystyle \frac{\unicode[STIX]{x2202}\widetilde{h}}{\unicode[STIX]{x2202}\widetilde{x}}}\right|^{-1}\right)^{(n+1)/n}\left(1+\left(\frac{n}{n+1}\right)\unicode[STIX]{x1D705}\left|{\displaystyle \frac{\unicode[STIX]{x2202}\widetilde{h}}{\unicode[STIX]{x2202}\widetilde{x}}}\right|^{-1}\right)\right],\end{eqnarray}$$

(3.4) $$\begin{eqnarray}\int _{0}^{\infty }\widetilde{h}\,\text{d}\widetilde{x}=\widetilde{t}^{\,\unicode[STIX]{x1D6FC}},\end{eqnarray}$$

where we have assumed that $\unicode[STIX]{x2202}\widetilde{h}/\unicode[STIX]{x2202}\widetilde{x}<0$ . Hereafter, the tilde is dropped. We can seek self-similar solutions of these equations by looking for solutions in the form $h=t^{\unicode[STIX]{x1D6FD}}f(\unicode[STIX]{x1D702})$ , where $\unicode[STIX]{x1D702}=x/t^{\unicode[STIX]{x1D6FE}}$ . Substitution into (3.3)–(3.4) yields

(3.5) $$\begin{eqnarray}\displaystyle & \hspace{-6.0pt}\displaystyle \unicode[STIX]{x1D6FD}t^{\unicode[STIX]{x1D6FD}-1}f-t^{\unicode[STIX]{x1D6FD}-1}\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D702}f^{\prime }=-\!\left[t^{\unicode[STIX]{x1D6FD}+(\unicode[STIX]{x1D6FD}-\unicode[STIX]{x1D6FE})/n}f|\,f^{\prime }|^{1/n}\!\left(\!1-\frac{\unicode[STIX]{x1D705}t^{\unicode[STIX]{x1D6FE}-\unicode[STIX]{x1D6FD}}}{|f^{\prime }|}\right)^{(n+1)/n}\!\!\left(\!1+\!\left(\frac{n}{n+1}\right)\!\frac{\unicode[STIX]{x1D705}t^{\unicode[STIX]{x1D6FE}-\unicode[STIX]{x1D6FD}}}{|f^{\prime }|}\right)\right]^{\prime }\!\!t^{-\unicode[STIX]{x1D6FE}}, & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

(3.6) $$\begin{eqnarray}\displaystyle & \displaystyle \int _{0}^{\unicode[STIX]{x1D702}_{e}}t^{\unicode[STIX]{x1D6FD}+\unicode[STIX]{x1D6FE}}f(\unicode[STIX]{x1D702})\,\text{d}\unicode[STIX]{x1D702}=t^{\unicode[STIX]{x1D6FC}}, & \displaystyle\end{eqnarray}$$

where primes denote differentiation with respect to $\unicode[STIX]{x1D702}$ . Comparison of powers of $t$ gives immediately $\unicode[STIX]{x1D6FD}+\unicode[STIX]{x1D6FE}=\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FE}-\unicode[STIX]{x1D6FD}=0,\unicode[STIX]{x1D6FD}-1=(\unicode[STIX]{x1D6FD}-\unicode[STIX]{x1D6FE})(1+1/n)$ . This has one consistent solution, namely $\unicode[STIX]{x1D6FC}=2$ and $\unicode[STIX]{x1D6FD}=\unicode[STIX]{x1D6FE}=1$ . A fluid injection rate scaling as $t$ allows a self-similar solution, for which both the height and the length of the current increase with time, similarly scaling as $t$ . However, for $\unicode[STIX]{x1D6FC}=2$ , the scales in (3.2) break down and an additional velocity scale embedded in the integral constraint of mass conservation given by (3.1) arises beyond $\unicode[STIX]{x1D6FA}$ , given by $(Q/L_{y})^{1/2}$ . A similar case is treated in Di Federico, Archetti & Longo (Reference Di Federico, Archetti and Longo2012a ,Reference Di Federico, Archetti and Longo b ). We define an arbitrary time scale $t^{\ast }$ , the velocity scale $u^{\ast }=(Q/L_{y})^{1/2}$ and the ratio between the two velocity scales $\unicode[STIX]{x1D6FF}=(L_{y}\unicode[STIX]{x1D6FA}^{2}/Q)^{1/2}$ , with $x^{\ast }=u^{\ast }t^{\ast }$ . Equations (3.3)–(3.4) become

(3.7) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}t}}=-\unicode[STIX]{x1D6FF}{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}}\left[h\left|{\displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}x}}\right|^{1/n}\left(1-\unicode[STIX]{x1D705}\left|{\displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}x}}\right|^{-1}\right)^{(n+1)/n}\left(1+\left(\frac{n}{n+1}\right)\unicode[STIX]{x1D705}\left|{\displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}x}}\right|^{-1}\right)\right], & \displaystyle\end{eqnarray}$$

(3.8) $$\begin{eqnarray}\displaystyle & \displaystyle \int _{0}^{\infty }h(x,t)\,\text{d}x=t^{2}. & \displaystyle\end{eqnarray}$$

We seek a self-similar solution of the form $h=tf(\unicode[STIX]{x1D702}),\text{where}\quad \unicode[STIX]{x1D702}=x/t$ ; substitution into (3.7)–(3.8) yields

(3.9) $$\begin{eqnarray}\displaystyle & \displaystyle f-\unicode[STIX]{x1D702}f^{\prime }=-\unicode[STIX]{x1D6FF}\left[f|\,f^{\prime }|^{1/n}\left(1-\frac{\unicode[STIX]{x1D705}}{|f^{\prime }|}\right)^{(n+1)/n}\left(1+\left(\frac{n}{n+1}\right)\frac{\unicode[STIX]{x1D705}}{|f^{\prime }|}\right)\right]^{\prime }, & \displaystyle\end{eqnarray}$$

(3.10) $$\begin{eqnarray}\displaystyle & \displaystyle \int _{0}^{\unicode[STIX]{x1D702}_{e}}f(\unicode[STIX]{x1D702})\,\text{d}\unicode[STIX]{x1D702}=1. & \displaystyle\end{eqnarray}$$

This system admits a simple solution, namely a linear profile for $f(\unicode[STIX]{x1D702})$ (Di Federico et al. Reference Di Federico, Archetti and Longo2012a ). Supposing a solution in the form $f(\unicode[STIX]{x1D702})=A(\unicode[STIX]{x1D702}_{e}-\unicode[STIX]{x1D702})$ , for some constants $A>0$ and $\unicode[STIX]{x1D702}_{e}>0$ , we substitute into (3.9)–(3.10) to obtain

(3.11) $$\begin{eqnarray}A\unicode[STIX]{x1D702}_{e}=\unicode[STIX]{x1D6FF}A^{(n+1)/n}\left(1-\frac{\unicode[STIX]{x1D705}}{A}\right)^{(n+1)/n}\left[1+\left(\frac{n}{n+1}\right)\frac{\unicode[STIX]{x1D705}}{A}\right],\quad \unicode[STIX]{x1D702}_{e}=\sqrt{\frac{2}{A}}.\end{eqnarray}$$

Elimination of $\unicode[STIX]{x1D702}_{e}$ gives one nonlinear equation to solve for $A$ .

Solutions for several values of the parameter $n$ are shown in figure 2, with all other parameters kept constant. The self-similar solution retains similarity by constraining the yield surface to be at a constant location along its length as the current spreads. Furthermore, the thickness of the plugged region simply grows as the fluid becomes more shear thinning.

Figure 2. (a) Shape of the similarity solution in a Hele-Shaw cell ( $\unicode[STIX]{x1D6FC}=2$ ) for different values of $n$ ; (b) plug regions.

3.1 Asymptotic analysis for $\unicode[STIX]{x1D6FC}\neq 2$

For $\unicode[STIX]{x1D6FC}\neq 2$ , no self-similar solution is predicted, but it is possible to find an approximate result starting from the self-similar solution for power-law fluids ( $\unicode[STIX]{x1D705}\rightarrow 0$ ). We briefly recall that for $\unicode[STIX]{x1D705}\rightarrow 0$ equation (3.3) becomes

(3.12) $$\begin{eqnarray}{\displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}x}}=-{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}}\left(h\left|{\displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}x}}\right|^{1/n}\right),\end{eqnarray}$$

while the integral mass balance given by (3.4) is unmodified. Equations (3.4)–(3.12) admit the similarity solution

(3.13a-d ) $$\begin{eqnarray}h=\unicode[STIX]{x1D702}_{N}^{n+2}t^{\,F_{2}}f(\unicode[STIX]{x1D701}),\quad \unicode[STIX]{x1D702}=xt^{\,-F_{1}},\quad \unicode[STIX]{x1D702}_{N}=\left(\int _{0}^{1}f\,\text{d}\unicode[STIX]{x1D701}\right)^{-1/(n+2)},\quad \unicode[STIX]{x1D701}=\unicode[STIX]{x1D702}/\unicode[STIX]{x1D702}_{N},\end{eqnarray}$$

where

(3.14a,b ) $$\begin{eqnarray}F_{1}={\displaystyle \frac{\unicode[STIX]{x1D6FC}+n}{2+n}},\quad F_{2}=\unicode[STIX]{x1D6FC}-F_{1},\end{eqnarray}$$

and where the shape function $f$ satisfies the following nonlinear ordinary differential equation:

(3.15) $$\begin{eqnarray}(f|f^{\prime }|^{1/n})^{\prime }+F_{2}\,f-F_{1}\unicode[STIX]{x1D701}f^{\prime }=0.\end{eqnarray}$$

The numerical integration of (3.15) for $\unicode[STIX]{x1D6FC}\neq 0$ and $\unicode[STIX]{x1D6FC}\neq 2$ requires two boundary conditions at $\unicode[STIX]{x1D701}\rightarrow 1$ . By assuming $f\approx a_{0}(1-\unicode[STIX]{x1D701})^{b}$ , substituting in (3.15) and balancing the lower-order terms, $b=1$ and $a_{0}=F_{2}^{n}$ are yielded. Hence, it follows that

(3.16a,b ) $$\begin{eqnarray}\left.f\right|_{\unicode[STIX]{x1D701}\rightarrow 1-\unicode[STIX]{x1D716}}=F_{2}^{n}\unicode[STIX]{x1D716},\quad \left.f^{\prime }\right|_{\unicode[STIX]{x1D701}\rightarrow 1-\unicode[STIX]{x1D716}}=-F_{2}^{n},\end{eqnarray}$$

with $\unicode[STIX]{x1D716}$ a small quantity. The two cases $\unicode[STIX]{x1D6FC}=0,2$ admit an analytical solution with $f$ represented by a parabola and a straight line respectively.

It is possible to extend the self-similar solution to $\unicode[STIX]{x1D705}>0$ with the following expansion in the term $\unicode[STIX]{x1D705}/|\unicode[STIX]{x2202}h/\unicode[STIX]{x2202}x|$ (see, e.g., Hogg, Ungarish & Huppert Reference Hogg, Ungarish and Huppert2000; Sachdev Reference Sachdev2000).

Upon assuming that $\unicode[STIX]{x1D705}/|\unicode[STIX]{x2202}h/\unicode[STIX]{x2202}x|$ is a small quantity, (3.3) becomes

(3.17) $$\begin{eqnarray}{\displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}t}}=-{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}}\left[h\left|{\displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}x}}\right|^{1/n}\left(1-{\displaystyle \frac{2n+1}{n(n+1)}}\unicode[STIX]{x1D705}\left|{\displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}x}}\right|^{-1}+{\displaystyle \frac{n+1-2n^{2}}{2n^{2}}}\unicode[STIX]{x1D705}^{2}\left|{\displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}x}}\right|^{-2}+O(\unicode[STIX]{x1D705}^{3})\right)\right].\end{eqnarray}$$

We propose the following expansion in the regime $\unicode[STIX]{x1D70E}\equiv \unicode[STIX]{x1D705}t^{n(2-\unicode[STIX]{x1D6FC})/(n+2)}\ll 1$ :

(3.18) $$\begin{eqnarray}\displaystyle & \displaystyle h=\unicode[STIX]{x1D702}_{N}^{n+1}t^{F_{2}}[f_{0}(\unicode[STIX]{x1D701})+\unicode[STIX]{x1D70E}f_{1}(\unicode[STIX]{x1D701})+\unicode[STIX]{x1D70E}^{2}f_{2}(\unicode[STIX]{x1D701})+\cdots \,], & \displaystyle\end{eqnarray}$$

(3.19) $$\begin{eqnarray}\displaystyle & \displaystyle x=\unicode[STIX]{x1D702}t^{F_{1}}(1+\unicode[STIX]{x1D70E}X_{1}+\unicode[STIX]{x1D70E}^{2}X_{2}+\cdots \,), & \displaystyle\end{eqnarray}$$

where $f_{0}(\unicode[STIX]{x1D701})$ and $\unicode[STIX]{x1D702}_{N}$ are given by the similarity solution for power-law fluids (3.13), and $X_{1},X_{2},\ldots$ are constants to be evaluated.

The variable $\unicode[STIX]{x1D70E}$ is selected in order to guarantee that at the zero order $O(\unicode[STIX]{x1D70E}^{0})$ the yield stress contribution is null and the solution is represented by (3.13). At the first order $O(\unicode[STIX]{x1D70E})$ there is a balance between the terms due to the yield stress and all other terms. The condition $\unicode[STIX]{x1D70E}\ll 1$ requires that $t\ll t_{c}\equiv \unicode[STIX]{x1D705}^{(n+2)/[n(\unicode[STIX]{x1D6FC}-2)]}$ if $\unicode[STIX]{x1D6FC}<2$ and $t\gg t_{c}$ if $\unicode[STIX]{x1D6FC}>2$ ; $t_{c}$ is defined as a ‘critical time’. Figure 3 shows the critical time $t_{c}(\unicode[STIX]{x1D705},n,\unicode[STIX]{x1D6FC})$ for $\unicode[STIX]{x1D705}=0.01$ and for different $n$ and $\unicode[STIX]{x1D6FC}$ . The critical time becomes infinite for $\unicode[STIX]{x1D6FC}=2$ ; notably, this coincides with the case that permits a self-similar solution without the necessity of an expansion. It is seen that $\unicode[STIX]{x2202}t_{c}/\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}>0$ for any $n$ . In addition, $\unicode[STIX]{x2202}t_{c}/\unicode[STIX]{x2202}n>0$ if $\unicode[STIX]{x1D6FC}>2$ and $\unicode[STIX]{x2202}t_{c}/\unicode[STIX]{x2202}n<0$ if $\unicode[STIX]{x1D6FC}<2$ ; hence, the domain of validity of the expansion is extended as the fluid becomes more shear thinning.

Figure 3. Dimensionless critical time for $\unicode[STIX]{x1D705}=0.01$ as a function of the fluid behaviour index $n$ and of $\unicode[STIX]{x1D6FC}$ . The vertical dashed line at $\unicode[STIX]{x1D6FC}=2$ is the asymptote; the hatched area indicates the domain where the condition $\unicode[STIX]{x1D70E}<1$ is satisfied for a fluid with $n=1.2$ .

The previous analysis implies that for a gravity current of a power-law fluid, all terms in the evolution equations evolve at a common rate (or, equivalently, a unique velocity scale exists). In contrast, for an HB fluid, the term arising due to the yield stress evolves at a different rate from other terms (i.e. it introduces a second velocity scale, for a given common length scale). In order to obtain an expansion to solve this equation, it should be noted that the correction is achieved at first order by computing a second function that evolves with a rate equal to the new one imposed by the yield stress term. The nonlinearity of the problem requires an increasing number of such terms in the series to improve the accuracy in the balance when extended to higher orders.

However, since $\unicode[STIX]{x1D705}$ appears always in powers of $\unicode[STIX]{x1D705}|\unicode[STIX]{x2202}h/\unicode[STIX]{x2202}x|^{-1}$ , we expect a reduction in the accuracy of the asymptotic solution for increasing time, if $\overline{|\unicode[STIX]{x2202}h/\unicode[STIX]{x2202}x|}\sim t^{(\unicode[STIX]{x1D6FC}-2)n/(n+2)}$ decreases in time, which happens for $\unicode[STIX]{x1D6FC}<2$ . Conversely, for $\unicode[STIX]{x1D6FC}>2$ , the asymptotic expansion becomes more accurate, but the current evolves to an increasing steepness which renders the thin-current assumption asymptotically invalid. The expansion will be uniform in $x$ as long as $|\unicode[STIX]{x2202}h/\unicode[STIX]{x2202}x|$ does not approach zero anywhere. In particular, the expansion will be non-uniform if $\unicode[STIX]{x1D6FC}=0$ , since in this case the current becomes flat at $x=0$ .

By inserting (3.18)–(3.19) in (3.17) and balancing the terms of equal power in $\unicode[STIX]{x1D70E}$ , at $O(\unicode[STIX]{x1D70E}^{0})$ we recover the fundamental balance, with $f_{0}\equiv f$ , and $f$ represented by (3.13).

At $O(\unicode[STIX]{x1D70E})$ , equation (3.17) becomes

(3.20) $$\begin{eqnarray}\left(f_{1}|f_{0}^{\prime }|^{1/n}-{\displaystyle \frac{1}{n}}f_{0}|f_{0}^{\prime }|^{1/n-1}f_{1}^{\prime }\right)^{\prime }+F_{2}f_{1}-F_{1}\unicode[STIX]{x1D701}f_{1}^{\prime }={\displaystyle \frac{1}{\unicode[STIX]{x1D702}_{N}^{n}}}{\displaystyle \frac{2n+1}{n(n+1)}}(f_{0}|f_{0}^{\prime }|^{1/n-1})^{\prime },\end{eqnarray}$$

which is an inhomogeneous linear ordinary differential equation for the unknown function $f_{1}$ , with a forcing term modulated by the fundamental solution $f_{0}$ . The numerical integration of equation (3.20) requires two boundary conditions for $\unicode[STIX]{x1D701}\rightarrow 1$ , obtained again by expanding in series near the front of the current. Assuming that $f_{1}\approx a_{1}(1-\unicode[STIX]{x1D701})^{b}$ (with $f_{0}\approx F_{2}^{n}(1-\unicode[STIX]{x1D701})$ , see (3.16)), substituting in (3.20) and equating the lower-order terms (corresponding to $b=1$ ) yields

(3.21a,b ) $$\begin{eqnarray}a_{1}={\displaystyle \frac{F_{2}}{(F_{2}-F_{1}n+F_{2}n)}}{\displaystyle \frac{2n+1}{\unicode[STIX]{x1D702}_{N}^{n}(n+1)}},\quad b=1.\end{eqnarray}$$

Hence, the function can be approximated by $f_{1}\approx a_{1}(1-\unicode[STIX]{x1D701})$ for $\unicode[STIX]{x1D701}\rightarrow 1$ , and the boundary conditions for equation (3.20) are

(3.22a,b ) $$\begin{eqnarray}\left.f_{1}\right|_{\unicode[STIX]{x1D701}\rightarrow 1-\unicode[STIX]{x1D716}}=a_{1}\unicode[STIX]{x1D716},\quad f_{1}^{\prime }|_{\unicode[STIX]{x1D701}\rightarrow 1-\unicode[STIX]{x1D716}}=-a_{1},\end{eqnarray}$$

where $\unicode[STIX]{x1D716}$ is a small quantity.

Even though the two functions $f_{0}$ and $f_{1}$ are defined in the domain $0\leqslant \unicode[STIX]{x1D701}\leqslant 1$ , the nose of the current is $\unicode[STIX]{x1D701}_{N}=1+\unicode[STIX]{x1D70E}X_{1}+\cdots \,,$ which is expected to be smaller than unity since the additional resistance supplied by the yield stress reduces the propagation rate compared with the case $\unicode[STIX]{x1D705}=0$ . The integral constraint given by (3.4) becomes

(3.23) $$\begin{eqnarray}\unicode[STIX]{x1D702}_{N}^{n+2}(1+\unicode[STIX]{x1D70E}X_{1}+\cdots \,)\int _{0}^{1+\unicode[STIX]{x1D70E}X_{1}+\cdots }[f_{0}(\unicode[STIX]{x1D701})+\unicode[STIX]{x1D70E}f_{1}(\unicode[STIX]{x1D701})]\,\text{d}\unicode[STIX]{x1D701}=1,\end{eqnarray}$$

which at $O(\unicode[STIX]{x1D70E})$ yields

(3.24) $$\begin{eqnarray}X_{1}=-{\displaystyle \frac{\displaystyle \int _{0}^{1}f_{1}(\unicode[STIX]{x1D701})\,\text{d}\unicode[STIX]{x1D701}}{\displaystyle \int _{0}^{1}f_{0}(\unicode[STIX]{x1D701})\,\text{d}\unicode[STIX]{x1D701}}}.\end{eqnarray}$$

Figure 4 shows the correction to the front position for waning inflow rate ( $\unicode[STIX]{x1D6FC}=0.5$ ), constant inflow rate ( $\unicode[STIX]{x1D6FC}=1$ ), waxing inflow rate ( $\unicode[STIX]{x1D6FC}=1.5$ ) and very waxing inflow rate ( $\unicode[STIX]{x1D6FC}=2.5$ ). The case $\unicode[STIX]{x1D6FC}=2$ is not shown since the correction is null. The smallest correction is for waxing inflow rates, with minimum effects for shear-thickening fluids, while for low values of $\unicode[STIX]{x1D6FC}$ , the corrections are minor for shear-thinning fluids. The first-order correction term may be valid for a limited time, as shown in figure 4(a) for shear-thickening and Newtonian fluids with $\unicode[STIX]{x1D705}=0.05$ ; a divergence from the first-order approximation solution appears for $t<10$ and $t<20$ respectively. Since the critical times are $t_{c}\approx 100$ and $t_{c}\approx 400$ for the two cases, we conclude that the limiting factor is the number of terms in the expansion. An extension of the range of validity can be achieved by increasing the number of these terms (see, e.g., Hogg et al. Reference Hogg, Ungarish and Huppert2000).

Figure 4. The effects of the first-order correction on the front position, for a Hele-Shaw cell: (a) $\unicode[STIX]{x1D6FC}=0.5$ , (b) $\unicode[STIX]{x1D6FC}=1$ , (c) $\unicode[STIX]{x1D6FC}=1.5$ , (d) $\unicode[STIX]{x1D6FC}=2.5$ . The thick, mid and thin curves refer to $n=1.5$ , $n=1$ and $n=0.5$ respectively. The continuous, dashed and dot-dashed curves refer to $\unicode[STIX]{x1D705}=0$ (power-law fluid), $\unicode[STIX]{x1D705}=0.01$ and $\unicode[STIX]{x1D705}=0.05$ respectively. The time decreasing curves, in grey, are unphysical. Variables are non-dimensional.

Figure 5 shows the profiles of the current at $t=5$ for constant volume ( $\unicode[STIX]{x1D6FC}=0$ ) and constant inflow rate ( $\unicode[STIX]{x1D6FC}=1$ ). The case $\unicode[STIX]{x1D6FC}=0,\unicode[STIX]{x1D705}=0$ has a closed-form solution (Ciriello et al. Reference Ciriello, Longo, Chiapponi and Di Federico2016) and is a parabola for a Newtonian fluid ( $n=1$ ). The presence of yield stress reduces the front position and increases the average steepness of the profile, without other significant variations as long as $\unicode[STIX]{x1D705}$ is a small quantity.

Figure 5. The effects of the first-order correction on the current profiles at time $t=5$ , for a Hele-Shaw cell: (a) $\unicode[STIX]{x1D6FC}=0$ , (b) $\unicode[STIX]{x1D6FC}=1$ . The thick, mid and thin curves refer to $n=1.5$ , $n=1$ and $n=0.5$ respectively. The continuous and dashed curves refer to $\unicode[STIX]{x1D705}=0$ (power-law fluid) and $\unicode[STIX]{x1D705}=0.01$ respectively. Variables are non-dimensional.

4 Two-dimensional flow in a porous medium

The case of flow through a porous medium requires the formulation of the equivalent Darcy’s law for an HB fluid, which may be written as (Chevalier et al. Reference Chevalier, Chevalier, Clain, Dupla, Canou, Rodts and Coussot2013)

(4.1) $$\begin{eqnarray}d\unicode[STIX]{x1D735}p=\unicode[STIX]{x1D712}\unicode[STIX]{x1D70F}_{p}+\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D707}_{0}\left({\displaystyle \frac{\bar{u}}{d}}\right)^{n},\end{eqnarray}$$

where $d$ is the diameter of grains, $\unicode[STIX]{x1D735}p$ is the pressure gradient, $\unicode[STIX]{x1D70F}_{p}$ is the yield stress, $\unicode[STIX]{x1D707}_{0}$ is the consistency index, $n$ is the flow behaviour index, $\bar{u}$ is the Darcian velocity, and $\unicode[STIX]{x1D712}$ and $\unicode[STIX]{x1D6FD}$ are coefficients. The coefficient $\unicode[STIX]{x1D712}$ is governed by the maximum width of the widest path of the flowing current while the coefficient $\unicode[STIX]{x1D6FD}$ depends on pore size distribution and structure. Their values should, in general, be determined experimentally, and theoretically they are related to the distribution of the second invariant of the strain tensor (Chevalier et al. Reference Chevalier, Rodts, Chateau, Chevalier and Coussot2014). Here, a pragmatic approach is adopted, and the values reported in Chevalier et al. (Reference Chevalier, Chevalier, Clain, Dupla, Canou, Rodts and Coussot2013), $\unicode[STIX]{x1D712}=5.5$ and $\unicode[STIX]{x1D6FD}=85$ , are used; the diameter of the glass beads employed in our experiments falls in the range adopted in their experiments (from 0.26 to 2 mm). Equation (4.1) indicates that the flow is possible only if $|\unicode[STIX]{x2202}p/\unicode[STIX]{x2202}x|>\unicode[STIX]{x1D712}\unicode[STIX]{x1D70F}_{p}/d$ ; otherwise a plug is formed. Indeed, the experiments indicate that percolation takes place even at a very low pressure gradient, with a progressive increment of the flow rate for increasing pressure drop. Chevalier et al. (Reference Chevalier, Rodts, Chateau, Chevalier and Coussot2014) have shown that even at very low Darcian velocity values, the region of fluid at rest is negligible and the velocity density distribution is similar to that obtained for a Newtonian fluid.

Under the relevant shallow water approximation, the pressure gradient becomes $\unicode[STIX]{x2202}p/\unicode[STIX]{x2202}x=\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}\,g\,(\unicode[STIX]{x2202}h/\unicode[STIX]{x2202}x)$ ; inverting (4.1), and under the constraint $|\unicode[STIX]{x2202}h/\unicode[STIX]{x2202}x|>\unicode[STIX]{x1D705}_{p}$ , the average velocity is equal to

(4.2) $$\begin{eqnarray}\bar{u}=-\text{sgn}\left({\displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}x}}\right)d^{(n+1)/n}\left({\displaystyle \frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}\,g}{85\unicode[STIX]{x1D707}_{0}}}\right)^{1/n}\left(1-\unicode[STIX]{x1D705}_{p}\left|{\displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}x}}\right|^{-1}\right)^{1/n}\left|{\displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}x}}\right|^{1/n},\end{eqnarray}$$

where $\unicode[STIX]{x1D705}_{p}=5.5\unicode[STIX]{x1D70F}_{p}/(d\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}\,g)$ . Insertion of the average velocity in the local mass conservation yields for $\unicode[STIX]{x2202}h/\unicode[STIX]{x2202}x<0$

(4.3) $$\begin{eqnarray}{\displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}t}}=-{\displaystyle \frac{d^{(n+1)/n}}{\unicode[STIX]{x1D719}}}\left(\frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}\,g}{85\unicode[STIX]{x1D707}_{0}}\right)^{1/n}{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}}\left[h\left|{\displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}x}}\right|^{1/n}\left(1-\unicode[STIX]{x1D705}_{p}\left|{\displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}x}}\right|^{-1}\right)^{1/n}\right],\end{eqnarray}$$

where $\unicode[STIX]{x1D719}$ is the porosity. Assuming again a current with time variable volume, the integral mass conservation reads

(4.4) $$\begin{eqnarray}L_{y}\unicode[STIX]{x1D719}\int _{0}^{\infty }h(x,t)\,\text{d}x=Qt^{\unicode[STIX]{x1D6FC}}.\end{eqnarray}$$

Upon introducing the velocity and length scales

(4.5) $$\begin{eqnarray}\displaystyle & \displaystyle u^{\ast }={\displaystyle \frac{d^{(n+1)/n}}{\unicode[STIX]{x1D719}}}\left(\frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}\,g}{85\unicode[STIX]{x1D707}_{0}}\right)^{1/n}, & \displaystyle\end{eqnarray}$$

(4.6) $$\begin{eqnarray}\displaystyle & \displaystyle x^{\ast }={\displaystyle \frac{d^{(n+1)/n}}{\unicode[STIX]{x1D719}}}\left({\displaystyle \frac{Q\unicode[STIX]{x1D719}}{L_{y}d^{(2n+2)/n}}}\right)^{1/(2-\unicode[STIX]{x1D6FC})}\left({\displaystyle \frac{85\unicode[STIX]{x1D707}_{0}}{\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}\,g}}\right)^{\unicode[STIX]{x1D6FC}/[n(2-\unicode[STIX]{x1D6FC})]}, & \displaystyle\end{eqnarray}$$

with the resulting time scale

(4.7) $$\begin{eqnarray}t^{\ast }={\displaystyle \frac{x^{\ast }}{u^{\ast }}}=\left({\displaystyle \frac{Q\unicode[STIX]{x1D719}}{L_{y}d^{(2n+2)/n}}}\right)^{1/(2-\unicode[STIX]{x1D6FC})}\left({\displaystyle \frac{85\unicode[STIX]{x1D707}_{0}}{\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}\,g}}\right)^{2/[n(2-\unicode[STIX]{x1D6FC})]},\end{eqnarray}$$

equations (4.3)–(4.4) may be written in non-dimensional form as

(4.8) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}t}}=-{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}}\left[h\left|{\displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}x}}\right|^{1/n}\left(1-\unicode[STIX]{x1D705}_{p}\left|{\displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}x}}\right|^{-1}\right)^{1/n}\right], & \displaystyle\end{eqnarray}$$

(4.9) $$\begin{eqnarray}\displaystyle & \displaystyle \int _{0}^{\infty }h(x,t)\,\text{d}x=t^{\unicode[STIX]{x1D6FC}}. & \displaystyle\end{eqnarray}$$

As for the flow in the fracture, equations (4.8)–(4.9) admit a self-similar solution only for $\unicode[STIX]{x1D6FC}=2$ , which breaks down the scales in (4.6)–(4.7). By defining an arbitrary time scale $t^{\ast }$ , the new velocity scale $u^{\ast }=(Q/L_{y}\unicode[STIX]{x1D719})^{1/2}$ and the coefficient

(4.10) $$\begin{eqnarray}\unicode[STIX]{x1D6FF}_{p}=\left({\displaystyle \frac{L_{y}d^{(2n+2)/n}}{Q\unicode[STIX]{x1D719}}}\right)^{1/2}\left({\displaystyle \frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}\,g}{85\unicode[STIX]{x1D707}_{0}}}\right)^{1/n},\end{eqnarray}$$

which is the ratio between the velocity scale (4.5) and the new velocity scale, (4.3) may be written in non-dimensional form (the tilde is dropped) as

(4.11) $$\begin{eqnarray}{\displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}t}}=-\unicode[STIX]{x1D6FF}_{p}{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}}\left[h\left|{\displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}x}}\right|^{1/n}\left(1-\unicode[STIX]{x1D705}_{p}\left|{\displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}x}}\right|^{-1}\right)^{1/n}\right],\end{eqnarray}$$

while (4.4) becomes

(4.12) $$\begin{eqnarray}\int _{0}^{\infty }h(x,t)\,\text{d}x=t^{2}.\end{eqnarray}$$

The self-similar solution is again $h=tf(\unicode[STIX]{x1D702}),\unicode[STIX]{x1D702}=x/t$ , and substitution into (4.11) gives

(4.13) $$\begin{eqnarray}f-\unicode[STIX]{x1D702}f^{\prime }=-\unicode[STIX]{x1D6FF}_{p}\left[f|f^{\prime }|^{1/n}\left(1-{\displaystyle \frac{\unicode[STIX]{x1D705}_{p}}{|f^{\prime }|}}\right)^{1/n}\right]^{\prime }.\end{eqnarray}$$

The system of equations (4.11)–(4.12) admits a solution $f(\unicode[STIX]{x1D702})=A_{p}(\unicode[STIX]{x1D702}_{ep}-\unicode[STIX]{x1D702})$ , which upon substitution yields

(4.14a,b ) $$\begin{eqnarray}{\displaystyle \frac{2}{A_{p}}}=\unicode[STIX]{x1D6FF}_{p}^{2}(A_{p}-\unicode[STIX]{x1D705}_{p})^{2/n},\quad A_{p}\unicode[STIX]{x1D702}_{ep}^{2}=2.\end{eqnarray}$$

For given values of $\unicode[STIX]{x1D705}_{p}$ and $\unicode[STIX]{x1D6FF}_{p}$ , it is possible to solve the first equation numerically in the unknown $A_{p}$ and then to compute $\unicode[STIX]{x1D702}_{ep}$ . The condition for flow requires that $A_{P}>\unicode[STIX]{x1D705}_{p}$ . The general case $\unicode[STIX]{x1D6FC}\neq 2$ can be treated with an expansion similar to that adopted for the flow in a Hele-Shaw cell.

There are many conceptual and formal similarities between the equations arising in the description of Hele-Shaw and 2D porous flows of an HB fluid, while the main point of difference is the treatment of the plug region. While in a Hele-Shaw flow the presence of a plug region is an explicit part of the model (possibly with wall slip), for a porous flow the model predicts the cessation of the flow below a threshold within the entire body of the granular medium. During the experiments described in the following, it was noted that below a threshold value of the pressure gradient ( $\unicode[STIX]{x2202}h/\unicode[STIX]{x2202}x$ in our approximation), a percolation develops and the flow never completely stops. This behaviour suggests that a biviscous model, able to smooth the transition from pre- to postyield behaviour, could better interpret the experiments.

5 The experiments

In order to validate the theoretical model, two series of experiments were conducted (i) in a Hele-Shaw cell with a small gap, simulating a fracture, and (ii) in the same cell with a larger gap and filled with glass beads of uniform size, reproducing a 2D porous medium. Figure 6 shows the experimental device and two different snapshots. The $75~\text{cm}$ long cell was made of two parallel plates of transparent plastic, the gap width between which could be varied as necessary. In order to limit slip for the experiments without glass beads, a commercial transparent antislip tape was used to line the inside of the plates. Fluid was injected with a syringe pump for the fracture tests, and with a vane pump controlled by an inverter for the porous flow tests, requiring higher flow rates.

Figure 6. A sketch of the experimental rectangular channel. (a) Front view, (b) side view, (c) a snapshot of the channel during experiment B1 (the shaded area is the advancing current) and (d) a snapshot of the channel filled with glass beads during experiment A1 (the dark area is the advancing current and the grey area is the porous medium not yet reached by the current).

The HB fluids were obtained by mixing deionized water, Carbopol 980 (0.05 %–0.14 %) and ink, with subsequent neutralization by adding NaOH. The mass density was measured with a hydrometer (STV3500/23, Salmoiraghi) or a pycnometer, with an accuracy of 1 %.

The rheologic behaviour of the fluid was analysed in a parallel-plate rheometer (dynamic shear rheometer, Anton Paar Physica MCR 101), conducting several different tests, both static and dynamic, to evaluate the fluid behaviour index, the consistency index and the yield stress. When dealing with non-Newtonian fluids, either Ostwald–de Waele (power law) or HB, the characterization of the correct parameters representing the rheological behaviour of the fluid can be challenging. However, the final aim of the experiments is clear, and that is to verify the proposed model by using independent measurements of both the flow field characteristics (front position and current thickness over time) and the fluid rheometric parameters. It is noteworthy that the models used to describe the fluid rheology and obtain the rheometric parameters do not have a general validity but are an approximation of the fluid behaviour in a limited shear range. Furthermore, none of the flow fields generated in the measuring instruments (mainly rheometers) are perfectly viscometric. The problem of correctly estimating the rheometric parameters is particularly significant for HB fluids, for which yield stress estimation (and definition) is far from trivial. See the review paper by Nguyen & Boger (Reference Nguyen and Boger1992) for a discussion on the topic and a description of the methodology adopted in our laboratory experiments. Therefore, in order to evaluate the accuracy of the measurements of the yield stress, we carried out numerous additional experiments with different methods, broadly classified as ‘direct’ and ‘indirect’ methods; see the further description in appendix A.

The glass beads forming the porous medium had nominal diameters of 3, 4 and 5 mm $\pm 0.1~\text{mm}$ , depending on the test. The profile of the current was detected with either a stills camera (Canon EOS 3D, $3456\times 2304$ pixels) operating at a rate of $0.5~\text{frames}\,\text{s}^{-1}$ or a high-resolution video camera (Canon Legria, $1920\times 1080$ pixels) operating at $25~\text{frames}\,\text{s}^{-1}$ . The images were then postprocessed with a proprietary software package in order to be referenced to a laboratory coordinate system and for the boundary between the (dark) intruding current and the empty cell or the (light) porous medium to be extracted and parameterized. With this set-up, the overall accuracy in detecting the profile of the current was approximately 1 mm, while the uncertainty in measuring timings was negligible. During all experiments, the natural packing prevented any movement of the beads, as demonstrated by the images used for extracting the fluid interface.

Tables 1 and 2 list the parameters for the two sets of experiments (fracture and 2D porous medium) and four series with increasing concentration of Carbopol (labelled from A to D, corresponding to $0.05\,\%$ , $0.08\,\%$ , $0.10\,\%$ and $0.14\,\%$ concentration respectively). Figure 7(a) shows the dimensionless front position of the currents versus time for 13 tests conducted in the fracture with three different HB fluids. The plotting variables have been chosen to collapse all of the experimental data into a single line. In order to separate the results for the three different fluids, the experimental front position was multiplied by 0.5 and 1.5 for experiments A and C respectively. The uncertainty in the model and experimental data was computed following the same procedure as reported in Di Federico et al. (Reference Di Federico, Longo, Chiapponi, Archetti and Ciriello2014), and is represented by the dashed lines and error bars corresponding to plus or minus one standard deviation (STD) for the data. The front propagation is generally linear for all tests, with some discrepancy at early times due to the effects of the injection geometry and the poor adherence of the experiments to model assumptions. We recall that the solution is an ‘intermediate asymptotic’ to the general solution, and is valid for times and distances from the boundaries large enough to forget the details of the boundary (or initial) conditions but far enough from the ultimate asymptotic state of the system. Hence, we expect that also for longer time the experimental results will deviate from the self-similar solution.

Table 1. Parameters for the experiments in the fracture. The injected volume scales with $t^{2}$ .

Table 2. Parameters of the experiments in a porous medium, with volume $\propto t^{\unicode[STIX]{x1D6FC}}$ . The experimental and theoretical front speed is not constant for the last four experiments.

Figure 7. A comparison of theory with experiments for the fracture. (a) The front position $x_{N}/\unicode[STIX]{x1D702}_{e}$ as a function of dimensionless time $t$ for all tests. The three bold lines represent the perfect agreement with theory for the three different fluids used in the experiments. (b) The dimensionless profile of the current at different times for experiment B3. The error bars refer to the profile at $t=66~\text{s}$ and correspond to $\pm$ one STD, the bold straight line indicates the theoretical profiles and the dashed lines are the confidence limits of the model. For clarity, in both (a) and (b), one point of every three is plotted.

Figure 7(b) shows the shape of the current at different times for experiment B3. The normalized profiles show a fairly good collapse to a single curve, even though the discrepancy between experiments and theory becomes evident for $\unicode[STIX]{x1D702}/\unicode[STIX]{x1D702}_{e}<0.30$ and near the front. This is due to the disturbances at the inlet (the inflow is located near the bottom in the experiments, while it is assumed to be evenly distributed along the vertical in the theory), to non-negligible vertical velocities near the inlet (see Longo & Di Federico Reference Longo and Di Federico2014) and to the bottom stress, which becomes relevant at the tip of the current. The experimental profiles are similar for all other tests.

Figure 8(a) shows the front position for the six experiments conducted in a porous medium with three different HB fluids. The profile was corrected for capillary effects following the procedure outlined in Longo et al. (Reference Longo, Di Federico, Chiapponi and Archetti2013b ). The front velocity shows a fairly good agreement with the theoretical constant velocity, even though asymptotically there are larger discrepancies due to the geometry of the current, with a very thin nose affected by the bottom boundary effects. The experimental profiles of the current at different times are shown in figure 8(b) for experiment A2. The profiles collapse to a single curve, which is significantly affected by the disturbances at the inlet. The experimental points are within the confidence limits of the theoretical models and show a clear linearity for $\unicode[STIX]{x1D702}>0.5\unicode[STIX]{x1D702}_{ep}$ . In order to also check the model for the asymptotic solution with $\unicode[STIX]{x1D6FC}\neq 2$ , some experiments were designed for constant ( $\unicode[STIX]{x1D6FC}=1$ ) and waning ( $\unicode[STIX]{x1D6FC}=0.6$ ) inflow rate, with two different HB fluids and three different bead diameters (see the last four experiments in table 2). Figure 9 shows the experimental front position (symbols) and the first-order expansion of the self-similar solution (solid curves). The overlap in each case is satisfactory, in particular for experiments at constant inflow rate.

Figure 8. A comparison of theory with experiments for 2D porous flow (rectangular channel filled with glass beads). For caption, see figure 7. (b) Refers to experiment A2.

Figure 9. A comparison of theory (solid curves) with experiments (symbols) for the 2D porous medium, showing $x_{N}$ as a function of dimensionless time $t$ for all tests. The injected volume scales with $t^{\unicode[STIX]{x1D6FC}}$ . The experimental parameters are listed in table 2.