2.1. Constitutive law

A generalized tensorial extension of the Herschel–Bulkley constitutive law within the framework of the shallow-flow asymptotic expansion has been proposed in our recent work (Denisenko et al. Reference Denisenko, Richard and Chambon2025). Namely, the extra-stress tensor $\boldsymbol{\tau }$ is expressed as

(2.1) \begin{align} \tau _{\textit{ij}} = {\tau }_{\textit{ij}}^Y \displaystyle + K{|\dot {\boldsymbol{\gamma }}|^{n-1}}\dot { \gamma }_{\textit{ij}} \qquad & |\dot {\boldsymbol{\gamma }}|\gt 0, \end{align}

(2.2) \begin{align} \dot {{\gamma }}_{\textit{ij}} = 0 \qquad & |\boldsymbol \tau | \leqslant \tau _c. \end{align}

The parameters $K$ and $n$ are the consistency and the power index of the material, respectively. The tensor $\boldsymbol{\tau }^Y$ is called the yield-stress tensor. The norm is defined as $ |\boldsymbol{T}| = (0.5\boldsymbol{T}_{\textit{ij}}\boldsymbol{T}_{\textit{ij}} )^{0.5}$ for any second-order tensor. Hereinafter, for $|\boldsymbol \tau | \gt \tau _c$ , we split the extra-stress in its yield-stress and viscous contributions as

(2.3) \begin{align} \tau _{\textit{ij}}=\tau _{\textit{ij}}^Y+\tau _{\textit{ij}}^{v}, \quad \tau _{\textit{ij}}^{v}=K\dot {\gamma }_{\textit{ij}}|\dot {{\boldsymbol \gamma }}|^{n-1}. \end{align}

Unlike in the classical tensorial extension of the Herschel–Bulkley rheology, here the yield-stress tensor is not assumed to be aligned with the strain-rate tensor via the expression $\tau _{\textit{ij}}^Y = \tau _c { {\dot {\gamma }_{\textit{ij}}}}/|{\boldsymbol{\dot \gamma }}|$ . We instead define the yield-stress tensor $\boldsymbol{\tau }^{Y}$ based on the following criteria.

1. (i) The norm of $\boldsymbol{\tau ^Y}$ is equal to the yield stress $\tau _c$ in the yielded regime: $|\boldsymbol{\tau ^Y}|=\tau _c$ . This corresponds to the von Mises yielding criterion.

2. (ii) The trace of $\boldsymbol{\tau ^Y}$ is zero in relation to the incompressibility assumption: $\mathrm{tr}\, \boldsymbol{\tau ^Y}=0$ .

3. (iii) There are no normal stress differences in simple shear, such that $\boldsymbol{\tau ^Y}$ is aligned with $\boldsymbol{\dot {\gamma }}/|\boldsymbol{\dot {\gamma }}|$ in the particular case of simple shear flows. Although questionable, this assumption is made to obtain the same leading-order solution as in models based on the classical formulation of the rheology.

4. (iv) The existence of a normal stress difference in the pseudoplug at leading order (see below), and all terms originating from these normal stresses in the asymptotic expansions at $O(\varepsilon )$ , are attributed to the yield-stress tensor and not to the viscous-stress tensor $\tau ^{v}=K\dot {\gamma }|\dot {{\boldsymbol \gamma }}|^{n-1}$ . All other terms in the asymptotic expansions are attributed to the viscous-stress tensor.

The constitutive law (2.1) together with these four conditions can be viewed as an alternative extension of the one-dimensional Herschel–Bulkley law for two-dimensional shallow flows, since the classical scalar expression $\tau _{xz}=\tau _c+\dot \gamma ^n$ is recovered in simple-shear flows. The fourth condition here can be interpreted as a closure assumption ensuring that the normal stress differences associated with the perfectly rigid regime do not influence the fluid-like dynamics. It avoids the spurious destabilizing effects of plasticity at large slope angles that are observed in the model based on the classical tensorial law (Denisenko et al. Reference Denisenko, Richard and Chambon2023). It is important to emphasize that the four conditions formulated above are meaningful only in the shallow-flow regime, where a dominant direction can be identified and the flow is close to one-dimensional. Although these conditions are not sufficient to fully determine the yield-stress tensor in a general flow, they are sufficient to work out a thin-layer asymptotic expansion up to $O(\varepsilon )$ .

The conditions on $\tau _{\textit{ij}}^{Y}$ formulated above were first proposed in Denisenko et al. (Reference Denisenko, Richard and Chambon2023) to avoid the singularity of the yield-stress tensor for $|\boldsymbol{\dot {\gamma }}|=0$ . This allowed us to derive a consistent three-equation shallow-flow model for Bingham fluids with a smooth velocity field and a physically correct instability threshold. Later these conditions were used in the derivation of a smooth asymptotic expansion for Herschel–Bulkley fluids (Denisenko et al. Reference Denisenko, Richard and Chambon2025). Here we build on these previous developments to derive a consistent three-equation shallow-flow model in the Herschel–Bulkley case.

2.2. Governing equations

The continuity equation (mass conservation balance) reads

(2.4) \begin{align} \frac {\partial u}{\partial x} + \frac {\partial w}{\partial y}=0. \end{align}

The Cauchy momentum equations in the $x$ and $z$ directions read

(2.5) \begin{align} \rho \left ( \frac {\partial u}{\partial t} + u\frac {\partial u}{\partial x} + w\frac {\partial u}{\partial z} \right ) = -\frac {\partial p}{\partial x} + \rho g \sin {\theta } +\frac {\partial \tau _{xz}}{\partial z} +\frac {\partial \tau _{xx}}{\partial x}, \end{align}

(2.6) \begin{align} \rho \left ( \frac {\partial w}{\partial t} + u \frac {\partial w}{\partial x} + w \frac {\partial w}{\partial z} \right ) = -\frac {\partial p}{\partial z} -\rho g \cos {\theta } +\frac {\partial \tau _{xz}}{\partial x} +\frac {\partial \tau _{xx}}{\partial z}. \end{align}

At the bottom, no-penetration and no-slip boundary conditions are assumed:

(2.7) \begin{align} u_{|_{z=0}}=w_{|_{z=0}}=0. \end{align}

At the free surface, the kinematic boundary condition is given by

(2.8) \begin{align} \frac {\partial h}{\partial t} + u_{|_{z=h(x)}}\frac {\partial h}{\partial x} =w_{|_{z=h(x)}}, \end{align}

while the dynamic boundary conditions are

(2.9) \begin{align} \left [ \displaystyle 1- \left ( \frac {\partial h}{\partial x} \right )^2 \right ]\tau _{xz|_{z=h(x)}} = 2 \displaystyle \left (\frac {\partial h}{\partial x}\right ) \tau _{xx|_{z=h(x)}},\end{align}

(2.10) \begin{align} \left [ \displaystyle 1- \left ( \frac {\partial h}{\partial x} \right )^2 \right ] p_{|_{z=h(x)}} = - \left [1+ \displaystyle \left (\frac {\partial h}{\partial x} \right )^2 \right ] \tau _{xx|_{z=h(x)}}.\end{align}

Note that capillarity is neglected in this study. The atmospheric pressure is taken constant and equal to zero.

2.3. Shallow-flow scaling

Let us denote by $h_0$ the characteristic length in $z$ direction, by $l_0$ the characteristic length in $x$ direction and by $u_0$ the characteristic velocity in $x$ direction. The shallow-flow hypothesis implies that the aspect ratio $\varepsilon = h_0/l_0$ is small. The relevant dimensionless groups of the problem are the Reynolds number $ \textit{Re}$ , the Froude number $ \textit{Fr}$ and the Bingham number $ \textit{Bi}$ :

(2.11) \begin{align} \textit{Re} = \frac {\rho u_0^{2-n}h_0^{n}}{K}, \qquad Fr = \frac {u_0}{\sqrt {gh_0\cos \theta }}, \qquad \textit{Bi} = \frac {\tau _c}{K} \left (\frac {h_0}{u_0}\right )^n. \qquad \end{align}

In this study, these parameters and the slope angle $\theta$ are assumed to be $O(1)$ with respect to $\varepsilon$ . To rewrite the problem (2.4)–(2.10) into a dimensionless form, the flow variables are expressed as

(2.12) \begin{align} x = l_0{x^\prime}; \quad z = h_0{z^\prime}; \quad u = u_0{u^\prime}; \quad w = \varepsilon u_0 {w^\prime}; \quad t = \frac {l_0}{u_0} {t^\prime}; \quad h = u_0 {h^\prime}; \nonumber \\ p = \rho g h_0\cos \theta {p^\prime}; \quad \tau _{\textit{ij}} = \tau _c {\tau ^\prime}_{\textit{ij}}; \quad |\boldsymbol \tau | =\tau _c|{\boldsymbol \tau }^\prime|; \quad |\dot {\boldsymbol{\gamma }}| = \frac {u_0}{h_0}|\dot {{\boldsymbol{\gamma }}}^\prime|. \end{align}

For convenience, we shall further omit the primes. The continuity (2.4) keeps the form

(2.13) \begin{align} \frac {\partial u}{\partial x} + \frac {\partial w}{\partial y}=0, \end{align}

while the scaled momentum (2.5)–(2.6) take the following form:

(2.14) \begin{align} \varepsilon \left (\frac {\partial u}{\partial t} + {u}\frac {\partial u}{\partial x} + {w}\frac {\partial u}{\partial z}\right ) = - \frac {\varepsilon }{Fr^2} \frac {\partial p}{\partial x} + \frac {\lambda }{ \textit{Re}} + \frac {\textit{Bi}}{ \textit{Re}} \frac {\partial \tau _{xz}}{\partial z} + \frac {\varepsilon \textit{Bi}}{\textit{Re}}\frac {\partial \tau _{xx}}{\partial x},\end{align}

(2.15) \begin{align} \varepsilon ^2 \left ( \frac {\partial w}{\partial t} + {u} \frac {\partial w}{\partial x} + {w} \frac {\partial w}{\partial z} \right ) = -\frac {1}{Fr^2} \left ( 1+\frac {\partial p}{\partial z} \right ) + \frac {\varepsilon \textit{Bi}}{\textit{Re}} \frac {\partial \tau _{xz}}{\partial x} + \frac {\textit{Bi}}{\textit{Re}}\frac {\partial \tau _{zz}}{\partial z},\end{align}

with the driving parameter $\lambda$ given by

(2.16) \begin{align} \lambda = \frac {\rho g h_0^{n+1} \sin {\theta }}{K{u}_0^n}=\frac {\textit{Re}}{Fr^2}\tan {\theta }. \end{align}

In the scaled variables, the yielding criterion reduces to $|\boldsymbol \tau |=1$ . For $|\boldsymbol \tau |\gt 1$ , the expressions for the viscous and yield-stress components of (2.3) are transformed to

(2.17) \begin{align} {\tau }_{xx} = {\tau }_{xx}^Y + {\tau }_{xx}^{v}, \qquad {\tau }_{xx}^{v} = \frac {2\varepsilon }{\textit{Bi}} \frac {\partial u}{\partial x}| {\boldsymbol{{\dot \gamma }}}|^{n-1},\end{align}

(2.18) \begin{align} {\tau }_{xz} = {\tau }_{xz}^Y + {\tau }_{xz}^{v}, \qquad {\tau }_{xz}^{v} = \frac {1}{\textit{Bi}}\left (\frac {\partial u}{\partial z} +\varepsilon ^2 \frac {\partial w}{\partial x} \right ) | {\boldsymbol{{\dot \gamma }}}|^{n-1}.\end{align}

The condition on the norm of the dimensionless yield-stress tensor is

(2.19) \begin{align} \left ({\tau }_{xx}^Y\right )^2 + \left ({\tau }_{xz}^Y\right )^2=1. \end{align}

The no-penetration and no-slip conditions (2.7) keep the same form,

(2.20) \begin{align} {u}_{|_{{z}=0}}={w}_{|_{{z}=0}}=0, \end{align}

while the scaled kinematic and dynamic boundary conditions (2.8)–(2.10) are rewritten as

(2.21) \begin{align} \frac {\partial h}{\partial x} + {u}_{|_{{z}={h}({x})}}\frac {\partial h}{\partial x} ={w}_{|_{{z}={h}({x})}},\end{align}

(2.22) \begin{align} \left [ \displaystyle 1- \varepsilon ^2\left ( \frac {\partial h}{\partial x} \right )^2 \right ]{\tau _{xz}}_{|_{{z}={h}({x})}} = 2 \varepsilon \left (\frac {\partial h}{\partial x}\right ) {\tau _{xx}}_{|_{{z}={h}({x})}},\end{align}

(2.23) \begin{align} \left [ 1-\varepsilon ^2 \left ( \frac {\partial h}{\partial x} \right )^2 \right ] {p}_{|_{{z}={h}({x})}} = -\frac {\textit{Bi} Fr^2}{\textit{Re}} \left [1+ \varepsilon ^2\left (\frac {\partial h}{\partial x} \right )^2 \right ] {{\tau }_{xx}}_{|_{{z}={h}({x})}}.\end{align}

Finally, the norms of the dimensionless stress and strain-rate tensors read

(2.24) \begin{align} |{\boldsymbol \tau }| = \sqrt {{\tau }_{xx}^2+{\tau }_{xz}^2}, \qquad |{\boldsymbol{{\dot \gamma }}}| = \sqrt { \left ( {\partial _z {u}} + \varepsilon ^2 {\partial _x {w}} \right )^2 + 4\varepsilon ^2 \left ( {\partial _x {u}} \right )^2 }. \end{align}

4.1. Mass conservation

Averaging the continuity (2.13) and taking into account the kinematic boundary condition (2.21) yields the following mass conservation equation:

(4.3) \begin{align} \frac {\partial h}{\partial t} +\frac {\partial hU}{\partial x}=0. \end{align}

Inserting the leading-order asymptotic expansion (4.2) for $U$ into (4.3) results in a kinematic wave equation:

(4.4) \begin{align} \frac {\partial h}{\partial t} +\lambda ^{\frac {1}{n}} h(h-h_p)^{\frac {1}{n}} \frac {\partial h}{\partial x}=\text{O}(\varepsilon ). \end{align}

This equation, together with the expressions of $u^{(1)}$ given by (A11)–(A17), allows us to write the first-order correction $U^{(1)}$ in the form given by (A18) in Appendix A.

4.2. Momentum balance

The momentum balance (2.14) is averaged over the depth, using the no-slip condition (2.20), the kinematic boundary condition (2.21) and the dynamic boundary conditions (2.22) and (2.23). Employing the leading-order expressions for the stresses (A1)–(A4) and for the pressure (A6), yields

(4.5) \begin{align} \frac {\partial }{\partial t} \left ( hU \right ) +\frac {\partial }{\partial x} \left ( h\langle u^2 \rangle +\frac {h^2}{2Fr^2} \right )= -\frac {\textit{Bi}\tau _{xz}^{(1)}(0)}{ \textit{Re}} + O(\varepsilon ), \end{align}

where $\tau _{xz}^{(1)}(0)$ denotes the first-order correction to the shear stress at the bottom $z=0$ .

In most shallow-flow models, the quadratic term $\langle u^2 \rangle$ is replaced by $U^2$ , so the mass conservation (4.3) and the momentum balance (4.5) yield a system of two coupled equations for the flow height $h$ and the depth-averaged velocity $U$ . This assumption is justified if the fluid is only slightly sheared. However, in general viscous fluids this is not the case. To treat properly the term $\langle u^2 \rangle$ , we introduce the deviation of the velocity from its averaged value:

(4.6) \begin{align} u^{\ast }(x,z,t) = u(x,z,t) - U(x,t). \end{align}

By definition $\langle u^{\ast } \rangle = 0$ , so that $\langle u^2 \rangle =U^2+\langle u^{\ast 2} \rangle$ in (4.5). The deviation $\langle u^{\ast 2} \rangle$ could be formally estimated by the asymptotic expansion at leading order. This would yield a closed system of two coupled equations for the variables $h$ and $U$ (Fernández-Nieto et al. Reference Fernández-Nieto, Noble and Vila2010; Noble & Vila Reference Noble and Vila2013). However, as shown by Richard et al. (Reference Richard, Gisclon, Ruyer-Quil and Vila2019) for viscous fluids, to guarantee the compatibility of the averaged mass and momentum equations with the depth-averaged energy equation, one should rather consider the deviation $\langle u^{\ast 2} \rangle$ as an independent variable (see also Teshukov Reference Teshukov2007). We adapt here this three-variable approach for Herschel–Bulkley fluids. Instead of $\langle u^{\ast 2} \rangle$ , it is more convenient to introduce the following third variable:

(4.7) \begin{align} \varphi = \frac {\langle {u}^{\ast 2} \rangle }{h(h-h_p)}. \end{align}

This variable $\varphi$ plays the role of an entropy for the system (Richard & Gavrilyuk Reference Richard and Gavrilyuk2012) and is called enstrophy. The negative correction term $-h\,h_p$ in the denominator of (4.7) is required only in the viscoplastic case, in order to prevent a non-physical relaxation of the enstrophy towards zero at $h=h_p$ later in the model.

The expansion of the enstrophy,

(4.8) \begin{align} \varphi = \varphi ^{(0)} + \varepsilon \varphi ^{(1)} + {O(\varepsilon ^2) } ,\end{align}

leads to $\varphi ^{(0)} = \big \langle (u^{\ast (0)} )^2 \big \rangle /h(h-h_p)$ and $\varphi ^{(1)} = 2\langle u^{\ast (0)}u^{\ast (1)} \rangle /h(h-h_p)$ . Thus, at leading order we have

(4.9) \begin{align} \varphi ^{(0)}= \frac {\lambda ^{\frac {2}{n}} h^{\frac {2}{n}}{n^3}}{(2n+1)^2(3n+2)} \left (1-\xi \right )^{\frac {2(1+n)}{n}} \left (1+\frac {3n+2}{(n+1)^2}{n}\xi \right ). \end{align}

Here $\xi =h_p/h$ . At order $O(\varepsilon )$ the expression $\varphi ^{(1)}$ is lengthy and given by (A19) in Appendix A.

Introducing the new variable $\varphi$ , (4.5) is rewritten as

(4.10) \begin{align} \frac {\partial }{\partial t} \left ( hU \right ) +\frac {\partial }{\partial x} \left (hU^2+h^3\varphi +\frac {h^2}{2Fr^2}\right )= -\frac {\textit{Bi}\tau _{xz}^{(1)}(0)}{ \textit{Re}} + \frac {\partial }{\partial x} \left ( h_p h^2\varphi \right ). \end{align}

To derive a model with a well-posed mathematical structure, the source terms on the right-hand side of (4.10) should be expressed as a sum of relaxation terms for $h$ , $U$ and $\varphi$ . Using the relation (4.4), the shear stress correction (A9) and the expression (A18) for $U^{(1)}$ , one can write

(4.11) \begin{align} \frac {\textit{Bi}{\tau _{xz}^{(1)}}{(0)}}{\textit{Re}} = \frac {3 n (n+1) (2 n+1)}{3 n \left [2 n \xi \left (n \xi +1\right )+n+1\right ] \left (1-\xi \right ){}^{\frac {1}{n}}+\zeta ^{\frac {1-n}{n}} (n+1) (2 n+1)}\frac {U^{(1)} }{\lambda ^{\frac {1-n}{n}} h^{\frac {1}{n}} \textit{Re}}\nonumber \\ -\lambda ^{\frac {2}{n}}h^{\frac {2(1+n)}{n}} C_1(n,\xi ,\zeta ) \frac {\partial h}{\partial x}, \end{align}

where the function $C_1$ is given in Appendix B and the dimensionless quantity $\zeta$ is defined as $\zeta = \delta /h$ . The second term on the right-hand side of (4.10) can be expressed at leading order using (4.9):

(4.12) \begin{align} \frac {\partial }{\partial x} \left ( h_p h^2 \varphi \right ) = \lambda ^{\frac {2}{n}}h^{\frac {2(1+n)}{n}} C_2(n,\xi )\frac {\partial h}{\partial x} + O(\varepsilon ). \end{align}

Here the function $C_2$ is given in Appendix B.

We can now express the different terms in (4.11) and (4.12) using forms that relax towards the equilibrium (leading-order) solution. For the first term in (4.11), using the asymptotic expansion $U = U^{(0)} + \varepsilon U^{(1)} + O(\varepsilon ^2)$ and expression (4.2) for $U^{(0)}$ , we can write

(4.13) \begin{align} \lambda h [1-q(n,\xi )]& - \textit{Bi} - \alpha _1(n)\frac { U}{h}\left |\frac {U}{h}\right |^{n-1} =\nonumber \\ &-\frac {\varepsilon }{\lambda ^{\frac {1-n}{n}}h^{\frac {1}{n}}}\left [(2 n+1) (1-\xi )^{n-\frac {1}{n}} \left (\frac {n+1}{n \xi +n+1}\right )^{1-n}\right ]U^{(1)} + O(\varepsilon ^2) \end{align}

where $\alpha _1$ and $q$ are given by

(4.14) \begin{align} \alpha _1(n)= \left (\frac {2 n+1}{n}\right )^n, \quad q(n,\xi )=(1-\xi ) \left [1-\left (\frac {(1-\xi ) (n \xi +n+1)}{n+1}\right )^n\right ]. \end{align}

Accordingly, the quantity $U^{(1)}$ can be written as a relaxation term in $U$ :

(4.15) \begin{align} U^{(1)} = -\frac {\lambda ^{\frac {1-n}{n}} h^{\frac {1}{n}}(1-\xi )^{\frac {1-n^2}{n}} \left ({n+1}\right )^{n-1}}{\varepsilon (2 n+1) ({n \xi +n+1})^{n-1} }\Bigg [ \lambda h [1-q(n,\xi )] - \textit{Bi} - \alpha _1(n)\frac {U}{h}\left |\frac {U}{h}\right |^{n-1} \Bigg ]. \end{align}

The derivative $\partial h/\partial x$ in (4.11) and (4.12) can be expressed consistently as the sum of two relaxation terms in $U$ and $\varphi$ . First, let us consider the following relation obtained from the asymptotic expansions for $U$ and $\varphi$ :

(4.16) \begin{align} \varphi ^{\frac {n-2}{2}}\Bigg [\varphi -\alpha _2(n,\xi )\frac {U^2}{h^2} \Bigg ] = \varepsilon \textit{Re}\lambda ^{\frac {2}{n}}h^{\frac {2(1+n)}{n}} C_3(n,\xi ,\zeta ) \frac {\partial h}{\partial x} + \varepsilon \frac {C_4(n,\xi ,\zeta )}{\lambda ^{\frac {1-n}{n}}h^{\frac {1}{n}}} U^{(1)} +O(\varepsilon ^2). \end{align}

Here the functions $C_3$ and $C_4$ are given in Appendix B, while the function $\alpha _2$ reads

(4.17) \begin{align} \alpha _2(n,\xi )=n\frac {n [(3 n+2) \xi +n+2]+1}{(3 n+2) (n \xi +n+1)^2}. \end{align}

From (4.2) and (4.9), the left-hand side of (4.16) vanishes at leading order, which means that this term has the structure of a relaxation term in $\varphi$ . The factor $\varphi ^{{(n-2)}/{2}}$ on the left-hand side of (4.16) is introduced to avoid the appearance of a factor $(1/\lambda )^{(2-n)/n}$ in front of a source term further in the model, that would lead to a divergence in the case $\theta =0$ . The second term on the right-hand side of (4.16) also has a relaxation structure owing to (4.15). Finally, the integrated momentum (4.10) can thus be rewritten as

(4.18) \begin{align} \frac {\partial }{\partial t} \left ( hU \right ) +\frac {\partial }{\partial x} \left (hU^2+h^3\varphi +\frac {h^2}{2Fr^2}\right ) = \nonumber \\ \frac {1}{\varepsilon \textit{Re}}\left [\lambda h [1-q(n,\xi ) ] - \textit{Bi} - \alpha _1(n)\frac {U}{h}\left |\frac {U}{h}\right |^{n-1} \right ]\beta _1(n,\xi ,\zeta ) \nonumber \\ +\frac {\varphi ^{\frac {n-2}{2}}}{\varepsilon \textit{Re}}\left [\varphi - \alpha _2(n,\xi )\frac {U^2}{h^2} \right ]\beta _2(n,\xi ,\zeta ), \end{align}

where the functions $\beta _1$ and $\beta _2$ are given in Appendix B. The first relaxation term in (4.18) interprets as the balance between gravity, yield stress and viscous friction along the $x$ -axis.

As will be shown later, writing the viscous friction term in the form of a power law in $U/h$ is important to get proper free-surface shapes for steady-state surges (see § 6). The additional friction term $\lambda h \, q(n,\xi )$ is specific to the viscoplastic case (since $q=0$ if $\tau _c=0$ ) and originates from the consistency requirement when expressing the first source term as a sum of weight, yield stress and power-law viscous friction.

The second term on the right-hand side of (4.18) corresponds to a relaxation for the enstrophy. The terms with $\lambda h q$ and $ \textit{Bi}$ in (4.18) formally oppose the flow and should, in general, depend on the sign of $U$ . However, in the present study, we restrict our analysis to flow regimes with positive velocity values only.

4.3. Enstrophy balance equation

The mass and momentum (4.3) and (4.18) involve three unknown variables: $h$ , $U$ and $\varphi$ . A third equation for an enstrophy balance is thus needed to close the problem. To that aim, we first consider the work–energy theorem (or kinetic energy equation), which in dimensional variables reads

(4.19) \begin{align} \frac {\partial }{\partial {t}} \left (\frac {1}{2}\rho {v}^2\right ) +\text{div} \left ( \frac {1}{2}\rho {v}^2{\boldsymbol{{v}} } \right )= \text{div}(\boldsymbol{{\sigma }} \boldsymbol{\cdot }\boldsymbol{ {v} } ) -{{\boldsymbol \sigma }}:{\frac {{\boldsymbol{{\dot \gamma }}}}{2}} +\rho {\boldsymbol g}{\boldsymbol{{v}} }, \end{align}

where $ {\boldsymbol{{\sigma } }} =-{p}\boldsymbol I+\boldsymbol{{\tau }}$ . Introducing the components of the vectors and tensors, we obtain the following expression in dimensionless form:

(4.20) \begin{align} \frac {\partial }{\partial t} \left ( \frac {u^2}{2}\!+\!\varepsilon ^2\frac {w^2}{2} \right ) \!+\!\frac {\partial }{\partial x}\! \left [ u \left ( \frac {u^2}{2} \!+\!\varepsilon ^2\frac {w^2}{2} -\frac {x\tan {\theta }}{\varepsilon Fr^2} \!+\!\frac {z}{Fr^2} \right ) \!+\!\frac {pu}{Fr^2} -\frac {\textit{Bi}}{\textit{Re}} (\tau _{xx}u+\varepsilon \tau _{xz}w) \right ] \nonumber \\ + \frac {\partial }{\partial z} \left [ w \left ( \frac {u^2}{2} +\varepsilon ^2\frac {w^2}{2} -\frac {x\tan {\theta }}{\varepsilon Fr^2} +\frac {z}{Fr^2}\right ) +\frac {pw}{Fr^2} -\frac {\textit{Bi}}{\varepsilon \textit{Re}}\tau _{xz} u -\frac {\textit{Bi}}{\textit{Re}}\tau _{zz} w \right ] \nonumber \\ =-\frac {2\textit{Bi}}{\textit{Re}} \left ( \tau _{xx}^{{Y}} +2\frac {\varepsilon }{\textit{Bi}} \frac {\partial u}{\partial x} \right ) \frac {\partial u}{\partial x} - \frac {\textit{Bi}}{\varepsilon \textit{Re}} \left ( \tau _{xz}^{{Y}} +\frac {1}{\textit{Bi}} \left ( \frac {\partial u}{\partial z} +\varepsilon ^2\frac {\partial w}{\partial z} \right )\right ) \left ( \frac {\partial u}{\partial z} +\varepsilon ^2\frac {\partial w}{\partial z} \right ). \end{align}

This equation is averaged over the depth using the boundary conditions and neglecting all terms smaller than $O(\varepsilon \delta ^{1/n-1})$ . We omit this lengthy but algebraically straightforward derivation. The conservative left-hand side matches the classical structure formulated in Richard & Gavrilyuk (Reference Richard and Gavrilyuk2012), whereas the right-hand side is given by relaxation terms similar to those introduced for the momentum (4.18). The final depth-averaged energy equation reads

(4.21) \begin{align} \frac {\partial }{\partial t} \left ( \frac {h U^2 }{2} +\frac {h^3 \varphi }{2} +\frac {h^2}{2Fr^2} \right ) +\frac {\partial }{\partial x} \left (\frac {hU^3}{2} +\frac {3h^3U\varphi }{2} +\frac {h^2U}{Fr^2} \right ) \nonumber \\ =\frac {U}{\varepsilon \textit{Re}}\left [\lambda h (1-q(n,\xi ) ) - \textit{Bi} - \alpha _1(n) \frac {U}{h}\left |\frac {U}{h}\right |^{n-1} \right ]\varGamma _1(n,\xi ,\zeta ) \nonumber \\ +\frac {U\varphi ^{\frac {n-2}{2}}}{\varepsilon \textit{Re}}\left [\varphi - \alpha _2(n,\xi )\frac {U^2}{h^2} \right ]\varGamma _2(n,\xi ,\zeta ). \end{align}

Here the functions $\varGamma _1$ and $\varGamma _2$ are given in Appendix B. The first relaxation term on the right-hand side of (4.21) interprets as the balance between the power of weight, of yield stress and of viscous friction. The second term corresponds to the relaxation for the enstrophy $\varphi$ .

By combining (4.3), (4.18) and (4.21), an evolution equation for the enstrophy can be derived. This equation can be written as

(4.22) \begin{align} \frac {h^2}{2} \left ( \frac {\partial h\varphi }{\partial t}+\frac {\partial hU\varphi }{\partial x}\right )= \frac {U}{\varepsilon \textit{Re}}\left [\lambda h (1-q(n,\xi ) ) - \textit{Bi} - \alpha _1(n)\frac {U}{h}\left |\frac {U}{h}\right |^{n-1} \right ]r_1(n,\xi ,\zeta ) \nonumber \\ +\frac {U\varphi ^{\frac {n-2}{2}}}{\varepsilon \textit{Re}}\left [\varphi - \alpha _2(n,\xi )\frac {U^2}{h^2} \right ]r_2(n,\xi ,\zeta ). \end{align}

Here the functions $r_1$ and $r_2$ are given in Appendix B. It can be checked that $r_2(n,\xi ,\zeta )\lt 0$ , indicating that the enstrophy relaxes properly towards its equilibrium value $\alpha _2(\xi ,n) U^2/h^2$ .

The system of conservation (4.3), (4.18), (4.21) and the system (4.3), (4.18), (4.22) are hyperbolic and create discontinuities – or shocks – in finite time. Apart from shocks, the two systems are equivalent. In the case of a shock, solving system (4.3), (4.18), (4.21) leads to conservation of energy and production of enstrophy. On the other hand, solving system (4.3), (4.18), (4.22) leads to conservation of enstrophy and dissipation of energy. A production of enstrophy at shocks is relevant in open-channel hydraulics to model the apparition of a turbulent roller in hydraulic jumps or turbulent roll waves with a mechanism similar to that of breaking waves (Richard & Gavrilyuk Reference Richard and Gavrilyuk2012, Reference Richard and Gavrilyuk2013; Richard Reference Richard2024). Conversely, the enstrophy is conserved – and the energy is dissipated – when the enstrophy is due to bottom friction or viscous shearing (see Richard (Reference Richard2024) and Deleage & Richard (Reference Denisenko, Richard and Chambon2025), for turbulent roll waves and granular flows, respectively). In the present case of a viscoplastic flow, enstrophy arises from viscous shearing associated with bottom friction, in the absence of turbulent rollers or wave-breaking mechanisms. Therefore, the enstrophy (4.22) is solved instead of the energy (4.21), along with the mass and momentum (4.3) and (4.18).

4.4. Structure of the model for material stoppage

An important requirement for viscoplastic flow models is the ability to describe the stoppage of the material. However, the structure of the system (4.3), (4.18) and (4.22) cannot account for this feature. Formally, the final shape of the material after stoppage should satisfy the condition $h\leqslant h_p$ for any $x$ , since the opposite situation would immediately generate yielding due to the gravity term $\lambda h$ being larger than the plastic term $ \textit{Bi}$ . However, the asymptotic expansions derived in § 3 are valid only for $h \gt h_p$ . Hence, it cannot be guaranteed that the coefficients $\beta _1, \beta _2$ , $r_1$ and $r_2$ remain well-behaved for $h\leqslant h_p$ . This drawback comes from the fact that the thickness of the pseudoplug is supposed to be constant in the asymptotic expansions. As we show below, the coefficients can nevertheless be modified to properly capture stoppage, although this regime is formally outside of the validity of the asymptotic approach. To that aim we use the following strategy: let us first consider that $h\leqslant h_p$ and that the material is highly sheared. In this case, the material is expected to behave similarly to a power-law fluid, so that instead of the coefficients $\beta _1, \beta _2$ , $r_1$ and $r_2$ , we expect to recover their power-law limits:

(4.23) \begin{align} \beta _1(0,0,n) = 1, \\[-28pt]\nonumber \end{align}

(4.24) \begin{align} r_1(0,0,n) = 0, \\[-28pt]\nonumber \end{align}

(4.25) \begin{align} \beta _2^*(n) = \beta _2(0,0,n) \equiv \frac {1}{2} n^{1-3 n/2} (2 n+1)^{n-1} (3 n+2)^{n/2} (4 n+3), \\[-28pt]\nonumber \end{align}

(4.26) \begin{align} r_2^*(n) = r_2(0,0,n) \equiv -\frac {1}{2} n^{1-3 n/2} (n+1) (2 n+1)^{n-2} (3 n+2)^{n/2} (4 n+3).\\[0pt]\nonumber \end{align}

To ensure a consistent transition from the original coefficients $\beta _1, \beta _2$ , $r_1$ and $r_2$ to their power-law limits (4.23)–(4.26) for $h\leqslant h_p$ , let us introduce the following function:

(4.27) \begin{align} f(n,\varphi ,\xi ) = \begin{cases} \exp \left (2-\frac {\displaystyle \varphi }{\displaystyle \varphi ^{(0)}(n,\xi )}-\frac {\displaystyle \varphi ^{(0)}(n,\xi )}{\displaystyle \varphi }\right ), \quad \text{for} \ \ \xi \lt 1 \ \ \text{and} \ \ \varphi \gt 0, \\ 0, \quad \text{otherwise}, \end{cases} \end{align}

where $\varphi ^{(0)}$ is given by (4.9). An important property of $f$ is that a multiplication of the source terms by this function does not change the order of consistency of the model, since $f(n,\varphi ,\xi ) = 1 + O(\varepsilon ^2)$ . This function can thus be used to cancel out unwanted terms when $h \leqslant h_p$ . Taking into account this property, the momentum (4.18) and enstrophy balance (4.22) can be rewritten consistently as

(4.28) \begin{align} &\frac {\partial }{\partial t} \left ( hU \right ) +\frac {\partial }{\partial x} \left (hU^2+h^3\varphi + \frac {h^2}{2Fr^2} \right ) = \nonumber \\& \quad + \!\frac {1}{\varepsilon \textit{Re}}\!\left [\!\lambda h (1-q(n,\xi ) ) - \textit{Bi} - \alpha _1(n)\frac { U}{h}\left |\frac { U}{h}\right |^{n-1} \right ] \Big (1+\big [\beta _1(n,\xi ,\zeta )-1\big ] f(n,\varphi ,\xi ) \Big ) \nonumber \\& \quad +\frac {\varphi ^{\frac {n-2}{2}}}{\varepsilon \textit{Re}}\left [\varphi - \alpha _1(n)\frac {U^2}{h^2} \right ]\Big (\beta _2^*(n)+\big [\beta _2(n,\xi ,\zeta )-\beta _2^*(n)\big ] f(n,\varphi ,\xi ) \Big ), \end{align}

(4.29) \begin{align} &\frac {h^2}{2} \left ( \frac {\partial h\varphi }{\partial t}+\frac {\partial hU\varphi }{\partial x}\right )= \nonumber \\& \quad \frac {U}{\varepsilon \textit{Re}}\left [\lambda h (1-q(n,\xi ) ) - \textit{Bi} - \alpha _1(n)\frac {U}{h}\left |\frac {U}{h}\right |^{n-1} \right ]r_1(\xi ,n,\zeta ) f(n,\varphi ,\xi ) \nonumber \\& \quad +\frac {U\varphi ^{\frac {n-2}{2}}}{\varepsilon \textit{Re}}\left [\varphi - \alpha _1(n,\xi )\frac {U^2}{h^2} \right ]\Big (r_2^*(n)+\big [r_2(n,\xi ,\zeta )-r_2^*(n)\big ] f(n,\varphi ,\xi ) \Big ). \end{align}

It should be noted that the function $q(n,\xi )$ is also not well-defined for $\xi \lt 1$ ( $h\lt h_p$ ). Hence, we consider its limiting value $q(n,1)=0$ for any $\xi \leqslant 1$ .

If now $h\leqslant h_p$ and the material is only slightly sheared, the fluid is expected to slow down and progressively stop. In this case, the fluid velocity and enstrophy are small, such that the dominant terms in (4.28) and (4.29) are the hydrostatic pressure, the gravity and the yield stress. In particular, when the material is stopped ( $U=0$ , $\varphi =0$ and $h\leqslant h_p$ ), the momentum balance (4.28) properly reduces to

(4.30) \begin{align} \lambda h\left (1 -\varepsilon \cot {\theta } \frac {\partial h}{\partial x}\right ) = \textit{Bi}, \end{align}

which is the classical equation for a static surge stopped on an inclined plane, as used in several previous studies (Balmforth et al. Reference Balmforth, Craster, Rust and Sassi2006; Hogg & Matson Reference Hogg and Matson2009; Dubash et al. Reference Dubash, Balmforth, Slim and Cochard2009). Note also that both the left- and right-hand sides of the enstrophy (4.29) are zero when the fluid is stopped. This is not the case for the original (4.22): here the first source term is never zero (except if $h=h_p$ ), which would prevent equality between the left- and right-hand sides at stoppage. Hence, the consistent modification of the source term coefficients in (4.28) and (4.29) properly captures the expected flow behaviour for $h\leqslant h_p$ both when the fluid is highly and slightly sheared: in the former case the fluid behaves as a power-law fluid, while in the latter case the material converges to a stoppage.

Lastly, note that (4.30) only refers to the relation between driving and resisting forces after stoppage. In general, the full stoppage criterion reads

(4.31) \begin{align} \left | \lambda h\left (1 -\varepsilon \cot {\theta } \frac {\partial h}{\partial x}\right ) -\varepsilon \textit{Re}\left [\frac {\partial }{\partial t}\left ( hU \right ) +\frac {\partial }{\partial x} \big (hU^2+h^3\varphi \big )\right ] \right | \leqslant \textit{Bi} \ \ \Rightarrow \ \ (U,\varphi )=0, \end{align}

meaning that the material moves only when the sum of the driving forces (weight, pressure gradient and inertia) overcomes the yield stress. Furthermore, (4.30) describes only one of the possible equilibrium shapes of the viscoplastic material at stoppage. This scenario happens when the fluid velocity decreases to zero as the driving forces approach the balance with the yield stress. In a more general case, it is also possible to have situations for which the fluid velocity is still large, while the sum of the gravity and hydrostatic pressure has already become smaller than the yield stress. In this case, according to the momentum (4.28), the source terms will become largely negative, slowing down the fluid and leading to an abrupt stoppage. Therefore, a more general condition at stoppage is

(4.32) \begin{align} \lambda h\left (1 -\varepsilon \cot {\theta } \frac {\partial h}{\partial x}\right ) \leqslant \textit{Bi}. \end{align}

This issue is further elaborated in § 5.3.2, where we simulate dam-break surges and study different shapes after stoppage.

4.5. The model in dimensional variables

In dimensional variables, the system of (4.3), (4.28) and (4.29) is rewritten as

(4.33) \begin{align}&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\frac {\partial h}{\partial t} +\frac {\partial hU}{\partial x}=0,\end{align}

(4.34) \begin{align} &\frac {\partial }{\partial t} \left ({h}{U}\right ) +\frac {\partial }{\partial {{x}}} \big ({h}{U}^2+ h^3\varphi \big ) = -\frac {\partial }{\partial x}\frac {g h^2\cos {\theta }}{2} \nonumber \\&\quad + \!\left [\! g h \sin {\theta } ( 1- q(n,\xi ) )-\frac {\tau _c}{\rho } -\alpha _1(n)\nu \frac {U}{h}\left |\frac {U}{h}\right |^{n-1} \right ]\!\Big (1\!+\!\big [\beta _1(n,\xi ,\zeta )-1\big ] f(n,\varphi ,\xi ) \Big ) \nonumber \\&\quad +\varphi ^{\frac {n-2}{2}} \Bigg [ \varphi -\alpha _2(n,\xi )\frac { U^2}{h^2} \Bigg ] \Big (\beta _2^*(n)+\big [\beta _2(n,\xi ,\zeta )-\beta _2^*(n)\big ] f(n,\varphi ,\xi ) \Big ), \end{align}

(4.35) \begin{align} &\frac {h^2}{2} \left ( \frac {\partial h\varphi }{\partial t}+\frac {\partial hU\varphi }{\partial x}\right )= \nonumber \\& \quad U\left [ g h \sin {\theta } ( 1- q(n,\xi ) )-\frac {\tau _c}{\rho } -\alpha _1(n)\nu \frac {U}{h}\left |\frac {U}{h}\right |^{n-1} \right ] r_1(n,\xi ,\zeta ) f(n,\varphi ,\xi ) \nonumber \\& \quad + U\varphi ^{\frac {n-2}{2}} \Bigg [ \varphi -\alpha _2(n,\xi )\frac { U^2}{h^2} \Bigg ] \Big (r_2^*(n)+\big [r_2(n,\xi ,\zeta )-r_2^*(n)\big ] f(n,\varphi ,\xi ) \Big ), \end{align}

where $\nu =K / \rho$ , $\xi =\tau _c/(\rho g h \sin {\theta })$ and $\zeta =\delta /h$ . The functions $\alpha _i$ ( $i=1,2$ ) and $q$ are defined in § 4.2. The functions $r_i$ and $\beta _i$ ( $i=1,2)$ are defined in Appendix B. The functions $\beta ^*_2$ , $r^*_2$ and $f$ are defined in § 4.4. Note that the parameter of pseudoplug shearing $\delta$ is here dimensional and can be interpreted as the thickness of the thin intermediate layer located between the sheared layer and the pseudoplug, in which the ‘Bingham plateau’ regularization contributes to the first-order correction of the velocity asymptotic expansion (A11). For completeness, the stoppage criterion (4.31) in dimensional variables reads

(4.36) \begin{align} \left | gh\sin {\theta }-\frac {\partial }{\partial x}\frac {gh^2\cos \theta }{2}- \frac {\partial }{\partial t} \left ( hU \right ) -\frac {\partial }{\partial x} \big (hU^2+h^3\varphi \big ) \right | \leqslant \frac {\tau _c}{\rho } \ \ \Rightarrow \ \ (U,\varphi ) = 0. \end{align}

The left-hand side of the system (4.33)–(4.35) has the structure of the Euler equations for compressible fluids with $h^3\varphi$ and $\varphi$ playing the roles of pressure and entropy, respectively. The model is therefore fully hyperbolic, enabling the implementation of an efficient well-balanced numerical scheme, described in detail in Appendix C. The well-balanced property, i.e. the ability of the scheme to capture non-trivial static solutions satisfying the inequality (4.32), is ensured by treating the hydrostatic pressure as a source term and by employing a discretized stoppage criterion (C11). Importantly, although the final model is highly nonlinear and complex, it can still be regarded as computationally efficient. The first step of the numerical scheme, which involves solving the homogeneous system, is handled using the standard Rusanov Riemann solvers developed for classical systems of gas dynamics (e.g. see Toro Reference Toro1999). This step is therefore straightforward and computationally inexpensive. The second step, which involves integrating the source terms, may appear more computationally demanding due to the length and complexity of the source-term expressions. Nevertheless, the use of a semi-implicit scheme allows the resulting algebraic system to be solved analytically with respect to the flow variables, thereby keeping the computational cost relatively low.

We recall that the model (4.33)–(4.35) is applicable only for $n\lt 1$ due to the specific regularization of the viscous stress tensor at $O(\varepsilon )$ , except in the case $\delta =0$ , which corresponds to the model with zero shearing in the pseudoplug and is valid for any value of $n$ . Although it is possible to adapt the regularization for $n\gt 1$ , the existing real viscoplastic fluids are generally characterized by $n\lt 1$ , and therefore we restricted our analysis to this case.

For the particular case of Bingham fluids ( $n=1)$ , the model is presented in Appendix D. This model is derived similarly to the model described above for Herschel–Bulkley fluids, but using the analytical asymptotic expansions obtained in Denisenko et al. (Reference Denisenko, Richard and Chambon2023) (without the regularization parameter $\delta$ ). Note that the model for Bingham fluids (D1)–(D3) cannot be obtained from the model (4.33)–(4.35) as $n\to 1$ . This issue is due to the divergence of the effective viscosity for $n\lt 1$ , that immediately leads to zero shearing in the pseudoplug for the asymptotic expansion at $O(\varepsilon )$ . The introduction of the ‘Bingham plateau’ regularization with the parameter $\delta$ , which is introduced to recover non-zero shearing in the pseudoplug for $n\lt 1$ , also slightly alters the expansion in the sheared layer ( $z\lt h-h_p$ ) at $O(\varepsilon )$ for any values of $n$ . Hence, at $n=1$ , the first-order correction (A11)–(A17) differs from that in Denisenko et al. (Reference Denisenko, Richard and Chambon2023) for any values of $\delta$ .

Note also that the model for Bingham fluids presented in Appendix D differs from the model previously derived in Denisenko et al. (Reference Denisenko, Richard and Chambon2023), since the latter was not adapted to describe flow stoppage. In particular, this new model includes the following modifications: the new definition of the enstrophy (4.7), which eliminates the need for the enstrophy to relax towards zero at $h=h_p$ ; the inclusion of the extra yield-stress friction term $q$ , allowing the removal of the dependence of the viscous friction on the yield stress and avoiding the divergence at $h=h_p$ ; the adaptation of the coefficients in front of the source terms to ensure proper behaviour for $h \leqslant h_p$ and to correctly capture stoppage; the use of the enstrophy balance equation instead of the energy balance equation in order to conserve the enstrophy variable associated with the viscous shearing (Richard Reference Richard2024).

For the case of power-law fluids (i.e. without yield stress), the model derived as a particular case of (4.33)–(4.35) with $\tau _c=0$ and $\delta =0$ is presented in Appendix E. Note that this model is applicable for any value of $n$ , since no regularization is required in this case. To the best of our knowledge, this is the first time that the three-equation approach has been applied to the power-law rheology.

5.1. Long-wave instability

Steady and uniform equilibrium flows down an inclined plane are unstable in a specific range of the flow parameters. Balmforth & Liu (Reference Balmforth and Liu2004) studied the linear stability of the equilibrium flow based on the linearized Cauchy equations (or generalized Orr–Sommerfeld equations) in the case of a Herschel–Bulkley fluid with zero shearing in the pseudoplug. In order to compare our results with this study, we use here the same characteristic velocity, namely $u_0=(\rho g {h}_0^{1+n} \sin {\theta } / K)^{1/n}$ . Note that this choice leads to $\lambda =1$ .

The system of (4.3), (4.28) and (4.29) is linearized around the leading-order solution (4.2) and (4.9). Let us consider $h=1+h^\prime$ , $U= U^{(0)} + U^\prime$ and $\varphi = \varphi ^{(0)} + \varphi ^\prime$ with small sinusoidal perturbations: $[h^\prime,U^\prime,\varphi ^\prime] = [A_1, A_2, A_3] \exp [\text{i}k(x-ct)]$ . The parameter $k$ is the wavenumber and $c$ is the phase velocity. The dispersion relation is derived by equating the determinant of the linearized system to zero, which yields

(5.1) \begin{align} c=\lambda ^{\frac {1}{n}}(1-h_p)^{\frac {1}{n}} + \text{i}\varepsilon k \, c_1(\textit{Re},n,h_p,\delta )+O(\varepsilon ^2). \end{align}

The expression of the first-order correction to the phase velocity $c_1$ is lengthy and does not need to be provided here. In the long-wave limit (large $k$ ), the base flow is stable if $c_1\lt 0$ . This leads to stability criterion $ \textit{Re}\lt \textit{Re}_c$ , with the critical Reynolds number $ \textit{Re}_c$ given by

(5.2) \begin{align} \textit{Re}_c = \frac {\cot {\theta } (n+1) (n+2) (3 n+1) (3 n+2) (4 n+1) }{\left (1-\textit{Bi}\right ){}^{\frac {2}{n}} D(n,\textit{Bi},\delta ) } \Big [\delta ^{1/n-1} (n+1) (2 n+1)\nonumber \\[5pt]+3 n \big (2 n \textit{Bi} (n \textit{Bi}+1)+n+1\big ) (1-\textit{Bi}){}^{\frac {1}{n}}\Big ], \end{align}

where

(5.3) \begin{align}& D(n,\textit{Bi},\delta )= 3 n (3 n+1) (4 n+1) \left (1-\textit{Bi}\right ){}^{\frac {1}{n}} \big (n \textit{Bi} \big (n (n (6 n+19)+11) \textit{Bi} \left (n \textit{Bi}+2\right )\nonumber \\[4pt]& \quad -\delta ^{1/n-1} (n+1) (n+2) (3 n+2) \big (n \textit{Bi} \big (3 n \textit{Bi} \big(4 n^2 \big(\textit{Bi}^2-3 \big)+n (4 \textit{Bi}-5)+1\big )\nonumber \\[4pt]& \quad -n (32 n\!+\!15)-1\!\big )\! +\!(n\!+\!2) (9 n\!+\!7)\big )\!+\!2 (n\!+\!2) (n\!+\!1)^2\big ) -(n\!+\!1) (3n\!+\!1) (5 n\!+\!1)\big ). \end{align}

Due to the term $(1-\textit{Bi})^{({2}/{n})}$ in the denominator of (5.2), the critical Reynolds number $ \textit{Re}_c$ increases as the Bingham number increases. This corresponds to a stabilizing effect of the plasticity, which was also obtained in former studies (Balmforth & Liu Reference Balmforth and Liu2004; Denisenko et al. Reference Denisenko, Richard and Chambon2023). Note that a different result would be obtained by implementing the non-zero shearing in the pseudoplug through asymptotic expansions built in the framework of the classical tensorial extension of the constitutive law. In this case, a destabilizing effect of plasticity is observed at large slope angles (see Denisenko et al. (Reference Denisenko, Richard and Chambon2023) for Bingham fluids), leading to an implausible prediction of absolutely unstable rigid solid under small perturbations. Although there are no experimental reports of instability of viscoplastic fluids at large slope angles, the existing results for smaller angles consistently show a stabilizing effect of plasticity (e.g. Noma et al. Reference Mounkaila Noma, Dagois-Bohy, Millet, Botton, Henry and Ben Hadid2021).

In the case $\delta = 0$ , i.e. for zero shearing in the pseudoplug, (5.2) reduces to

(5.4) \begin{align} \textit{Re}_c=\frac {\cot {\theta }(n+1) (n+2) (3 n+2) \left (1-\textit{Bi}\right ){}^{-2/n} \left (2 n \textit{Bi} \left (n \textit{Bi}+1\right )+n+1\right )}{n \textit{Bi} \left (n (n (6 n+19)+11) \textit{Bi} \left (n \textit{Bi}+2\right )+(n+2) (9 n+7)\right )+2 (n+2) (n+1)^2}, \end{align}

which exactly coincides with the expression obtained by Balmforth & Liu (Reference Balmforth and Liu2004). This result is made possible by the consistency of our three-equation model and the fact that for $\delta =0$ the velocity expansion (A11)–(A17) is the same as that used by Balmforth & Liu (Reference Balmforth and Liu2004).

Figure 2.

Instability thresholds predicted by the three-equation model: critical Reynolds number $ \textit{Re}_c \tan \theta$ as a function of the Bingham number $ \textit{Bi}$ . For $n=0.4$ , the red and blue curves correspond to the threshold (5.2) with $\delta =0$ and $\delta =0.2$ , respectively. For $n=1$ , the red curve corresponds to (5.2) with $\delta =0$ , while the blue curve corresponds to the threshold (D9) obtained for the Bingham fluid model (D1)–(D3).

Figure 2 presents the evolution of the critical Reynolds number $ \textit{Re}_c$ as a function of $ \textit{Bi}$ for different values of $n$ and $\delta$ . For power-law fluid ( $ \textit{Bi}=0)$ , $ \textit{Re}_c$ decreases as $n$ decreases. An opposite trend is observed for viscoplastic fluids if the Bingham number is large enough: in this case a decrease in the power-law index $n$ leads to stabilization of the flow. This effect was already highlighted in Balmforth & Liu (Reference Balmforth and Liu2004). Conversely, variations in $\delta$ produce hardly visible differences: the curves corresponding to $\delta =0$ and to $\delta =0.2$ are almost identical. Hence, the influence of a slight shearing of the pseudoplug on the stability criterion is small. In detail, the curves for $\delta =0.2$ are slightly below those for $\delta =0$ , showing a slight destabilizing effect of the shearing in the pseudoplug. This finding supports the intuitive expectation that shearing in the pseudoplug, by suppressing fluid rigidity, makes the flow more prone to instability.

Based on comparisons with experimental data, Noma et al. (Reference Mounkaila Noma, Dagois-Bohy, Millet, Botton, Henry and Ben Hadid2021) showed that (5.4) predicts relatively well the instability thresholds for the real viscoplastic fluids. As shown in figure 3, our criterion (5.2) yields the same agreement with experiments, irrespective of whether $\delta$ is zero or non-zero.

Figure 3.

Instability thresholds obtained from our model with $\delta =0$ (red curve) and with $\delta =0.2$ (blue curve) compared with the experimental results of Noma et al. (Reference Mounkaila Noma, Dagois-Bohy, Millet, Botton, Henry and Ben Hadid2021) (grey dots) for Carbopol (a) and kaolin (b).

5.2. Simulation of roll waves

We study here the development of roll waves from equilibrium flows under small perturbations. Such waves correspond to flow regimes that remain relatively close to equilibrium and therefore lie within the validity domain of the assumptions underlying the model derivation. An initially uniform flow satisfying the leading-order solution (4.2) and (4.9) is perturbed by applying a small sinusoidal disturbance of fixed frequency at the entrance of the system. The length of the computational domain is $L=3.0$ m, the mesh size is ${\rm d}x=0.5\,$ mm and Courant–Friedrichs–Lewy number is $0.5$ . The depth of the uniform flow is $h_0=10\,$ mm and the angle of the slope is $\theta =20^{\circ }$ . At the entrance, the flow height $h$ is perturbed at a frequency of $2.0$ Hz, the discharge $hU$ is preserved and the enstrophy $\varphi$ is taken equal to $\varphi ^{(0)}(h)$ with the perturbed value of $h$ . At the outlet, a Neumann boundary condition is considered. We explore the influence of the shear-thinning and the plasticity on the simulated roll waves by varying the values of $n$ and $ \textit{Bi}$ . Additionally, we explore the influence of the shearing in the pseudoplug by varying the parameter $\delta$ . The values of Reynolds and Froude numbers are fixed for all simulations ( $ \textit{Re}=37.28$ and $ \textit{Fr}=3.68$ ), in order to isolate the influence of variations in the other parameters.

Figure 4.

Numerical simulations of roll waves with the three-equation model: evolution of flow depth $h$ as a function of distance from inlet $x$ for different values of $n$ and $ \textit{Bi}$ . For $n=0.4$ , the red and blue curves are obtained from (4.33)–(4.35) with $\delta =0$ and $\delta =0.2$ , respectively. For $n=1$ , the red curve is computed from (4.33)–(4.35) with $\delta =0$ , while the blue curve is computed from the Bingham fluid model (D1)–(D3). Other parameters are $ \textit{Re}=37.28$ , $ \textit{Fr}=3.68$ , $\lambda =1$ and $\theta = 20^{\circ }$ .

Figure 5.

Comparison between simulated roll waves and the experimental results of Fiorot et al. (Reference Fiorot, da Rocho, Möller, Pereira, da Cunha and de Maciel2024): variation of the depth $h$ at a function of time at $x=1.5$ from the disturbance – grey dots corresponds to the experimental data Fiorot et al. (Reference Fiorot, da Rocho, Möller, Pereira, da Cunha and de Maciel2024), blue and red curves correspond to numerical results with $\delta =0$ and $\delta =0.2$ , respectively. Other parameters are indicated in the text.

Figure 4 shows the simulated roll waves for all parameter cases in dimensional form. It is clearly observed that an increase of plasticity, i.e. of the Bingham number, leads to a decrease of the roll wave amplitude and wavelength. This result has already been highlighted by Denisenko et al. (Reference Denisenko, Richard and Chambon2023) and can be related to the stabilizing effect of plasticity. The influence of shear-thinning, i.e. of the power-law index $n$ , is less straightforward. For $ \textit{Bi}=0.35$ , a decrease in $n$ leads to a decrease of the roll wave amplitude and their wavelength. For $ \textit{Bi}=0.05$ , a decrease in $n$ leads to an increase of the roll wave amplitude and to a decrease of their wavelength. These different effects can be related to the evolution of the instability criterion and the crossover observed for $ \textit{Bi} \approx 0.1$ in figure 2: in general, the larger the distance $ \textit{Re} - \textit{Re}_c$ between the fixed Reynolds number and the critical value, the larger the wave amplitude. Generally, however, this rule seems to be valid only in the vicinity of the critical values of the Reynolds number.

Figure 4 also illustrates the influence of shearing in the pseudoplug. Recall that, for $n=1$ (Bingham fluid), the model with non-zero shearing in the pseudoplug is derived without any parameter $\delta$ (see Appendix D), while, for $n \neq 1$ , the model with zero shearing in the pseudoplug is derived from the general model (4.33)–(4.35) with $\delta =0$ . Although the difference $ \textit{Re}-\textit{Re}_c$ slightly increases as $\delta$ increases (see figure 2), an increase in $\delta$ generally leads to a slight decrease in the roll wave amplitudes, for both values of $n$ . This behaviour is presumably due to the fact that shearing in the pseudoplug layer introduces an extra energy dissipation, which slightly damps the growth of the roll waves.

Fiorot et al. (Reference Fiorot, da Rocho, Möller, Pereira, da Cunha and de Maciel2024) recently reported on experimental roll waves generated with a viscoplastic fluid. Figure 5 shows a comparison between these experimental results and the prediction of the three-equation model. For this simulation, we use the same parameters as in the experiments: $n=0.39$ , $\theta =20^{\circ }$ , $\rho =1003.2$ kg m $^{-3}$ , $\tau _c=9.56$ Pa, $K=5.73$ Pa s $^n$ and $h_0=12.56$ mm. The length of the channel is 2.5 m, and the frequency applied at the inlet is $1.5$ Hz.

It is observed that the agreement with the experiments is relatively good behind the wavefronts. However, the amplitude of the waves in the simulation is notably larger than in the experiments. This discrepancy is presumably related to the fact that the model at $O(\varepsilon )$ does not include any diffusion effect, such that the simulated roll waves are in fact shock waves. Diffusion effects, proportional to second-order derivatives, would smooth out these shocks and lead to a decrease in wave amplitude. In practice, considering such effects in a consistent manner would require the derivation of higher-order terms by continuing the asymptotic expansions up to $O(\varepsilon ^2)$ . Note also that, in this case, the shearing in the pseudoplug zone (parameter $\delta$ ) has essentially no effect on the free-surface shape.

5.3. Dam-break simulations

We consider here the classical dam-break problem for a viscoplastic material. In this case, due to the existence of a dry front and to the possibility of a flow stoppage, the flow is far from the equilibrium flow and therefore lies outside the assumptions underlying the model derivation. This case can thus be used to assess the validity of the strategy presented in § 4.4.

An initial volume of material is filled into a rectangular box of length $l$ , which is located upstream of the inclined plane. The rear end of this reservoir is chosen as the origin of the coordinate system. The box is equipped with a gate perpendicular to the slope, which is suddenly opened at $t = 0$ . The initial flow depth is thus given by

(5.5) \begin{align} h_0(x)= \begin{cases} h_g+(x-l)\tan {\theta }, \quad &0\leqslant x \leqslant l, \\ 0, &l\leqslant x, \end{cases} \end{align}

with $h_g$ the fluid height at the gate.

Figure 6.

Simulation of a viscoplastic dam-break: front position $x_f$ as a function of time $(a)$ and flow depth profiles at different times $(b)$ . Red curves, numerical scheme with hydrostatic pressure taken as a source term; blue curves, numerical scheme with hydrostatic pressure taken as the part of the momentum flux. Parameters of the flow are $\tau _c=350\,$ Pa, $K=200\,$ Pa s $^n$ , $n=0.4$ , $\rho =1000\,$ kg m $^{-3}$ , $\theta =20^{\circ }$ and $\delta =0$ .

The typical dynamics of the flow is presented in figure 6. Two different numerical solutions of the three-equation model are compared, which differ by treating the hydrostatic pressure as a source term (red curve) or as a part of the momentum flux (blue line) (see Appendix C.5). Two main flow regimes separated by a relatively abrupt transition are observed. At early times, an inertially dominated regime lasts for approximately 1 s and is characterized by a rapid propagation of the front. The front then slows down in the second regime, which corresponds to a balance between gravity, yield stress, viscous friction and hydrostatic pressure. The two numerical schemes provide similar results for the first regime, both in terms of front velocity and free-surface shape. In particular, they reproduce a contact angle close to $90^{\circ }$ at the front, as expected for such viscous fluids. In the second regime, however, the non-well-balanced numerical scheme with the pressure in the flux tends to overestimate the front speed and loses the 90 $^{\circ }$ contact angle. This ‘smearing’ of the front is due to the fact that, with this scheme, the numerical diffusion becomes dominant in the mass equation at small velocities (for more details see Appendix C.5). In contrast, the well-balanced numerical scheme with the pressure in the source terms preserves a realistic front shape as the flow progressively slows down.

Similar differences between well-balanced and not well-balanced schemes were already highlighted by Fernández-Nieto et al. (Reference Fernández-Nieto, Garres-Díaz and Vigneaux2023), who used another well-balanced scheme based on hydrostatic reconstruction of the depth variable. However, the authors reported that their scheme becomes mesh-dependent in the slowly creeping regime, leading to an overestimation of the front velocity. In contrast, the scheme proposed in the present work does not show such a drawback.

5.3.1. Comparison with experiments

To evaluate the predictions of the model, we simulate the experiments reported by Ancey & Cochard (Reference Ancey and Cochard2009). These authors used a Carbopol material characterized by the following parameters: $\tau _c=32.1\,$ Pa, $K=78\,$ Pa s $^n$ , $n=0.388$ and $\rho =1000\,$ kg m $^3$ . The length of the reservoir is $l=0.51\,$ m and four values of slope angle are considered: $\theta =6^{\circ };\,12^{\circ };\,18^{\circ };\,24^{\circ }$ . For these slopes, the gate aperture $h_g$ is 0.26 m, 0.23 m, 0.20 m and 0.177 m, respectively. For the numerical simulations, we use the consistent three-equation model (4.33)–(4.35) with $\delta =0$ . We checked that varying the parameter $\delta$ (shearing in the pseudoplug) has essentially no influence in this case. To illustrate the influence of inertia, we also plot the results of the inertialess lubrication model presented by Ancey & Cochard (Reference Ancey and Cochard2009).

Figure 7 compares the evolution of front position for the two models and the experiments. The experimental results exhibit two main regimes separated by a relatively abrupt transition. At early times, an inertially dominated regime persists for approximately one second and is characterized by rapid front propagation. Next, the flow enters a second, slow creeping regime characterized by a deceleration of the front. It is observed that the lubrication model, which does not include inertial terms, cannot capture the first regime. In the second regime, the lubrication model predicts a trend parallel to the experimental data but with a slight shift. In contrast, the three-equation model captures the first regime and the transition relatively well. In the second regime, the three-equation model coincides almost exactly with the lubrication model, as expected, and the slight shift with the experimental data is again observed. This shift may relate to the sidewall friction, as suggested by Andreini, Epely-Chauvin & Ancey (Reference Andreini, Epely-Chauvin and Ancey2012). These authors pointed out that similar shifts appear in the modelling of dam-break problems for Newtonian fluids and can be eliminated by incorporating a sidewall friction term into the model. However, to the best of our knowledge, there is no general method for modelling sidewall friction in open-channel flows of viscoplastic fluids (e.g. Coussot Reference Coussot1994).

Figure 7.

Comparison between dam-break simulations and experiments: front position $x_f$ as a function of time for different slope angles. Experimental data and results of the lubrication model come from Ancey & Cochard (Reference Ancey and Cochard2009). The three-equation model (4.33)–(4.35) is simulated with $\delta =0.0$ .

Another discrepancy with the experiments is observed at long times for the smaller angles ( $\theta =12^{\circ }$ and $6^{\circ }$ ). In the predictions of both the lubrication and the three-equation model, the yield stress term becomes dominant, so that the material progressively approaches full stoppage. In the experiment, however, the material keeps moving without any noticeable slowdown. These deviations at long times could be due to some residual slip of the Carbopol on the bottom plane, such that the no-slip boundary condition (2.20) may not be completely valid in the experimental flows. A simple Navier-slip condition could potentially be incorporated into the model within the framework of the considered constitutive law, although this would lead to lengthy algebraic modifications in the asymptotic expansion.

Flow depth profiles predicted by the three-equation model are compared with experiments in figure 8. A satisfactory agreement is found despite the differences in front position. The fact that the maximum flow depth is not observed at the front in the experimental data at $t=1.0$ s (for $\theta =24^{\circ }$ and $\theta =18^{\circ }$ ) is likely a consequence of the finite opening time of the gate.

Figure 8.

Comparison between dam-break simulations and experiments: flow-depth profiles for different slope angles. The red curve corresponds to the three-equation model (4.33)–(4.35); the black dashed line corresponds to the experimental data of Ancey & Cochard (Reference Ancey and Cochard2009).

Overall, the new three-equation shallow-flow model appears to capture relatively well the main trends observed in the experiments. Only few previous studies attempted such comparisons between shallow-flow models and experiments for viscoplastic dam breaks. Ancey et al. (Reference Ancey, Andreini and Epely-Chauvin2012) considered an inconsistent shallow-flow model and captured the inertia-dominated regime. However, they got an even larger shift with the experiments in the second regime, probably because these experiments were performed in a narrower channel. Later, comparisons with the same experiments were carried out by Fernández-Nieto et al. (Reference Fernández-Nieto, Garres-Díaz and Vigneaux2023), who used a multilayer shallow-flow model. Their result slightly improves the agreement with the data, in particular due to using a well-balanced scheme based on hydrostatic reconstruction. However, in both these studies, the predictions of the shallow-flow models deviate from those of the lubrication model at large times. Such deviations can be attributed to the numerical diffusion of the employed schemes and are similar to the behaviour observed with our model when the hydrostatic pressure is kept in the flux (see figure 6). Recently, Muchiri et al. (Reference Muchiri, Hewett, Sellier, Moyers-Gonzalez and Monnier2024) proposed a new inconsistent shallow-flow model and reported good agreement with several experiments of Ancey & Cochard (Reference Ancey and Cochard2009). However, only large slope angles were considered. Additional comparisons at smaller slope angles would be useful. Furthermore, this model cannot describe stoppage of the material.

5.3.2. Flow stoppage

We focus here on the long-time dynamics and final shape of the dam-break flow. Again, let us recall that the essential ingredient to properly capture full stoppage in our simulations is to consider the hydrostatic pressure as a source term in the numerical scheme (see Appendix C.5 for more details). As discussed in § 4.4, the final shape of the deposits should satisfy inequality (4.32) on the balance between hydrostatic pressure, gravity force and yield stress. In case of equality, (4.32) yields the ordinary differential (4.30) that was considered in several previous studies (Balmforth et al. Reference Balmforth, Craster, Rust and Sassi2006; Hogg & Matson Reference Hogg and Matson2009) to describes static surges stopped on an inclined plane. The free-surface shape predicted by this equation was shown to be in good agreement with experimental deposits of kaolin fluids (Coussot Reference Coussot1994). However, as was shown by Balmforth et al. (Reference Balmforth, Craster, Rust and Sassi2006), deposit shapes do not systematically correspond to solutions of (4.30). Below we illustrate such phenomena in the framework of the three-equation model as well.

Figure 9.

Dam-break simulation and final stoppage. (a) Flow depth profiles at different times: red solid curves are the numerical solutions of model (4.33)–(4.35); red dotted curve is the final shape at stoppage; black dashed line is the solution of (4.30) for the given volume. (b) Front position $x_f$ as a function of time. (c) Velocity maximum $U_{\textit{max}}$ within the flow as a function of time. Parameters are $\tau _c=300$ Pa, $K=150$ Pa s $^n$ , $n=0.4$ , $\theta =20^{\circ }$ and $\rho =1000$ kg m $^{-3}$ .

Figure 10.

Dam-break simulation and final stoppage. Same legend as figure 9. Parameters are $\tau _c=1000$ Pa, $K=2$ Pa s $^n$ , $n=0.4$ , $\theta =20^{\circ }$ and $\rho =1000$ kg m $^{-3}$ .

Figures 9 and 10 present the results of two dam-break simulations characterized by different rheological parameters, including the final stoppage. In figure 9 $(a)$ , it is observed that the material slowly reaches a final shape satisfying (4.30). In the dynamics of the surge front (figure 9 b), the two regimes already described can be distinguished, namely the inertially dominated regime ( $0\lesssim t \lesssim 10^{0}$ s) and the lubrication regime ( $10^{0}\lesssim t \lesssim 10^{2}$ s). A plasticity-dominated regime ( $10^{2}\lesssim t \lesssim 10^{9}$ s) then takes over at long times and slowly leads to stoppage, which is not reached before $t \gtrsim 10^{8}$ s in this case. These different regimes are also visible in the evolution of the velocity maximum (figure 9 c). Note that at long times, the velocity maximum also demonstrates marked fluctuations. These fluctuations are associated with stop-and-go motions of the front. Typically, the front stops, while rear parts of the surge continue to move, causing a accumulation of material behind the front. This leads to a progressive increase of the pressure gradient until the yielding threshold is reached, resulting in a new displacement of the front.

In figure 10 the dynamics of the flow is markedly different. The material propagates quickly at the beginning but then suddenly stops. In this case, the shape after stoppage does not satisfy (4.30). This behaviour is due to the small value of consistency $K$ but relatively large value of yield stress $\tau _c$ . The reduction of $K$ decreases the viscous friction, so that the surges quickly reach a high velocity during the inertial regime. In turn, this leads to the sum of pressure gradient and weight becoming less than the yield stress within the flow, resulting in a negative source term in the momentum (4.34). Once inertia becomes negligible, the flow exhibits a rapid deceleration, ending with a non-trivial shape at stoppage. In this flow scenario, only two regimes are observed, namely the inertially dominated regime and the plasticity-dominated regime with a quick stoppage. The intermediate lubrication regime is not visible. This result illustrates the ability of our model to predict various stoppage shapes, which can be important for applications. Further studies could be conducted to systematically investigate the final deposit shapes as a function of the different flow parameters.