<figure><img src="image_1.png" /></figure>

HAL
open science

# A consistent three-equation shallow-flow model for Bingham fluids

Danila Denisenko, Gaël Loïc Richard, Guillaume Chambon

► To cite this version:

Danila Denisenko, Gaël Loïc Richard, Guillaume Chambon. A consistent three-equation shallow-flow model for Bingham fluids. Journal of Non-Newtonian Fluid Mechanics, 2023, 321, pp.105111. 10.1016/j.jnnfm.2023.105111. hal-04736101v3

HAL Id: hal-04736101

https://hal.inrae.fr/hal-04736101v3

Submitted on 2 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A consistent three-equation shallow-flow model for Bingham fluids

Danila Denisenko*, Gaël Loïc Richard and Guillaume Chambon

*Univ. Grenoble Alpes, INRAE, CNRS, IRD, Grenoble INP, IGE, Grenoble, 38000, France

## ARTICLE INFO

**Keywords:**

viscoplastic fluids
depth-averaged model
roll waves
instability
entrophy

## ABSTRACT

We derive a model for Bingham fluid flows down an inclined plane with a consistent asymptotic method in the shallow-flow approximation. The variables are expanded up to the first order of accuracy both in the sheared and pseudo-plug layers. The divergence of the strain rate, which is obtained in classical approaches, is here avoided by a specific regularization of the rheology allowing to implement a regular perturbation method in the whole fluid domain. Unlike classical regularization methods, the material is here characterized by a true yield stress. Below the yield point, the behavior is perfectly rigid. An alternative tensor expression of the constitutive law is proposed. In particular, the assumption of an alignment between the yield-stress tensor and the strain-rate tensor is removed. The model is derived by averaging the mass, momentum and energy balance equations over the depth. This yields a hyperbolic model of three equations for the fluid depth, the average velocity and a third variable, called entrophy, related to the variance of the velocity. The model features new relaxation source terms and admits an exact balance energy equation. The velocity field in the depth is consistently reconstructed using only the variables of the depth-averaged model without any derivative. The physical relevance of the entrophy is related to the shape of the velocity profile. The linear stability of a uniform solution is investigated for this model, showing a stabilizing effect of the plasticity. Roll waves are simulated numerically using a classical Godunov's scheme. The model for a Newtonian fluid is presented as a particular case.

## 1. Introduction

Viscoplastic materials behave like solid bodies when the exerted stress is less than a certain threshold (the yield stress), and flow like viscous fluids above this threshold. Such materials are encountered in various contexts including biological fluids (blood clots, mucus), industrial processes (cement, waxy crude oil), and geophysical flows (avalanches, debris and mud flows) [1]. The development of accurate models for describing free-surface flows of such viscoplastic materials is of great importance for applications such as ink-jet printing or to better predict natural hazards [2]. In the present work, we consider gravity-driven free-surface flows of idealized viscoplastic fluids propagating down an inclined plane. By idealized viscoplasticity, we refer to a perfectly rigid behaviour in the solid-like regime [3].

The mathematical modeling of idealized viscoplastic fluids generally relies on the Herschel-Bulkley or Bingham constitutive laws. Combined with the Cauchy momentum equations, such constitutive laws can be used to compute fluid flows. However, direct numerical simulation (DNS) of viscoplastic flows is not a straightforward task, notably due to the complexity involved in identifying the yield surfaces separating unyielded (solid-like) from yielded (fluid-like) regions. Two main methods have been developed to treat this issue. The regularization method consists in replacing the rigid behaviour in the unyielded zones by a highly viscous flow [4]. The Augmented Lagrangian method introduces a reformulation of the Cauchy momentum equations into a variational form to compute the flow as the solution of an optimization problem [5-8]. In either cases, accurate DNS of viscoplastic flows generally requires large computing times and the use of very fine meshes in the vicinity of yield surfaces.

For free-surface flows, an alternative to DNS is to derive models of reduced dimensionality. The most common approach is based on a thin-layer approximation, which together with averaging the Cauchy momentum equations over the depth of the flow forms the basis of numerous reduced-order models used in hydraulics [9, 10]. Another benefit of the depth-averaged approach is that the boundary conditions are directly incorporated into the model, thus allowing for easier and faster numerical resolution. Formally, the derivation of thin-layer models is generally based on two steps. The first step consists in obtaining long-wave asymptotic expansions of the flow variables with respect to the flow aspect ratio $\epsilon = h_0/l_0$, where $h_0$ and $l_0$ denote the typical depth and length of the flow, respectively. The second step consists in averaging the governing equations over the depth of the flow and rewriting the resulting system in terms of

<footnote>*Corresponding author (danila.denisenko@inrae.fr)
ORCID(s): 0000-0003-1293-2633 (D. Denisenko)</footnote>
averaged quantities. To capture the right physics and properly account for the fluid rheology, the derived models should be consistent at least at order 1. A model is said to be consistent at order *n* if the leading terms in the model equations are of *O*(1), and if all terms vanish except for a remainder of *O*(ε^(n+1)) after having inserted the asymptotic expansions obtained in the first step into the model equations. Inconsistent reduced models lead to inaccurate predictions of, e.g., instability thresholds [11, 12].

Depth-averaged models can be formulated as systems of one, two or three equations. One-equation models are obtained from the mass conservation by enslaving the fluid velocity to the fluid height [13–18]. However, consistent one-equation models generally produce diverging or inaccurate solutions when the instability threshold for uniform flows is exceeded [19–21]. Two-equation models introduce the averaged velocity as a second independent variable. The second equation can be based on either the momentum balance or the work-energy theorem. Such consistent two-equation models have been derived for a variety of Newtonian and power-law fluids and were shown to provide accurate results in many applications [11, 14, 22–25]. Although relatively rare, a few studies also considered the case of viscoplastic fluids [12, 26]. However, the mathematical structure and numerical resolution of consistent two-equation models can be complicated, notably in the case of sheared flows. Moreover, Richard et al. [27] showed that two-equation models based on the depth-averaged momentum equation are not compatible with the work-energy theorem, and vice-versa. To ensure Galilean invariance and compatibility between the depth-averaged momentum and energy equations, Richard et al. [28] derived a three-equation model for Newtonian fluids by introducing a third variable, called enstrophy, related to deviation of the velocity from its averaged value. An important benefit of this three-equation approach is that the resulting system has the mathematical structure of the Euler equations of compressible fluids, which ensures the well-posedness of the problem and guarantees an efficient computational resolution with reliable numerical schemes.

Although the derivation of long-wave asymptotic expansions is relatively straightforward for Newtonian or power-law fluids, for viscoplastic fluids difficulties arise from the possible coexistence within the flows of yielded and unyielded regions, the latter corresponding to rigid plugs. However, in many flow configurations, the region that appears to be a rigid plug is in fact very slightly above the yield point and flows almost imperceptibly. Such a region is called a pseudo-plug (Walton & Bittleston (1991) [29]). At leading order with respect to ε, the asymptotic expansion of longitudinal velocity describes a yielded layer at the base of the flow, overlaid by an unyielded plug zone close to the free surface [30, 31]. At the next order of approximation, Balmforth and Craster [32] showed that to have a consistent long-wave theory, this plug layer has to be treated as a pseudo-plug in which the strain-rate is of order O(ε). These authors derived a first-order correction for the longitudinal velocity profile in inertia-less limit. Later, Chambon et al. [33] constructed the full expressions for the longitudinal velocity up to the order O(ε). However, the obtained asymptotic solution shows two main drawbacks: (1) the strain rate in the pseudo-plug becomes infinite at the fake yield surface (i.e., the interface between the pseudo-plug and the sheared layer), leading to an unphysical kink in the velocity profile; (2) the viscous contribution of the rheology does not contribute to the shearing of the pseudo-plug at first order. For Bingham fluids, Balmforth and Craster [32] and Fernández-Nieto et al. [26] proposed to avoid the divergence of the strain rate by introducing a transition layer between the pseudo-plug and the sheared layer. Fernández-Nieto et al. [26] also derived a consistent two-equation model. However, the extra terms arising from the transition layer are of order O(ε^4/3) and cannot be given by analytical expressions, such that the expansions providing the smooth transition at the fake yield surface were not considered in deriving this model. At the order O(ε), the strain rate at the bottom of the pseudo-plug layer remains singular, such that the possibility to reconstruct smooth velocity profiles in a consistent manner from this model is not warranted. Moreover, the shallow-flow model of [26] does not admit an equation of energy conservation. This is a common problem for all two-equation models involving a coefficient other than 1 in front of the term in *hU*² in the momentum flux (see Richard et al. [27] for the proof in the Newtonian case). Its hyperbolicity is also not completely guaranteed, due to a term proportional to |∂h/∂x| (*h* being the flow depth). As a consequence, the numerical resolution of this model is not straightforward, and no proof of stability can be given for the numerical scheme. A three-equation model would solve these issues, since the mathematical structure of the Euler equations of compressible fluids with source terms without derivatives can be recovered. This guarantees the energy conservation, the hyperbolicity and the existence of a mathematical entropy, which is needed for a hyperbolic system of equations in conservative form. With this structure, the numerical resolution is straightforward and well-known numerical schemes can be implemented.

The goal of this paper is to generalize the three-equation approach of Richard et al. [28] to the case of Bingham fluid flows propagating down an inclined plane, in order to derive a consistent depth-averaged model with a well-posed mathematical structure. We construct a new asymptotic solution up to the first order in ε, based on a specific regularized version of the tensorial constitutive law. This allows us to eliminate the issue with the diverging strain rate in the
pseudo-plug without the need to introduce a third layer in the model. A consistent depth-averaged model is then derived by averaging the mass, momentum, and energy conservation equations over the fluid depth, introducing an entrophy variable. The resulting model is a fully hyperbolic system with relaxation source terms, whose computational resolution can be handled by robust numerical schemes. An analysis of the linear stability of the derived model demonstrates the stabilizing effect of plasticity and shows good agreement with the instability criterion obtained by Balmforth and Liu [12] from the generalized Orr-Sommerfeld equations. Another important advantage of the three-equation approach is also highlighted, namely that the full velocity field can be consistently reconstructed at order $O(\epsilon)$ directly from the variables (flow height, averaged velocity and entrophy) of the model, without any derivatives. In particular, we analyze the physical relevance of the entrophy in terms of shapes of the velocity profiles within a roll wave.

In §2, we present the regularization of the constitutive law and formulate the equations for the fluid flow. In §3, we construct the new shallow-flow asymptotic expansion up to $O(\epsilon)$ order. In §4, the consistent three-equation model is derived by averaging the mass, momentum and energy balances. In §4, the velocity field is reconstructed from the variables of the model. Finally, in §5 we investigate the linear stability of the uniform flow, present numerical simulations of roll waves and discuss the physical relevance of the entrophy variations predicted by the model.

## 2. Formulation of the problem

We consider a two-dimensional flow of a viscoplastic fluid propagating down an inclined plane under gravity $g$ (Figure 1). The angle of the slope with respect to the horizontal is $\theta$. The directions $Ox$ and $Oz$ are parallel and normal to the plane, respectively. The corresponding components of the velocity field $v$ are denoted by $u$ and $w$, and the components of the strain-rate tensor $\dot{\gamma}$ are defined as: $\dot{\gamma}_{xx} = 2\partial u/\partial x$, $\dot{\gamma}_{xz} = \partial u/\partial z + \partial w/\partial x$, $\dot{\gamma}_{zz} = 2\partial w/\partial z$. With this definition, the strain-rate tensor is related to the symmetrical part $D$ of the velocity gradient by $\dot{\gamma} = 2D$. The fluid depth is denoted by $h(x, t)$. Lastly, the fluid is assumed to be incompressible (tr $\dot{\gamma} = 0$) with a density $\rho$.

<figure><img src="image_3.png" /><figcaption>Figure 1: Definition sketch</figcaption></figure>

### 2.1. Constitutive law

The fluid is assumed to obey the Bingham constitutive law. The 1D Bingham law, assuming a simple shear, is written $\tau = \tau_c + K\dot{\gamma}$ above the yield point, where $\tau_c$ and $K$ are the yield stress and the Bingham viscosity of the material, respectively. Hohenemser & Prager (1932) [34] were the first to propose a 3D tensorial extension of this law, which is commonly used nowadays. Above the yield point, the deviatoric part $\tau$ of the stress tensor is written

$$ \tau = \tau_c \frac{\dot{\gamma}}{|\dot{\gamma}|} + K\dot{\gamma}. \qquad (2.1) $$

The tensor norm is defined as $|T| = (T : T/2)^{1/2}$ for any second-order tensor $T$, where the colon denotes the double dot product. Oldroyd (1947) [35] modelled the solid behaviour below the yield point as a linear elastic solid and the liquid behaviour above the yield point as a viscous fluid with the same tensor extension as Hohenemser & Prager (1932) [34]. The yield point is defined with a von Mises criterion, namely $|\tau| = \tau_c$ and $\dot{\gamma} = 0$. In this approach the deviatoric stress tensor is related to the strain tensor $\gamma$ below the yield point, and to the strain-rate tensor $\dot{\gamma}$ above the yield point. Consequently, there is a discontinuity of the stress tensor at the yield point [36, 37].
To avoid the discontinuity of the stress tensor at the yield point, a perfectly rigid behaviour of the material below the yield point is commonly assumed. In this case, the law (2.1) is assumed if $|\tau| > \tau_c$, and $\dot{\gamma} = 0$ if $|\tau| \le \tau_c$. This implies that the stress tensor is indeterminate below the yield point, and determinate above the yield point with an infinite apparent viscosity at the yield point, since $\tau_c/|\dot{\gamma}| \to \infty$. Various approaches, such as regularizing methods, have been proposed to avoid both this indeterminacy and the divergence of the viscosity [1, 4].

It is important to highlight that the expression (2.1) for the deviatoric stress tensor, and in particular the form of the first term $\tau_c \dot{\gamma}/|\dot{\gamma}|$, called hereafter the yield-stress tensor, are in no way necessary consequences of the 1D law. This tensorial expression results from an additional assumption, which was clearly formulated by Oldroyd (1947) [35]. This author assumed that the constitutive law, at any point where flow is occurring, can be reduced to three 1D equations corresponding to three superposed simple shear flows, in a well-chosen system of axes. This assumption is actually equivalent to stating that the stress tensor is aligned with the strain-rate tensor.

From a physical point of view, this assumption of alignment between the two tensors, and the validity of the constitutive law (2.1) for general 3D flows, remain to be proved (Coussot [38]). Experiments performed in pure elongational flows tend to display contrasted conclusions regarding the validity of the yield point predicted by (2.1) in this configuration [39, 40]. The assumption of alignment between the stress and strain-rate tensors is clearly called into question by, e.g., rheometrical measurements obtained in large amplitude oscillatory shear [41] or measurements of normal stress components at yielding [42, 43]. However, the interpretation of those results is generally complicated by the fact that most real yield-stress fluids present additional rheological properties that are not captured by idealized viscoplastic models, such as elasticity or thixotropy.

In the frame of asymptotic methods, using the constitutive law (2.1) also leads to several mathematical complications, such as the divergence of the strain rate in the pseudo-plug layer at order $O(\epsilon)$ and the need to introduce a transition layer based on matched asymptotic expansions (see introduction). Furthermore, due to the term in $1/|\dot{\gamma}|$, the strain rate in the pseudo-plug at order $O(\epsilon)$ is given solely by the plastic contribution to the stress and does not involve the viscous contribution, which appears counterintuitive (see §3.2.1 below) and leads to a nonphysical destabilizing effect of plasticity at large slopes (see §6.2 below).

Instead of the classical formulation, and in order to derive a more tractable asymptotic model, we thus chose to propose an alternative tensor expression of the constitutive equation. Namely, we assume that

$$ \tau = \tau^Y + K\dot{\gamma} \qquad (2.2) $$

where $\tau^Y$ is a yield-stress tensor that is not, in general, aligned with the strain-rate tensor $\dot{\gamma}$. Hence, the strain rate $\dot{\gamma}$ is aligned with the “overstress” $\tau - \tau^Y$ and not with $\tau$ itself. In general, the tensor $\tau^Y$, which is analogous to the backstress introduced in certain generalized viscoplastic models [37], depends on internal variables of the material. For the purpose of our model, it is not necessary to fully specify this tensor, although this would of course be required to compute general flows. We only make here the minimum assumptions necessary to derive our model, namely:

- The norm of $\tau^Y$ is always equal to the yield stress $\tau_c$: $|\tau^Y| = \tau_c$. This allows us to recover the von Mises yielding criterion.

- The trace of $\tau^Y$ is zero, due to the incompressibility hypothesis: $\text{tr}\tau^Y = 0$.

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

- The existence of a normal stress difference in the pseudo-plug at leading order, and all terms originating from these normal stresses in the asymptotic expansions, are attributed to the yield-stress tensor and not to the viscous-stress tensor.

As will be shown below, these assumptions on $\tau^Y$ are sufficient to work out a consistent asymptotic model at $O(\epsilon)$ order, which has the benefit of alleviating the drawbacks of the classical formulation.

This alternative tensor expression can be seen as a way to regularize the constitutive law at the yield point. In particular, the divergence of the apparent viscosity and the singularity of the yield-stress tensor for a vanishing strain rate are removed. This regularization is however very different from the classical regularizations used for viscoplastic fluids, in which the apparent viscosity is assumed to saturate at a very large value [4]. Such regularizations imply that the material behaves as a viscous liquid both above and below the yield point. In our approach, on the contrary, the suppression of the infinite effective viscosity is formally achieved by removing the alignment of the yield-stress tensor
with the strain-rate tensor, which enables to keep a true yield point. It is important to note that the classical formulation (2.1) cannot be recovered as a particular case of our approach, since we do not assume that the yield-stress tensor $\tau^Y$ can be expressed through an apparent viscosity inversely proportional to $|\dot{\gamma}|$.

## 2.2. Governing equations

The fluid motion is governed by the Cauchy mass and momentum conservation equations, completed by boundary conditions on the bottom wall and at the free surface. The total stress tensor can be written $\sigma = -p\mathbf{I} + \tau$, where $p$ is the pressure, $\mathbf{I}$ is the identity tensor, and the deviatoric part $\tau$ is given by (2.2). The continuity equation is

$$ \frac{\partial u}{\partial x} + \frac{\partial w}{\partial z} = 0. \qquad (2.3) $$

The Cauchy momentum equations in the $Ox$ and $Oz$ directions are

$$ \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}, \qquad (2.4) $$

$$ \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_{zz}}{\partial z}. \qquad (2.5) $$

According to the constitutive equation (2.2) and the conditions on the yield-stress tensor indicated in §2.1, the stress components, for $|\tau| > \tau_c$, can be written as

$$ \tau_{xx} = -\tau_{zz} = \tau_{xx}^Y + \tau_{xx}^v, \qquad \tau_{xx}^v = 2K \frac{\partial u}{\partial x}, \qquad (2.6) $$

$$ \tau_{xz} = \tau_{xz}^Y + \tau_{xz}^v, \qquad \tau_{xz}^v = K \left( \frac{\partial u}{\partial z} + \frac{\partial w}{\partial x} \right), \qquad (2.7) $$

with the relation for the norm of the yield stress tensor:

$$ (\tau_{xx}^Y)^2 + (\tau_{xz}^Y)^2 = \tau_c^2. \qquad (2.8) $$

At the bottom we consider the no-penetration and the no-slip conditions:

$$ u_{|z=0} = w_{|z=0} = 0. \qquad (2.9) $$

At the free surface $z = h(x)$, the following kinematic boundary condition holds:

$$ \frac{\partial h}{\partial t} + u_{|z=h(x)} \frac{\partial h}{\partial x} = w_{|z=h(x)}. \qquad (2.10) $$

Lastly, capillarity is neglected and the atmospheric pressure is assumed to be constant and taken equal to zero. Accordingly, the free surface is stress-free and the dynamic boundary conditions can be written

$$ \left[ 1 - \left( \frac{\partial h}{\partial x} \right)^2 \right] \tau_{xz|_{z=h(x)}} = 2 \frac{\partial h}{\partial x} \tau_{xx|_{z=h(x)}}, \qquad (2.11) $$

$$ \left[ 1 - \left( \frac{\partial h}{\partial x} \right)^2 \right] p_{|_{z=h(x)}} = - \left[ 1 + \left( \frac{\partial h}{\partial x} \right)^2 \right] \tau_{xx|_{z=h(x)}}. \qquad (2.12) $$

## 2.3. Shallow-flow scaling

Let us define $h_0$ the characteristic depth of the flow in the $Oz$ direction, and $u_0$ the characteristic velocity in the $Ox$ direction. The characteristic length in the $Ox$ direction is denoted by $l_0$. The shallow-flow hypothesis corresponds to assuming that the aspect ratio $\epsilon = h_0/l_0$ is small. The main dimensionless groups of this problem are the Reynolds number $Re$, the Froude number $Fr$ and the Bingham number $Bi$, which are defined as

$$ Re = \frac{\rho u_0 h_0}{K}, \qquad Fr = \frac{u_0}{\sqrt{gh_0 \cos \theta}}, \qquad Bi = \frac{\tau_c h_0}{K u_0}. \qquad (2.13) $$
These parameters, as well as the slope angle $\theta$, are assumed to be of order $O(1)$ with respect to aspect ratio $\epsilon$. In order to reformulate the problem (2.3)-(2.12) into a dimensionless form, let us rescale the variables as follows:

$$ 
\begin{aligned}
x &= l_0 \bar{x}; & z &= h_0 \bar{z}; & u &= u_0 \bar{u}; & w &= \epsilon u_0 \bar{w}; & t &= \frac{l_0}{u_0} \bar{t}; \\
h &= h_0 \bar{h}; & p &= \rho g h_0 \cos \theta \bar{p}; & \tau_{ij} &= \tau_c \bar{\tau}_{ij}; & |\boldsymbol{\tau}| &= \tau_c |\bar{\boldsymbol{\tau}}|; & |\dot{\boldsymbol{\gamma}}| &= \frac{u_0}{h_0} |\dot{\bar{\boldsymbol{\gamma}}}|
\end{aligned} 
\quad (2.14) 
$$

Omitting the bars, the dimensionless continuity equation (2.3) keeps the same form

$$ \frac{\partial u}{\partial x} + \frac{\partial w}{\partial z} = 0, \quad (2.15) $$

while the dimensionless momentum equations (2.4)-(2.5) now are rewritten as

$$ \epsilon \left( \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + w \frac{\partial u}{\partial z} \right) = - \frac{\epsilon}{Fr^2} \frac{\partial p}{\partial x} + \frac{\lambda}{Re} + \frac{Bi}{Re} \frac{\partial \tau_{xz}}{\partial z} + \frac{\epsilon Bi}{Re} \frac{\partial \tau_{xx}}{\partial x}, \quad (2.16) $$

$$ \epsilon^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{\epsilon Bi}{Re} \frac{\partial \tau_{xz}}{\partial x} + \frac{Bi}{Re} \frac{\partial \tau_{zz}}{\partial z}, \quad (2.17) $$

with the driving parameter $\lambda$ given by

$$ \lambda = \frac{\rho g h_0^2 \sin \theta}{K u_0} = \frac{Re}{Fr^2} \tan \theta. \quad (2.18) $$

The yielding criterion becomes $|\boldsymbol{\tau}| = 1$, and the expressions for the viscous and yield-stress components of the stress (2.7)-(2.8), for $|\boldsymbol{\tau}| > 1$, are transformed to

$$ \tau_{xx} = \tau_{xx}^{Y} + \tau_{xx}^{v}, \qquad \tau_{xx}^{v} = 2 \frac{\epsilon}{Bi} \frac{\partial u}{\partial x}, \quad (2.19) $$

$$ \tau_{xz} = \tau_{xz}^{Y} + \tau_{xz}^{v}, \qquad \tau_{xz}^{v} = \frac{1}{Bi} \left( \frac{\partial u}{\partial z} + \epsilon^2 \frac{\partial w}{\partial x} \right), \quad (2.20) $$

with the condition on the norm of the dimensionless yield-stress tensor:

$$ (\tau_{xx}^Y)^2 + (\tau_{xz}^Y)^2 = 1. \quad (2.21) $$

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

$$ u_{|z=0} = w_{|z=0} = 0, \quad (2.22) $$

while the dimensionless kinematic and dynamic boundary conditions at the free surface (2.10)-(2.12) are rewritten as

$$ \frac{\partial h}{\partial t} + u_{|z=h(x)} \frac{\partial h}{\partial x} = w_{|z=h(x)}, \quad (2.23) $$

$$ \left[ 1 - \epsilon^2 \left( \frac{\partial h}{\partial x} \right)^2 \right] \tau_{xz|_{z=h(x)}} = 2\epsilon \frac{\partial h}{\partial x} \tau_{xx|_{z=h(x)}}, \quad (2.24) $$

$$ \left[ 1 - \epsilon^2 \left( \frac{\partial h}{\partial x} \right)^2 \right] p_{|_{z=h(x)}} = - \frac{Bi Fr^2}{Re} \left[ 1 + \epsilon^2 \left( \frac{\partial h}{\partial x} \right)^2 \right] \tau_{xx|_{z=h(x)}}. \quad (2.25) $$

Lastly, the norms of the dimensionless stress and strain-rate tensors express as

$$ |\boldsymbol{\tau}| = \sqrt{\tau_{xx}^2 + \tau_{xz}^2} \qquad |\dot{\boldsymbol{\gamma}}| = \sqrt{\left(\frac{\partial u}{\partial z} + \epsilon^2 \frac{\partial w}{\partial x}\right)^2 + 4\epsilon^2 \left(\frac{\partial u}{\partial x}\right)^2}. \quad (2.26) $$

### 3. Asymptotic expansions

Let us assume the existence of regular expansions of the form

$$ f = f^{(0)} + \epsilon f^{(1)} + \dots \quad (3.1) $$

for all variables of the problem (2.15)-(2.25), namely the longitudinal and normal velocities $u$ and $w$, the pressure $p$ and the stress components $\tau_{ij}$.
## 3.1. Leading-order expansion

### 3.1.1. Pseudo-plug thickness

At $O(1)$ with respect to $\epsilon$, the integration of the momentum equation (2.16) with the boundary condition (2.24) gives the following shear-stress profile:

$$ \tau_{xz}^{(0)} = \frac{\lambda}{Bi}(h-z) \qquad (3.2) $$

If the material is in the yielded regime at this order, then $|\tau|^{(0)} > 1$. Following the assumptions discussed in §2.1, the normal stress is zero in the yielded case at the leading-order since the flow is a simple shear flow (see below):

$$ \tau_{xx}^{(0)} = 0. \qquad (3.3) $$

This implies that $|\tau|^{(0)} = \tau_{xz}^{(0)}$ (since $z \le h$). This behaviour is possible as long as $\tau_{xz}^{(0)} > 1$, which is equivalent to $z < h - Bi/\lambda$. This condition defines two layers: a sheared layer for $z < h - Bi/\lambda$, where the material is above the yield point, and a plug layer for $z \ge h - Bi/\lambda$, where, at order 0, $|\dot{\gamma}^{(0)}| = 0$.

However, as long as the flow is not uniform and stationary, the existence of a true plug is not possible. Following the approach of Walton and Bittleston [29] and Balmforth and Craster [32], we assume a non-zero strain rate at order 1 with respect to $\epsilon$. Hence, the fluid is actually slightly above the yield point for $z \ge h - Bi/\lambda$, because of the first-order correction. This layer is hereafter called a pseudo-plug layer. Its thickness is given by

$$ h_p = \frac{Bi}{\lambda}. \qquad (3.4) $$

In practice, the asymptotic expansions have to be constructed in the sheared and the pseudo-plug layers separately. In the following, tilde notations will be used for the expansions of the variables in the pseudo-plug layer.

In the pseudo-plug layer, since $|\dot{\gamma}^{(0)}| = 0$, the leading-order longitudinal velocity $\tilde{u}^{(0)}$ is independent of the normal coordinate $z$ and the expansion can be written

$$ u = \tilde{u}^{(0)}(x, t) + \epsilon \tilde{u}^{(1)}(x, z, t) + O(\epsilon^2). \qquad (3.5) $$

The strain rate thus reads

$$ |\dot{\gamma}| = \epsilon \sqrt{\left(\frac{\partial \tilde{u}^{(1)}}{\partial z}\right)^2 + 4 \left(\frac{\partial \tilde{u}^{(0)}}{\partial x}\right)^2} + O(\epsilon^2) \qquad (3.6) $$

while the shear stress at order 0 has the same form as in the sheared layer

$$ \tilde{\tau}_{xz}^{(0)} = \frac{\lambda}{Bi}(h-z). \qquad (3.7) $$

It is important to note that this definition of a pseudo-plug layer of constant thickness $h_p$ is obtained only in the framework of the asymptotic expansion. It defines a fake yield surface at which the asymptotic expansions in both layers must be matched. This pseudo-plug layer is completely different from the physical pseudo-plug layer that could be measured in an experimental investigation. In experiments, as there is no true plug, a pseudo-plug layer can only be defined with respect to some critical value $\dot{\gamma}_c$ of the strain rate, which can be for example a typical level of experimental noise [33]. In a similar way, the model derived below allows to calculate an effective pseudo-plug layer from the reconstructed velocity field (see §5) if a critical strain rate $\dot{\gamma}_c$ is chosen. The thickness of this effective pseudo-plug layer is certainly not constant, and is different from the formal pseudo-plug thickness $h_p$ defined from the leading order of the asymptotic expansion.

In some previous studies [12, 33], the pseudo-plug used for the asymptotic method was considered to have a varying thickness. This resulted from assuming that the Froude number is small ($Fr^2 = O(\epsilon)$) and that the slope angle is also small ($\tan \theta = O(\epsilon)$) [33]. These assumptions imply that the pressure term in the momentum balance equation in the $Ox$-direction (2.16) is of $O(1)$, which results in a dependence of the pseudo-plug thickness $h_p$ with $\partial h/\partial x$. In contrast,
in the approach presented here, the pressure is taken into account only at $O(\epsilon)$. Moreover, in the case of a small Froude number and of a small slope, there is a corrective term to the pressure at $O(\epsilon)$, which is not found in the present approach (where it would be found at $O(\epsilon^2)$). This corrective term includes a second derivative of the fluid depth $\partial^2 h/\partial x^2$ and a nonlinear term $(\partial h/\partial x)(\partial^2 h/\partial x^2)$, which significantly complicate the expansions and prevent from deriving a hyperbolic model. Note that, while the leading order term of the pressure is certainly of major importance, this may not be the case for the corrective term. In the present approach, the reconstructed velocity field (see §5) and the effective (variable) pseudo-plug thickness can only be derived by taking into account all corrections of $O(\epsilon)$, since the pressure appears at this order of accuracy. As it is anyway necessary to include the terms of $O(\epsilon)$ to derive a consistent model, we expect that both approaches should give at the end very similar results in terms of thickness of the effective pseudo-plug layer. The advantage of the present approach, without the corrective term of the pressure, is that the asymptotic expansions and the derivation of a well-posed model are much simpler.

### 3.1.2. In the pseudo-plug ($z \ge h - h_p$)

Owing to the expansion (3.5), the constitutive law (2.19)-(2.20) reduces to

$$ \tilde{\tau}_{xx}^{(0)} = \tilde{\tau}_{xx}^{Y(0)}, \quad (3.8) $$

$$ \tilde{\tau}_{xz}^{(0)} = \tilde{\tau}_{xz}^{Y(0)}. \quad (3.9) $$

From the condition on the norm of the yield-stress tensor (2.21), we then obtain

$$ \tilde{\tau}_{xx}^{Y(0)} = \delta \sqrt{1 - \left(\frac{h-z}{h_p}\right)^2} \quad (3.10) $$

with $\delta = \operatorname{sgn}(\tau_{xx})$. Hence, it is found that the normal stresses contribute at leading-order in the pseudo-plug, to ensure that the layer is just at the verge of yielding, $|\tilde{\tau}|^{(0)} = 1$. Note that these normal stresses vanish at the fake yield surface $z = h - h_p$. Further, the pressure profile is obtained from integration of the momentum equation (2.17) with the dynamic condition (2.25):

$$ \tilde{p}^{(0)} = h - z - \frac{BiFr^2}{Re} \tilde{\tau}_{xx}^{Y(0)}. \quad (3.11) $$

The expression for the longitudinal velocity $\tilde{u}^{(0)}(x, t)$ will be given later from matching with the sheared zone. The normal velocity $\tilde{u}^{(0)}$ can then be derived by integration of the continuity equation (2.15):

$$ \tilde{u}^{(0)} = -z \frac{\partial \tilde{u}^{(0)}}{\partial x} + w_{+}^{(0)}, \quad (3.12) $$

where, again, the term $w_{+}^{(0)}$ will be obtained from matching.

### 3.1.3. In the sheared zone ($z < h - h_p$)

In this zone, the constitutive law expresses as

$$ \tau_{xx}^{(0)} = \tau_{xx}^{Y(0)}, \quad (3.13) $$

$$ \tau_{xz}^{(0)} = \tau_{xz}^{Y(0)} + \frac{1}{Bi} \frac{\partial u^{(0)}}{\partial z}. \quad (3.14) $$

The condition on the smallness of the normal stress in the sheared layer (3.3) leads to $\tau_{xx}^{Y(0)} = 0$. Then the von Mises criterion (2.21) implies that $\tau_{xz}^{Y(0)} = 1$ and integration of equation (3.14) with the no-slip condition (2.22) results in a parabolic longitudinal velocity profile:

$$ u^{(0)} = \lambda z \left(h - h_p - \frac{z}{2}\right). \quad (3.15) $$

Note that $\partial u^{(0)}/\partial z = 0$ at $z = h - h_p$, consistently with the expansion (3.5) assumed in the pseudo-plug. Note also that in the sheared layer, the assumption $\tau_{xx}^{Y(0)} = 0$ implies that the yield-stress tensor is aligned with the tensor $\dot{\gamma}/|\dot{\gamma}|$. This allows us to recover the 1D Bingham law in this particular case of a simple shear flow.
The normal velocity $w^{(0)}$ is then found from the integration of equation (2.15) with the no-slip boundary condition (2.22):

$$w^{(0)} = \frac{-\lambda z^2}{2} \frac{\partial h}{\partial x}. \qquad (3.16)$$

Finally, the integration of the momentum equation (2.17) and the stress continuity at the fake yield surface $z = h - h_p$ provide the pressure profile:

$$p^{(0)} = h - z. \qquad (3.17)$$

Note that, unlike in the pseudo-plug, a classical hydrostatic pressure distribution is recovered in the sheared layer.

### 3.1.4. Matching

The continuity condition $\tilde{u}^{(0)} = u^{(0)}$ at the fake yield surface $z = h - h_p$ provides the expression of the longitudinal velocity in the pseudo-plug:

$$\tilde{u}^{(0)} = \frac{\lambda}{2} (h - h_p)^2. \qquad (3.18)$$

It is worth noting that the leading-order longitudinal velocity profile given by Eqs. (3.15) and (3.18) is identical to the profile that would be obtained in a steady uniform flow of height $h$, namely a parabolic profile overlaid by a "true" unsheared plug. As will be shown in the next section, however, the plug effectively becomes a slightly-sheared pseudo-plug at order $O(\epsilon)$.

Finally, matching the solutions $w^{(0)}$ and $\tilde{w}^{(0)}$ at $z = h - h_p$ gives the following expression for the normal velocity in the pseudo-plug:

$$\tilde{w}^{(0)} = -\lambda(h - h_p) \left(z - \frac{h - h_p}{2}\right) \frac{\partial h}{\partial x}. \qquad (3.19)$$

## 3.2. $O(\epsilon)$ expansion

Here we construct the expansions at order $O(\epsilon)$ for the shear stress $\tau_{xz}^{(1)}$ and the longitudinal velocity $u^{(1)}$. The expansions of the other variables will not be needed for the derivation of a depth-averaged model consistent at first order. Note that the leading-order solution derived above does not depend on whether one considers the classical extension (2.1) of Bingham's rheology, or the formulation (2.2) with the assumptions detailed in §2.1. This does not remain true, however, at $O(\epsilon)$.

### 3.2.1. Determination of the yield-stress tensor

With the classical constitutive law (2.1), equations (3.6) and (3.9) lead to the following relation in the pseudo-plug:

$$\tilde{\tau}_{xz}^{(0)} = \tilde{\tau}_{xz}^{Y(0)} = \frac{\frac{\partial \tilde{u}^{(1)}}{\partial z}}{\sqrt{\left(\frac{\partial \tilde{u}^{(1)}}{\partial z}\right)^2 + 4 \left(\frac{\partial \tilde{u}^{(0)}}{\partial x}\right)^2}} \qquad (3.20)$$

and thus:

$$\frac{\partial \tilde{u}^{(1)}}{\partial z} = \frac{2(h-z)/hp}{\sqrt{1 - ((h-z)/h_p)^2}} \left| \frac{\partial \tilde{u}^{(0)}}{\partial x} \right|. \qquad (3.21)$$

Expression (3.21) was used in former studies to derive the first-order velocity correction $\tilde{u}^{(1)}$ in the pseudo-plug [26, 32, 33]. However, it is easily seen that $\partial \tilde{u}^{(1)} / \partial z$ diverges at the fake-yield surface $z = h - h_p$, which contradicts the assumption that the strain-rate should remain small in the pseudo-plug. As already suggested by Balmforth and Craster [32] and further explored by Fernández-Nieto et al. [26], this inconsistency could be alleviated by introducing a
transition layer between the sheared zone and the pseudo-plug using asymptotic matching. However, this significantly complicates the process of constructing the solution, due in particular to the appearance of terms of $O(\epsilon^{4/3})$ in the velocity profile [26]. Furthermore, the terms arising from the transition layer are not amenable to analytical expressions.

As already mentioned, another drawback of expressions (3.20) (3.21) is that the first-order correction $\partial\tilde{u}^{(1)}/\partial z$ is controlled solely by terms related to the yield-stress tensor. Hence, the implementation of the asymptotic method is completely different in the sheared layer and in the pseudo-plug. In the sheared layer, the first-order correction to $\tau_{xz}$ is found from the Ox-momentum equation written at order 1. The first-order correction $u^{(1)}$ to the velocity is then obtained from the constitutive law by integration of the viscous part of $\tau_{xz}^{(1)}$, since the yield-stress component is equal to $1 + O(\epsilon^2)$ and therefore does not appear in this calculation. In contrast, in the pseudo-plug zone, the component $\tau_{xz}^{(0)}$ of order 0 is used, and the first-order correction to the velocity is found from the constitutive law by a calculation involving only the yield-stress term. The viscous term does not appear in this calculation, because the effective viscosity associated to the yield-stress term becomes infinite at the yield point, and thus dominates over the Bingham viscosity. However, above the yield point, the material behaves as a liquid, and it might thus seem more physical to assume that this viscous term should be involved in the determination of the first-order correction to the velocity, as in the sheared layer. In other words, it could be expected that the slight shearing of the pseudo-plug be associated to viscous stresses, rather than to yield-stress terms.

The method proposed here, based on the 3D extension to the Bingham rheology discussed in §2.1, avoids any divergence of the strain rate and dispenses with the need to introduce a transition layer. Thus, the implementation of the asymptotic method follows the same algorithm in the sheared layer and in the pseudo-plug layer. Removing the assumption that the yield-stress and strain-rate tensors are aligned implies that the yield-stress tensor is well defined and can be regularly expanded even if the strain rate is zero. In particular, its xz-component can be expanded in the pseudo-plug as

$$ \tau_{xz}^{Y} = \tilde{\tau}_{xz}^{Y(0)} + \epsilon \tilde{\tau}_{xz}^{Y(1)} + O(\epsilon^2) \qquad (3.22) $$

where the term $\tilde{\tau}_{xz}^{Y(0)}$ is obtained from the zero-order solution in the pseudo-plug, i.e. for $\dot{\gamma}^{(0)} = 0$, and given by the formula (3.9). In the pseudo-plug, when the strain rate is small and of $O(\epsilon)$, there is a small correction of $O(\epsilon)$ to the xz-component of the yield-stress tensor, and it is $\tilde{\tau}_{xz}^{Y(1)}$. Consequently, this term, added to the viscous part, is equal to the first-order correction in the asymptotic expansion of $\tau_{xz}$ as

$$ \tilde{\tau}_{xz}^{(1)} = \tilde{\tau}_{xz}^{Y(1)} + \frac{1}{Bi} \frac{\partial \tilde{u}^{(1)}}{\partial z} \qquad (3.23) $$

Therefore, it is seen that the first-order correction to the velocity is obtained from the first-order shear stress $\tilde{\tau}_{xz}^{(1)}$, and not from the zero-order shear stress $\tilde{\tau}_{xz}^{(0)}$ as in the classical approach (see (3.20)). Moreover, contrary to the classical approach, the viscous stress contribution is involved in this calculation.

Again, note that the classical expansion (3.20) and (3.21) cannot be recovered as a subset of our approach since, as soon as the infinite apparent viscosity for a zero strain rate is removed, the first-order correction to the velocity is obtained from the first-order correction to the shear stress, rather than from the zero-order term. This also explains why our approach can be seen as a specific type of regularization of the constitutive law, albeit very different from conventional regularizations (see §2.1). A common feature of all regularized approaches is that the first-order correction to the velocity is obtained from the first-order correction to the shear stress instead of from the zero-order shear stress.

### 3.2.2. In the pseudo-plug ($z > h - h_p$)

The momentum balance equation in the Ox-direction is used to find the first-order correction to the xz-component of the stress tensor. At order $O(\epsilon)$, the momentum equation along Ox (2.16) leads to

$$ \frac{Bi}{Re} \frac{\partial \tilde{\tau}_{xz}^{(1)}}{\partial z} = \frac{\partial \tilde{u}^{(0)}}{\partial t} + \tilde{u}^{(0)} \frac{\partial \tilde{u}^{(0)}}{\partial x} + \frac{1}{Fr^2} \frac{\partial \tilde{p}^{(0)}}{\partial x} - \frac{Bi}{Re} \frac{\partial \tilde{\tau}_{xx}^{(0)}}{\partial x}. \qquad (3.24) $$

Using the expression (3.11) of the zero-order pressure, this expression can be written

$$ \frac{Bi}{Re} \frac{\partial \tilde{\tau}_{xz}^{(1)}}{\partial z} = \frac{\partial \tilde{u}^{(0)}}{\partial t} + \tilde{u}^{(0)} \frac{\partial \tilde{u}^{(0)}}{\partial x} + \frac{1}{Fr^2} \frac{\partial h}{\partial x} - 2 \frac{Bi}{Re} \frac{\partial \tilde{\tau}_{xx}^{(0)}}{\partial x}. \qquad (3.25) $$
Since $\tilde{\tau}_{xx}^{Y(0)}$ is a function of the single variable $h-z$, we have a relation between its derivatives with respect to $x$ and $z$, which is

$$ \frac{\partial \tilde{\tau}_{xx}^{Y(0)}}{\partial x} = - \frac{\partial \tilde{\tau}_{xx}^{Y(0)}}{\partial z} \frac{\partial h}{\partial x} \qquad (3.26) $$

The first dynamic boundary condition (2.24) at $O(\epsilon)$ is

$$ \tilde{\tau}_{xz}^{(1)}|_{z=h(x)} = 2\delta \frac{\partial h}{\partial x} = 2 \frac{\partial h}{\partial x} \tilde{\tau}_{xx}^{Y(0)}|_{z=h(x)}. \qquad (3.27) $$

Integration of equation (3.24) coupled with (3.27) thus leads to the following expression for the shear stress correction:

$$ \tilde{\tau}_{xz}^{(1)} = \frac{Re}{Bi}(z-h) \left[ \lambda (h-h_p) \frac{\partial h}{\partial t} + \frac{\lambda^2}{2} (h-h_p)^3 \frac{\partial h}{\partial x} \right] + \frac{Re}{Bi Fr^2} (z-h) \frac{\partial h}{\partial x} + 2 \tilde{\tau}_{xx}^{Y(0)} \frac{\partial h}{\partial x}. \qquad (3.28) $$

Using the assumptions detailed in §2.1 to determine the viscous and yield-stress parts of the stress tensor, we obtain the viscous part by taking $\tilde{\tau}_{xx}^{Y(0)} = 0$ in the above expression. This leads to

$$ \tilde{\tau}_{xz}^{Y(1)} = 2 \tilde{\tau}_{xx}^{Y(0)} \frac{\partial h}{\partial x}, \qquad (3.29) $$

while the first-order correction to the viscous contribution is

$$ \tilde{\tau}_{xz}^{v(1)} = \frac{Re}{Bi}(z-h) \left[ \lambda (h-h_p) \frac{\partial h}{\partial t} + \frac{\lambda^2}{2} (h-h_p)^3 \frac{\partial h}{\partial x} \right] + \frac{Re}{Bi Fr^2} (z-h) \frac{\partial h}{\partial x}. \qquad (3.30) $$

The expression of $\tilde{\tau}_{xx}^{Y(1)}$ could then be found from the expression of $\tilde{\tau}_{xz}^{Y(1)}$ and from the condition $|\tilde{\tau}^Y| = 1$.

The first-order correction to the velocity is obtained using the constitutive law (2.20), which becomes at order 1 $\partial \tilde{u}^{(1)}/\partial z = Bi \tilde{\tau}_{xz}^{v(1)}$. We thus obtain the following expression for the first-order correction of the velocity:

$$ \tilde{u}^{(1)} = Re \left( \frac{z^2}{2} - hz \right) \left[ \lambda (h-h_p) \frac{\partial h}{\partial t} + \frac{\lambda^2}{2} (h-h_p)^3 \frac{\partial h}{\partial x} \right] + \frac{Re}{Fr^2} \left( \frac{z^2}{2} - hz \right) \frac{\partial h}{\partial x} + u_+(x, t), \qquad (3.31) $$

where $u_+$ is a matching term to be determined later. The first term on right side of (3.31) corresponds to an inertial contribution, while the second term corresponds to a contribution due to the hydrostatic pressure.

It should be noted that, in contrast to equation (3.21) derived from the classical constitutive law (2.1), the velocity derivative $\partial \tilde{u}^{(1)}/\partial z$ remains here bounded everywhere in the pseudo-plug zone. Moreover, the first-order correction $\tilde{u}^{(1)}$ is here controlled by the viscous stress $\tilde{\tau}_{xz}^{v(1)}$ and not by the yield-stress term $\tilde{\tau}_{xz}^{Y(0)}$. As a consequence, the expression (3.31) for $\tilde{u}^{(1)}$ includes the inertial terms, which was not the case in the previous approach.

### 3.2.3. In the sheared layer ($z \le h - h_p$)

The momentum equation along $Ox$ (2.16) at $O(\epsilon)$ writes here

$$ \frac{Bi}{Re} \frac{\partial \tau_{xz}^{(1)}}{\partial z} = \frac{\partial u^{(0)}}{\partial t} + u^{(0)} \frac{\partial u^{(0)}}{\partial x} + w^{(0)} \frac{\partial u^{(0)}}{\partial z} + \frac{1}{Fr^2} \frac{\partial p^{(0)}}{\partial x} = \lambda z \frac{\partial h}{\partial x} + \frac{\lambda^2 z^2}{2} (h - h_p) \frac{\partial h}{\partial x} + \frac{1}{Fr^2} \frac{\partial h}{\partial x}. \qquad (3.32) $$

Integration of this equation yields

$$ \tau_{xz}^{(1)} = \frac{Re}{Bi} \left[ \frac{\lambda z^2}{2} \frac{\partial h}{\partial t} + \frac{\lambda^2 z^3}{6} (h - h_p) \frac{\partial h}{\partial x} \right] + z \frac{Re}{Bi Fr^2} \frac{\partial h}{\partial x} + \tau_{xz}^{(1)}|_{z=0}, \qquad (3.33) $$

where the unknown term $\tau_{xz}^{(1)}|_{z=0}$ can be found from the stress continuity condition at the fake yield surface $\tau_{xz}^{(1)}|_{z=h-h_p} = \tilde{\tau}_{xz}^{(1)}|_{z=h-h_p}$. This leads to the following expression for the shear stress correction:

$$ \tau_{xz}^{(1)} = \frac{Re\lambda}{2Bi} (z^2 - h^2 + h_p^2) \frac{\partial h}{\partial t} + \frac{Re\lambda^2}{6Bi} (h - h_p) [z^3 - (h - h_p)^2 (h + 2h_p)] \frac{\partial h}{\partial x} + \frac{Re}{Bi Fr^2} (z - h) \frac{\partial h}{\partial x}. \qquad (3.34) $$
<figure><img src="image_4.png" /><figcaption>**Figure 2:** Velocity profiles for the leading order solution (black dashed-dot curve), for the classical O(ε)-solution (dashed blue curve) and for the new O(ε)-solution (red curve), obtained with h = 1, ∂h/∂x = -0.1, Re = 1, Fr = 0.84, Bi = 0.3 and θ = 35° (λ = 1): (a) dimensionless velocity profiles; (b) velocity profiles normalized by the depth-averaged value U.</figcaption></figure>

Using the assumptions of §2.1, we obtain that the yield-stress contribution $\tau_{xz}^{Y(1)}$ is zero. Then from (2.20) and the no-slip condition (2.22), the following velocity profile is obtained:

$$ u^{(1)} = \frac{Re\lambda}{2} z \left( \frac{z^2}{3} - h^2 + h_p^2 \right) \frac{\partial h}{\partial t} + \frac{Re\lambda^2}{6} (h-h_p) z \left[ \frac{z^3}{4} - (h-h_p)^2 (h+2h_p) \right] \frac{\partial h}{\partial x} + \frac{Re}{Fr^2} z \left( \frac{z}{2} - h \right) \frac{\partial h}{\partial x}. \quad (3.35) $$

### 3.2.4. Matching

The matching term $u_+$ in (3.31) is obtained from the continuity condition $u^{(1)}|_{z=h-h_p} = \tilde{u}^{(1)}|_{z=h-h_p}$. Finally, the correction of the velocity profile in the pseudo-plug thus expresses as

$$ \tilde{u}^{(1)} = \frac{Re\lambda}{6}(h-h_p) \left[ 3z^2 - 6hz + (h-h_p)^2 \right] \frac{\partial h}{\partial t} + \frac{Re\lambda^2}{8}(h-h_p)^3 \left[ 2z^2 - 4hz + (h-h_p)^2 \right] \frac{\partial h}{\partial x} + \frac{Re}{Fr^2} z \left( \frac{z}{2} - h \right) \frac{\partial h}{\partial x}. \quad (3.36) $$

### 3.3. Velocity profile comparisons

In this section we compare the asymptotic expansion of the longitudinal velocity $u$ derived by using the classical expression of the yield-stress tensor (2.1) with the new asymptotic expansion obtained when this relation is relaxed. Figure 2a presents the profiles of the dimensionless velocity for the leading order (black dash-dotted curve), for the classical O(ε)-solution (red dashed curve) and for the new O(ε)-solution (red curve). Figure 2b shows the profiles of the velocity normalized by the depth-averaged value in each case. As already noted, at leading order all variables (longitudinal and normal velocities $u$ and $w$, pressure $p$ and stress components $\tau_{ij}$) are the same in both cases (see Figure 2a). At order O(ε), the differences appear only in the pseudo-plug zone. In our approach, shearing in the pseudo-plug zone is related to the viscous contribution given by the first-order correction to the shear stress (3.30), which naturally includes the inertial terms. On the contrary, with the classical approach, the shear rate is controlled by the yield-stress contribution to the stress. As seen in Figure 2b, the leading-order solution clearly features a plug zone. The O(ε) correction is expected to introduce a slight shearing in this zone. However, for the classical approach, the derivative (3.21) diverges at $z = h - h_p$, which corresponds to an infinite strain rate and leads to a non-physical kink in the velocity profile at the fake yield surface. In contrast, with our approach, the velocity profile shows a smooth transition from the sheared layer to the pseudo-plug zone. This property is obtained without the need to introduce a transition layer and, furthermore, all terms are given by analytical expressions.

## 4. Depth-averaged model

In this section, we derive a model by averaging the Cauchy mass and momentum equations over the fluid depth, taking into account the boundary conditions and the formal asymptotic expansions at order O(ε) obtained in the
previous section. An important benefit of this approach is that the dimensionality of the final model is reduced by one compared to the initial governing equations (2.15)-(2.17), and that the boundary conditions (2.22)-(2.25) are directly included into the equations of the model. Accordingly, this approach is expected to make numerical solutions faster and easier to compute, provided that the final model has a proper mathematical structure.

For any variable A of the flow, let us define the depth-averaged value $\langle A \rangle$ by

$$ \langle A \rangle = \frac{1}{h} \int_{0}^{h} Adz. \qquad (4.1) $$

For the averaged longitudinal velocity, it is convenient to use further the special notation $\langle u \rangle = U$. In particular, the expression for the leading order term $U^{(0)}$ is readily obtained from the equations (3.15) and (3.18):

$$ U^{(0)} = \frac{\lambda h^2}{3} \left(1 - \frac{h_p}{h}\right)^2 \left(1 + \frac{h_p}{2h}\right). \qquad (4.2) $$

## 4.1. Mass conservation

Averaging the continuity equation (2.15), taking into account the kinematic boundary condition (2.23), yields the following exact equation for mass conservation:

$$ \frac{\partial h}{\partial t} + \frac{\partial hU}{\partial x} = 0. \qquad (4.3) $$

Introducing the leading-order asymptotic solution for $U$, i.e. $U = U^{(0)}$ into (4.3) results in a kinematic wave equation:

$$ \frac{\partial h}{\partial t} + \lambda h(h - h_p) \frac{\partial h}{\partial x} = O(\epsilon), \qquad (4.4) $$

which indicate that a perturbation of depth is propagated at the speed $c_0 = \lambda h(h - h_p)$. Equation (4.4) provides a useful expansion for $\partial h/\partial t$:

$$ \frac{\partial h}{\partial t} = -\lambda h(h - h_p) \frac{\partial h}{\partial x} + O(\epsilon). \qquad (4.5) $$

In particular, this expansion allows us to write the first-order correction $U^{(1)}$ in the following form:

$$ \frac{3U^{(1)}}{Reh} = \left[ \frac{2}{5}\lambda^2 h^3 \left(1 - \frac{h_p}{h}\right)^2 \left(1 + \frac{h_p}{h} + \frac{h_p^2}{h^2} - \frac{h_p^3}{4h^3} - \frac{h_p^4}{4h^4}\right) - \frac{1}{Fr^2} \right] h \frac{\partial h}{\partial x}. \qquad (4.6) $$

## 4.2. Momentum balance

Averaging the momentum balance equation in the $Ox$-direction (2.16), together with the no-slip condition (2.22), the kinematic boundary condition (2.23), and the dynamic boundary conditions (2.24) and (2.25), provides

$$ \frac{\partial}{\partial t}(hU) + \frac{\partial}{\partial x} \left\{ h\langle u^2 \rangle + \frac{1}{Fr^2} \int_0^h pdz - \frac{Bi}{Re} \int_{h-h_p}^h \tau_{xx} dz \right\} = \frac{\lambda h - Bi\tau_{xz}(0)}{\epsilon Re}, \qquad (4.7) $$

where the quantity $\tau_{xz}(0)$ denotes the shear stress at the bottom $z = 0$. To express the integral of the pressure in (4.7) we use the equations (3.11) and (3.17):

$$ \frac{1}{Fr^2} \int_0^h pdz = \frac{h^2}{2Fr^2} - \frac{Bi}{Re} \int_{h-h_p}^h \tau_{xx}^{(0)} dz + O(\epsilon). \qquad (4.8) $$
Taking into account the expression for the normal stress at leading-order (3.10), the integral term on the right-hand side of (4.8) writes

$$ \int_{h-h_p}^{h} \tau_{xx}^{(0)} dz = \frac{\pi h_p}{4}. \qquad (4.9) $$

As a result, the integral terms involving the normal stress $\tau_{xx}^{(0)}$ are constant and disappear from equation (4.7) after differentiation. Further, using the leading-order representation for the shear stress (3.2), the averaged momentum equation (4.7) reduces to

$$ \frac{\partial}{\partial t}(hU) + \frac{\partial}{\partial x}\left(h\langle u^2 \rangle + \frac{h^2}{2Fr^2}\right) = -\frac{Bi\tau_{xz}^{(1)}(0)}{Re} + O(\epsilon^2). \qquad (4.10) $$

The quantity $\langle u^2 \rangle$ in (4.10) can be expressed by considering the velocity as the sum of its average value $U$ and a deviation $u^*$:

$$ u(x, z, t) = U(x, t) + u^*(x, z, t). \qquad (4.11) $$

By definition $\langle u^* \rangle = 0$, so that $\langle u^2 \rangle = U^2 + \langle u^{*2} \rangle$. In equation (4.10), the deviation $\langle u^{*2} \rangle$ could be estimated at order 0 either as a function of $h$ or as a function of $U$. This would lead to a closed system of two coupled equations for the flow height $h$ and the depth-averaged velocity $U$ [23, 26]. However, as shown by Richard et al. [27] for the Newtonian fluid, while for the primitive Cauchy equations the momentum balance equation and the kinetic energy equation are equivalent, this is not the case for the depth-averaged equations. A more robust approach, initially proposed by Teshukov [44] and expanded by Richard and Gavrilyuk [45], is thus to define an independent variable related to $\langle u^{*2} \rangle$. The introduction of this new variable guarantees the compatibility of the averaged mass and momentum equations with the averaged energy equation (see §4.3). We adopt here this three-variable approach, and characterize the flow based on its depth $h$, its average velocity $U$ and the variance of its velocity. In fact, it is more convenient to use

$$ \varphi = \frac{\langle u^{*2} \rangle}{h^2} \qquad (4.12) $$

as the third variable of the model, since it plays the role of an entropy for the system [45]. In the particular case of a constant vorticity, the quantity $\varphi$ is proportional to the square of the vorticity [44]. For this reason, $\varphi$ is called entrophy. The expansion of this new variable

$$ \varphi = \varphi^{(0)} + \varepsilon\varphi^{(1)} + O(\varepsilon^2) \qquad (4.13) $$

leads to $\varphi^{(0)} = \langle (u^{*(0)})^2 \rangle / h^2$ and $\varphi^{(1)} = 2\langle u^{*(0)}u^{*(1)} \rangle / h^2$. The calculation at leading order gives

$$ \varphi^{(0)} = \frac{\lambda^2 h^2}{45} \left(1 - \frac{h_p}{h}\right)^5 \left(1 + \frac{5h_p}{4h}\right). \qquad (4.14) $$

For the expansion at order $O(\varepsilon)$, we obtain

$$ \frac{\varphi^{(1)}}{Re} = \frac{2\lambda h^2}{45} \left[ \frac{3}{7}\lambda^2 h^3 \left(1 - \frac{h_p}{h}\right)^2 \left(1 + \frac{101h_p}{48h} + \frac{29h_p^2}{12h^2} + \frac{43h_p^3}{48h^3} - \frac{7h_p^4}{12h^4}\right) \right. \\ \left. - \frac{1}{Fr^2} \left(1 + \frac{9h_p}{8h} + \frac{3h_p^2}{8h^2}\right) \right] \left(1 - \frac{h_p}{h}\right)^3 \frac{\partial h}{\partial x}. \qquad (4.15) $$

Finally, with the introduction of the entrophy $\varphi$, the equation (4.10) can be written

$$ \frac{\partial}{\partial t}(hU) + \frac{\partial}{\partial x}\left(hU^2 + h^3\varphi + \frac{h^2}{2Fr^2}\right) = -\frac{Bi\tau_{xz}^{(1)}(0)}{Re}. \qquad (4.16) $$
To keep a proper mathematical structure, the main idea now is to remove all derivatives from the right-hand side and to express the source terms of equation (4.16) as a sum of relaxation terms. Using the relation (4.5), the stress at the bottom $\tau_{xz}^{(1)}(0)$ given by (3.34) reduces to

$$ \frac{Bi\tau_{xz}^{(1)}}{Re} = \left[ \frac{1}{3}\lambda^2 h^3 \left(1 - \frac{h_p}{h}\right)^2 \left(1 + \frac{h_p}{h} + \frac{h_p^2}{h^2}\right) - \frac{1}{Fr^2} \right] h \frac{\partial h}{\partial x}, \quad (4.17) $$

which can also be rewritten in terms of $U^{(1)}$:

$$ \frac{Bi\tau_{xz}^{(1)}}{Re} = \frac{3U^{(1)}}{Reh} - \frac{\lambda^2 h^4}{15} \left(1 - \frac{h_p}{h}\right)^3 \left(1 + \frac{2h_p}{h} + \frac{3h_p^2}{h^2} + \frac{3h_p^3}{2h^3}\right) \frac{\partial h}{\partial x}. \quad (4.18) $$

Note that the quantity $U^{(1)}$ can be expressed as:

$$ U^{(1)} = \frac{U - U^{(0)}}{\epsilon} + O(\epsilon) \quad (4.19) $$

which has the structure of a relaxation term in $U$. To express the derivative $\partial h / \partial x$ in (4.18), let us consider the following relation obtained from the expansions for $U$ and $\varphi$:

$$ \begin{aligned} \varphi - \frac{U^2}{5h^2} \left(1 - \frac{h_p}{h}\right) \left(1 + \frac{5h_p}{4h}\right) \left(1 + \frac{h_p}{2h}\right)^{-2} & \\ = \frac{2\epsilon}{35} Re \lambda \varphi h^3 \left(1 - \frac{h_p}{h}\right)^2 \left(1 + \frac{5h_p}{4h}\right)^{-1} & \left(1 + \frac{61h_p}{16h} + \frac{63h_p^2}{8h^2} + \frac{63h_p^3}{8h^3} + \frac{21h_p^4}{16h^4}\right) \frac{\partial h}{\partial x} \\ & + \frac{\lambda h_p}{20h} \left(1 - \frac{h_p}{h}\right)^3 \left(1 + \frac{5h_p}{2h} + \frac{h_p^2}{2h^2}\right) \left(1 + \frac{h_p}{2h}\right)^{-1} \left[U - U^{(0)}\right] + O(\epsilon^2). \end{aligned} \quad (4.20) $$

The left-hand side of (4.20) has the structure of a relaxation term in $\varphi$. Indeed, using expressions (4.2) and (4.14) we obtain

$$ \varphi^{(0)} - \frac{(U^{(0)})^2}{5h^2} \left(1 - \frac{h_p}{h}\right) \left(1 + \frac{5h_p}{4h}\right) \left(1 + \frac{h_p}{2h}\right)^{-2} = 0. \quad (4.21) $$

As a result, the derivative $\partial h / \partial x$ can be expressed as a sum of two relaxation terms for $U$ and $\varphi$. Finally, the integrated momentum equation (4.16) can thus be written as

$$ \begin{aligned} \frac{\partial}{\partial t}(hU) + \frac{\partial}{\partial x}\left(hU^2 + h^3\varphi + \frac{h^2}{2Fr^2}\right) = \frac{1}{\epsilon Re}\left[\lambda h - Bi - \frac{3U}{h\alpha_1(\xi)}\right]\left[\alpha_1(\xi) + \frac{7}{360}\frac{\lambda^2 h^2}{\varphi}\beta_1(\xi)\right] \\ + \frac{1}{\epsilon Re}\frac{7}{6}\frac{\lambda h}{\varphi}\left[\varphi - \frac{U^2}{5h^2}\alpha_2(\xi)\right]\beta_2(\xi) \end{aligned} \quad (4.22) $$

where $\xi = h_p / h$ and the functions $\alpha_1, \alpha_2, \beta_1$ and $\beta_2$ are defined as:

$$ \alpha_1(\xi) = (1 - \xi) \left(1 + \frac{\xi}{2}\right) \quad (4.23) $$

$$ \alpha_2(\xi) = (1 - \xi) \left(1 + \frac{5}{4}\xi\right) \left(1 + \frac{\xi}{2}\right)^{-2} \quad (4.24) $$

$$ \beta_1(\xi) = \xi(1-\xi)^5 \left(1 + \frac{5}{4}\xi\right) \frac{\left(1 + \frac{5}{2}\xi + \frac{1}{2}\xi^2\right) \left(1 + 2\xi + 3\xi^2 + \frac{3}{2}\xi^3\right)}{1 + \frac{61}{16}\xi + \frac{63}{8}\xi^2 + \frac{63}{8}\xi^3 + \frac{21}{16}\xi^4} \quad (4.25) $$
$$ \beta_2(\xi) = (1 - \xi) \left(1 + \frac{5}{4}\xi\right) \frac{1 + 2\xi + 3\xi^2 + \frac{3}{2}\xi^3}{1 + \frac{61}{16}\xi + \frac{63}{8}\xi^2 + \frac{63}{8}\xi^3 + \frac{21}{16}\xi^4} \quad (4.26) $$

The first relaxation term on right-hand side of (4.22) interprets as the balance between the gravity, the yield stress and the viscous friction forces along the Ox-axis. The second term corresponds to a relaxation for the entrophy.

### 4.3. Kinetic energy equation

The proposed depth-averaged model involves three unknown variables, namely h, U and φ. To close the problem, a third equation is thus needed, which is provided by the energy conservation. In dimensional form, the kinetic energy equation (or work-energy theorem) can be written as

$$ \frac{\partial}{\partial t} \left( \frac{1}{2} \rho v^2 \right) + \operatorname{div} \left( \frac{1}{2} \rho v^2 \boldsymbol{v} \right) = \operatorname{div}(\boldsymbol{\sigma} \cdot \boldsymbol{v}) - \boldsymbol{\sigma} : \frac{\dot{\boldsymbol{v}}}{2} + \rho g \boldsymbol{v} \quad (4.27) $$

where we recall that v denotes the velocity field and $\boldsymbol{\sigma} = -p\boldsymbol{I} + \boldsymbol{\tau}$. Introducing the components of the vectors and tensors involved in (4.27), we get the following expression written in dimensionless form:

$$ \begin{aligned} & \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 F r^2} + \frac{z}{Fr^2} \right) + \frac{p u}{Fr^2} - \frac{B i}{R e} (\tau_{xx} u + \varepsilon \tau_{xz} w) \right] \\ & + \frac{\partial}{\partial z} \left[ w \left( \frac{u^2}{2} + \varepsilon^2 \frac{w^2}{2} - \frac{x \tan \theta}{\varepsilon F r^2} + \frac{z}{Fr^2} \right) + \frac{p w}{Fr^2} - \frac{B i}{\varepsilon R e} \tau_{xz} u - \frac{B i}{R e} \tau_{zz} w \right] \\ & = -\frac{2 B i}{R e} \left( \tau_{xx}^Y + 2 \frac{\varepsilon}{B i} \frac{\partial u}{\partial x} \right) \frac{\partial u}{\partial x} - \frac{B i}{\varepsilon R e} \left( \tau_{xz}^Y + \frac{1}{B i} \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{aligned} \quad (4.28) $$

This equation has to be averaged over the depth. The details of this calculation are presented in appendix A. Taking into account the boundary conditions and dropping all second-order terms, the depth-averaged work-energy theorem can finally be expressed as:

$$ \frac{\partial}{\partial t} \left( \frac{hU^2}{2} + \frac{h^3\varphi}{2} + \frac{h^2}{2Fr^2} \right) + \frac{\partial}{\partial x} \left( \frac{h\langle u^3 \rangle}{2} + \frac{h^2U}{Fr^2} \right) = -\frac{\lambda h U^{(1)}}{Re} + \frac{Bi}{Re} \int_0^h \tau_{xz}^{Y(0)} \frac{\partial u^{(1)}}{\partial z} dz. \quad (4.29) $$

From the equation (4.11) and the definition of entrophy φ, we have $\langle u^2 \rangle = U^2 + h^2\varphi$ and

$$ \langle u^3 \rangle = U^3 + 3h^2 U \varphi + \langle u^{*3} \rangle. \quad (4.30) $$

Accordingly, (4.29) can be rewritten as:

$$ \begin{aligned} & \frac{\partial}{\partial t} \left( \frac{hU^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) = \\ & = -\frac{\lambda h U^{(1)}}{Re} + \frac{Bi}{Re} \int_0^h \tau_{xz}^{Y(0)} \frac{\partial u^{(1)}}{\partial z} dz - \frac{\partial}{\partial x} \left( \frac{h\langle u^{*3} \rangle}{2} \right). \end{aligned} \quad (4.31) $$

As for the momentum equation, we wish to express the right-hand side of (4.31) under the form of relaxation terms. For turbulent water flows, Richard and Gavrilyuk [45] considered a weakly-sheared flow assumption, which allowed them to neglect the term involving the cubic deviation $\langle u^{*3} \rangle$ and thus to close the model. In the present work, since the solution must satisfy the no-slip condition, this assumption is not appropriate. On the other hand, treating $\langle u^{*3} \rangle$ as a fourth variable would lead to an infinite hierarchy of equations. To close the system, we thus calculate $\langle u^{*3} \rangle$ as a function of h by employing the asymptotic expansions derived above. Namely, using expression (4.30), we obtain the following expression at leading order:

$$ \langle u^{*3} \rangle = -\frac{2}{945} \lambda^3 h^6 \left(1 - \frac{h_p}{h}\right)^7 \left(1 + \frac{49h_p}{16h} + \frac{35h_p^2}{8h^2}\right) + O(\varepsilon). \quad (4.32) $$
The derivative of $\langle u^{*3} \rangle$ involved in (4.31) has thus the form:

$$ -\frac{\partial}{\partial x} \left( \frac{h \langle u^{*3} \rangle}{2} \right) = \frac{\lambda^3 h^6}{135} \left( 1 - \frac{h_p}{h} \right)^6 \left( 1 + \frac{21 h_p}{8 h} + \frac{57 h_p^2}{16 h^2} + \frac{5 h_p^3}{4 h^3} \right) \frac{\partial h}{\partial x} + O(\epsilon). \quad (4.33) $$

The remaining term to be calculated in (4.31) is the integral:

$$ \frac{Bi}{Re} \int_{0}^{h} \tau_{xz}^{Y(0)} \frac{\partial u^{(1)}}{\partial z} dz = \frac{\lambda h_p}{2} \left[ \frac{5}{12} \lambda^2 h^5 \left(1 - \frac{h_p}{h}\right)^2 \left(1 + \frac{h_p}{h} + \frac{3h_p^2}{5h^2} - \frac{h_p^3}{h^3}\right) - h^2 Fr^{-2} \left(1 - \frac{h_p^2}{3h^2}\right) \right] \frac{\partial h}{\partial x} \quad (4.34) $$

In the above expression, the term involving Froude number can be expressed in terms of $U^{(1)}$ through equation (4.6). Finally, employing also expression (4.19), the averaged energy equation (4.31) can be rewritten as follows:

$$ \frac{\partial}{\partial t} \left( \frac{hU^2}{2} + \frac{h^3 \varphi}{2} + \frac{h^2}{2Fr^2} \right) + \frac{\partial}{\partial x} \left( \frac{hU^3}{2} + \frac{3h^3 U \varphi}{2} + \frac{h^2 U}{Fr^2} \right) = \frac{U}{\epsilon Re} \left[ \lambda h - Bi - \frac{3U}{h\alpha_1(\xi)} \right] \left[ \alpha_1(\xi) + \frac{7}{1080} \frac{\lambda^2 h^2}{\varphi} \beta_1(\xi) r(\xi) \right] + \frac{U}{\epsilon Re} \frac{7}{18} \frac{\lambda h}{\varphi} \left[ \varphi - \frac{U^2}{5h^2} \alpha_2(\xi) \right] \beta_2(\xi) r(\xi) \quad (4.35) $$

where again $\xi = h_p/h$ and the function $r$ is defined by:

$$ r(\xi) = \frac{1 + \frac{11}{4}\xi + \frac{87}{16}\xi^2 + \frac{107}{16}\xi^3 + \xi^4}{\left(1 + \frac{\xi}{2}\right)\left(1 + 2\xi + 3\xi^2 + \frac{3}{2}\xi^3\right)}. \quad (4.36) $$

The first relaxation term on right-hand side of (4.35) interprets as the balance between the power of the component of the weight along the $Ox$ axis and the power of the yield stress and the viscous friction forces. The second term corresponds to the relaxation for the entrophy $\varphi$.

### 4.4. Equation of entrophy

From equations (4.3), (4.22) and (4.35), the following evolution equation for the entrophy $\varphi$ can be derived:

$$ \frac{h^2}{2} \left( \frac{\partial h \varphi}{\partial t} + \frac{\partial h U \varphi}{\partial x} \right) = \frac{U}{\epsilon Re} \left[ \lambda h - Bi - \frac{3U}{h \alpha_1(\xi)} \right] \left[ \frac{7}{360} \frac{\lambda^2 h^2}{\varphi} \beta_1(\xi) \left( \frac{r(\xi)}{3} - 1 \right) \right] + \frac{U}{\epsilon Re} \frac{7}{6} \frac{\lambda h}{\varphi} \left[ \varphi - \frac{U^2}{5h^2} \alpha_2(\xi) \right] \beta_2(\xi) \left( \frac{r(\xi)}{3} - 1 \right). \quad (4.37) $$

Note that, since $r(\xi)/3 - 1 < 0$, the entrophy properly relaxes toward its equilibrium value $\alpha_2(\xi)U^2/5h^2$.

Formally, any set of three equations from (4.3), (4.22), (4.35) and (4.37) can be used to describe the motion of the fluid. However, to model shock waves (see Section 6), it is more natural to use the mass, momentum and energy conservation equations, while the entrophy, which plays the role of an entropy of the system, should increase according to the Rankine–Hugoniot relations (for more details see [45]). Hence, we shall only consider the system consisting of equations (4.3), (4.22), (4.35) in what follows.

### 4.5. Structure of the model

Here we show that the system of equations (4.3), (4.22) and (4.35) is equivalent to Euler equations for compressible fluids with relaxation terms. This property ensures that the model is fully hyperbolic and can be handled by efficient numerical schemes. Returning to dimensional variables, let us denote

$$ \Pi = h^3 \varphi + \frac{gh^2 \cos \theta}{2}, \quad (4.38) $$

$$ e = \frac{1}{2} (U^2 + h^2 \varphi + gh \cos \theta). \quad (4.39) $$
Equations (4.3), (4.22) and (4.35) can then be written as (in dimensional form):

$$ \frac{\partial h}{\partial t} + \frac{\partial hU}{\partial x} = 0, \qquad (4.40) $$

$$ \frac{\partial}{\partial t}(hU) + \frac{\partial}{\partial x}(hU^2 + \Pi) = \left[gh \sin \theta - \frac{\tau_c}{\rho} - 3\nu \frac{U}{h\alpha_1(\xi)}\right] \left[\alpha_1(\xi) + \frac{7}{360} \frac{(gh \sin \theta)^2}{\nu^2 \varphi} \beta_1(\xi)\right] + \frac{7}{6} \frac{gh \sin \theta}{\varphi} \left[\varphi - \frac{U^2}{5h^2} \alpha_2(\xi)\right] \beta_2(\xi), \quad (4.41) $$

$$ \frac{\partial}{\partial t}(he) + \frac{\partial}{\partial x}\left[hU\left(e + \frac{\Pi}{h}\right)\right] = U\left[gh \sin \theta - \frac{\tau_c}{\rho} - 3\nu \frac{U}{h\alpha_1(\xi)}\right] \left[\alpha_1(\xi) + \frac{7}{1080} \frac{(gh \sin \theta)^2}{\nu^2 \varphi} \beta_1(\xi)r(\xi)\right] + \frac{7}{18} \frac{Ugh \sin \theta}{\varphi} \left[\varphi - \frac{U^2}{5h^2} \alpha_2(\xi)\right] \beta_2(\xi)r(\xi) \quad (4.42) $$

with $v = K/\rho$, $\xi = \tau_c/(\rho gh \sin \theta)$, and the functions $r$, $\alpha_i$, $\beta_i$ ($i = 1, 2$) defined in (4.23)-(4.26) and (4.36). Hence, the left-hand side of system (4.40)-(4.42) indeed has the form of the Euler equations for compressible fluids with $\Pi$, $e$ and $\varphi$ playing the roles of pressure, energy and entropy, respectively.

The reduced system obtained for the Newtonian case ($\tau_c = 0$) is given in Appendix B. It should be noted that, in this Newtonian case, the present model improves on the former model derived by Richard et al. [28] in that the enstrophy relaxes toward a more physical value $U^2/5h^2$ (in particular, this implies that $\varphi = 0$ if $U = 0$).

## 5. Reconstruction of the velocity field

The solution of the system (4.40)-(4.42) provides us with the free-surface $h$, the averaged velocity $U$ and the enstrophy $\varphi$ at order $O(\varepsilon)$. For a comprehensive analysis of the flows, it can also be useful to reconstruct the velocity field $(u, w)$ from these computed variables $h$, $U$ and $\varphi$. Let us recall that, unlike the previous studies based on the classical expression (2.1) for the yield-stress tensor, a notable benefit of our model is the possibility to consistently reconstruct velocity profiles that are smooth at order $O(\varepsilon)$ (see Section 3.3).

At leading order, the longitudinal velocity $u$ can be expressed as

$$ u^{(0)} = U^{(0)} f_{sh}(z) \quad \text{for } z < h - h_p, \qquad (5.1) $$

$$ \tilde{u}^{(0)} = U^{(0)} f_{pl} \quad \text{for } z \ge h - h_p, \qquad (5.2) $$

with

$$ f_{sh}(z) = \frac{3z}{h} \left(1 - \frac{z}{2(h - h_p)}\right) \left(1 - \frac{h_p}{h}\right)^{-1} \left(1 + \frac{h_p}{2h}\right)^{-1}, \quad f_{pl} = \frac{3}{2} \left(1 + \frac{h_p}{2h}\right)^{-1}. \qquad (5.3) $$

Note that since the quantity $U^{(0)}$ is the leading-order solution for the average velocity, we have

$$ \frac{1}{h} \int_{0}^{h-h_p} f_{sh}(z) \, dz + \frac{1}{h} \int_{h-h_p}^{h} f_{pl} \, dz = 1. \qquad (5.4) $$

At order $O(\varepsilon)$, the expansions for the longitudinal velocity are written

$$ u = U f_{sh} + \varepsilon [u^{(1)} - U^{(1)} f_{sh}] + O(\varepsilon^2) \quad \text{for } z < h - h_p, \qquad (5.5) $$

$$ u = U f_{pl} + \varepsilon [\tilde{u}^{(1)} - U^{(1)} f_{pl}] + O(\varepsilon^2) \quad \text{for } z \ge h - h_p, \qquad (5.6) $$
where $U f_{sh}$ and $U f_{pl}$ give the velocity profile at order 0 and the second terms in the right-hand side of these expressions give the correction of order 1. These expressions can be used to reconstruct the velocity profile with the values of $h$, $U$ and $\varphi$, excluding their derivatives, since the derivatives $\partial h/\partial t$ and $\partial h/\partial x$, which are involved in the first-order corrections $u^{(1)}$ and $\tilde{u}^{(1)}$ (see equations (3.35) and (3.36)), can be consistently expressed in terms of $h$, $U$ and $\varphi$ using (4.5), (4.14), (4.19) and (4.20). Finally, this leads to

$$ u = U f(\eta; \xi) + \left[ \chi(\eta; \xi) \frac{h^2}{2Fr^2} + \psi(\eta; \xi) \frac{\lambda^2 h^5}{15} \right] \times \left( \frac{7}{24\varphi h^2} \left[ \lambda h - Bi - \frac{3U}{h\alpha_1(\xi)} \right] \zeta_1(\xi) + \frac{7}{18}\frac{\lambda}{\varphi^2 h} \left[ \varphi - \frac{U^2}{5h^2}\alpha_2(\xi) \right] \zeta_2(\xi) \right) \quad (5.7) $$

with $\xi = h_p/h$ and $\eta = z/h$. As above, the function $f$, which corresponds to $f_{sh}$ and $f_{pl}$, gives the profile of order 0 and the other terms give the corrections of order 1. For the sheared layer, $z < h - h_p$, the functions $f$, $\chi$ and $\psi$ are defined by

$$ f(\eta; \xi) = \frac{3\eta}{(1-\xi)^2} \left(1 - \xi - \frac{\eta}{2}\right) \left(1 + \frac{1}{2}\xi\right)^{-1}, \quad (5.8) $$

$$ \chi(\eta; \xi) = \frac{\xi\eta}{(1-\xi)^2} \left(1 - \frac{3}{2}\eta - \xi^2 + \frac{1}{2}\eta\xi^2\right) \left(1 + \frac{1}{2}\xi\right)^{-1}, \quad (5.9) $$

$$ \psi(\eta; \xi) = -\eta \left( 1 + \frac{5}{2}\xi - 3\eta + \frac{5}{2}\xi^2 - 3\xi\eta + \frac{5}{2}\eta^2 - \frac{5}{2}\xi^3 - 3\xi^2\eta - \frac{5}{4}\xi\eta^2 - \frac{5}{8}\eta^3 \right. \\ \left. - \frac{5}{2}\xi^4 + \frac{3}{4}\xi^3\eta - \frac{5}{4}\xi^2\eta^2 + \frac{5}{16}\xi\eta^3 - \xi^5 + \frac{3}{4}\xi^4\eta + \frac{5}{16}\xi^2\eta^3 \right) \left(1 + \frac{1}{2}\xi\right)^{-1} \quad (5.10) $$

while for the pseudo-plug zone, $z \ge h - h_p$, they are defined as

$$ f(\eta; \xi) = \frac{3}{2} \left(1 + \frac{1}{2}\xi\right)^{-1}, \quad (5.11) $$

$$ \chi(\eta; \xi) = \left(1 - 2\eta + \eta^2 - \xi\eta + \frac{1}{2}\xi\eta^2\right) \left(1 + \frac{1}{2}\xi\right)^{-1}, \quad (5.12) $$

$$ \psi(\eta; \xi) = -\frac{29}{8} \left( 1 + \frac{63}{58}\xi - \frac{60}{29}\eta + \frac{3}{58}\xi^2 - \frac{90}{29}\xi\eta + \frac{30}{29}\eta^2 - \frac{30}{29}\xi^2\eta \right. \\ \left. - \frac{7}{58}\xi^3 + \frac{45}{29}\xi\eta^2 + \frac{3}{58}\xi^4 + \frac{15}{29}\xi^2\eta^2 \right) \left(1 - \xi\right)^2 \left(1 + \frac{1}{2}\xi\right)^{-1}. \quad (5.13) $$

The functions $\zeta_1$ and $\zeta_2$ are given by

$$ \zeta_1(\xi) = \xi \left(1 - \xi\right)^2 \left(1 + \frac{5}{4}\xi\right) \left(1 + \frac{5}{2}\xi + \frac{1}{2}\xi^2\right) \left(1 + \frac{61}{16}\xi + \frac{63}{8}\xi^2 + \frac{63}{8}\xi^3 + \frac{21}{16}\xi^4\right)^{-1}, \quad (5.14) $$

$$ \zeta_2(\xi) = \left(1 - \xi\right)^3 \left(1 + \frac{5}{4}\xi\right)^2 \left(1 + \frac{61}{16}\xi + \frac{63}{8}\xi^2 + \frac{63}{8}\xi^3 + \frac{21}{16}\xi^4\right)^{-1}. \quad (5.15) $$

Lastly, the functions $\alpha_1$ and $\alpha_2$ are defined by (4.23)–(4.24).

Regarding the normal velocity $w$, the derivative $\partial h/\partial x$ in (3.16) can be expressed consistently in a similar way. Since $w$ is a first-order quantity compared to $u$, it is sufficient to express $w^{(0)}$ as a sum of relaxation terms to obtain an accuracy at order 1. In the end, the following expression is obtained:

$$ w = \frac{\lambda h^2}{\varepsilon Re} \left( \frac{7}{24\varphi h^2} \left[ \lambda h - Bi - \frac{3U}{h\alpha_1(\xi)} \right] \zeta_1(\xi) + \frac{7}{18}\frac{\lambda}{\varphi^2 h} \left[ \varphi - \frac{U^2}{5h^2}\alpha_2(\xi) \right] \zeta_2(\xi) \right) q(\eta; \xi) \quad (5.16) $$

The function $q$ for the sheared layer, $z < h - h_p$, is given by

$$ q(\eta; \xi) = -\frac{\eta^2}{2}, \quad (5.17) $$
and for the pseudo-plug, $z \ge h - h_p$, by

$$q(\eta; \xi) = -(1 - \xi) \left( \eta - \frac{1 - \xi}{2} \right). \qquad (5.18)$$

In dimensional variables the reconstructed longitudinal velocity expresses as

$$u = U f(\eta; \xi) + \left[ \chi(\eta; \xi) \frac{gh^2 \cos \theta}{2} + \psi(\eta; \xi) \frac{g^2 h^5 \sin^2 \theta}{15\nu^2} \right] \times \left( \frac{7}{24\nu\varphi h^2} \left[ gh \sin \theta - \frac{\tau_c}{\rho} - \frac{3\nu U}{h\alpha_1(\xi)} \right] \zeta_1(\xi) + \frac{7}{18} \frac{g \sin \theta}{\nu \varphi^2 h} \left[ \varphi - \frac{U^2}{5h^2} \alpha_2(\xi) \right] \zeta_2(\xi) \right) \quad (5.19)$$

with $\xi = \tau_c/(\rho g h \sin \theta)$. The dimensional normal velocity expresses as

$$w = \frac{gh^2 \sin \theta}{\nu} \left( \frac{7}{24\varphi h^2} \left[ gh \sin \theta - \frac{\tau_c}{\rho} - \frac{3\nu U}{h\alpha_1(\xi)} \right] \zeta_1(\xi) + \frac{7}{18} \frac{g \sin \theta}{\nu \varphi^2 h} \left[ \varphi - \frac{U^2}{5h^2} \alpha_2(\xi) \right] \zeta_2(\xi) \right) q(z/h; \xi). \quad (5.20)$$

Using this reconstruction of the velocity field, the strain rate $|\dot{\gamma}|$ can also be reconstructed within the fluid. If a critical value $\dot{\gamma}_c$ is chosen (see §3.1.1), an effective pseudo-plug layer can be found with the condition $|\dot{\gamma}| < \dot{\gamma}_c$. As explained in §3.1.1, this effective pseudo-plug layer differs from the formal pseudo-plug layer used for the asymptotic expansions. It includes all corrections of order 1, notably the effect of the pressure, and its thickness is generally not constant in the fluid.

# 6. Applications

## 6.1. Characteristic velocities

The left-hand side of the system (4.40)–(4.42) is identical as in Richard and Gavrilyuk [45]. The right-hand side involves only relaxation terms. Accordingly, the model is hyperbolic [44]. The system can be rewritten in the matrix form:

$$\frac{\partial \mathbf{V}}{\partial t} + \mathbf{A} \frac{\partial \mathbf{V}}{\partial x} = \mathbf{S}, \qquad (6.1)$$

where $\mathbf{V} = [h, U, \varphi]^T$, $\mathbf{S}$ is the matrix of the source terms and the matrix $\mathbf{A}$ is given by

$$\mathbf{A} = \begin{bmatrix} U & h & 0 \\ g \cos \theta + 3h\varphi & U & h^2 \\ 0 & 0 & U \end{bmatrix} \qquad (6.2)$$

The three characteristic velocities of the model are given by the eigenvalues of $\mathbf{A}$, namely:

$$U, \quad U - \sqrt{gh \cos \theta + 3h^2 \varphi}, \quad U + \sqrt{gh \cos \theta + 3h^2 \varphi} \qquad (6.3)$$

As it is clear from these expressions, the shearing effect contributes to the characteristic velocities through the term $3h^2\varphi$.

## 6.2. Long-wave instability

As explained in the previous paragraph, the characteristic velocities of the system depend only on the conservative part of the equations, and not on the relaxation terms. On the contrary, the linear instability threshold and the expression of the phase velocity of perturbations depend strongly on the relaxation source terms and, in particular, on the yield-stress and viscous friction terms.

To establish the dispersion relation, let us linearize the system of equations (4.3), (4.22) and (4.35) around the base solution (4.2) and (4.14). We write $h = 1 + h'$, $U = U^{(0)} + U'$ and $\varphi = \varphi^{(0)} + \varphi'$, where $h'$, $U'$ and $\varphi'$ are small sinusoidal perturbations. Namely, we take the perturbations of the form $[h', U', \varphi']^T = [A_1, A_2, A_3]^T \exp[ik(x - ct)]$, where $k$ is the wavenumber and $c$ is the phase velocity. The dispersion relation is found by equating the determinant
<figure><img src="image_5.png" /><figcaption>**Figure 3:** Instability threshold obtained by Balmforth and Liu [12] (black dashed curve), by the present model (red curve), and by Fernández-Nieto et al. [26] (blue curves) for different slope angles (30°, 55°, 60° and 70°): critical Reynolds number $\tan\theta Re_c$ as a function of the Bingham number $Bi$.</figcaption></figure>

of the linearised system to zero. The details of this lengthy calculation are given in Appendix C. Up to the first order in $\epsilon$, the relation is written

$$c = \lambda(1 - h_p) + ik\epsilon \frac{Re}{3} \left[ \frac{2\lambda^2}{5}(1 - h_p)^2 \left( 1 + h_p + h_p^2 - \frac{h_p^3}{4} - \frac{h_p^4}{4} \right) - \frac{1}{Fr^2} \right] + O(\epsilon^2). \quad (6.4)$$

The base flow is stable if $Im(c) < 0$.

Balmforth and Liu [12] studied the linear stability of the base flow for the linearized Cauchy equations (or generalized Orr-Sommerfeld equations) in the case of a Herschel-Bulkley fluid with power flow index $n$. In order to compare our result with the result obtained by these authors in the particular case of the Bingham fluid ($n = 1$), we choose the same characteristic velocity, namely $u_0 = gh_0^2 \sin\theta/v$. Note that this choice imposes $\lambda = 1$. As a result, in the long-wave limit ($\epsilon \to 0$) stability occurs for $Re < Re_c$, where the critical Reynolds number $Re_c$ is given by

$$Re_c = \frac{5 \cot \theta}{2(1 - Bi)^2} \frac{1}{1 + Bi + Bi^2 - Bi^3/4 - Bi^4/4} \quad (6.5)$$

Analysis of the expression (6.5) shows that the critical Reynolds number increases as the Bingham number increases. This reveals the stabilizing effect of the plasticity, which was also highlighted in former studies [12, 26]. Figure 3 compares our instability criterion (6.5) (red curve) with the one obtained by Balmforth and Liu [12] (black dashed curve) from the generalized Orr-Sommerfeld equations in the particular case of the Bingham fluid. Using our notations, this latter criterion can be written

$$Re_c = \frac{5 \cot \theta}{2(1 - Bi)^2} \frac{1 + Bi + Bi^2}{1 + 2Bi + 3Bi^2 + 3Bi^3/2} \quad (6.6)$$

As shown, the two model predicts almost identical instability thresholds. The slight difference (almost not visible in the figure) can be attributed to the fact that the slight shear in the pseudo-plug was neglected by Balmforth and Liu [12]. In contrast, in our approach the shear in the pseudo-plug is consistently taken into account. However, since the corresponding strain rate remains very small close to the equilibrium flow, its effect on the instability threshold is also very small. Note that the functions $(1 + Bi + Bi^2 - Bi^3/4 - Bi^4/4)^{-1}$ and $(1 + Bi + Bi^2)/(1 + 2Bi + 3Bi^2 + 3Bi^3/2)$ are equal for $Bi = 0$ and $Bi = 1$, and give almost the same values in the range $0 < Bi < 1$.

Figure 3 also presents the instability criterion obtained by Fernández-Nieto et al. [26] (blue curves) using a consistent shallow-flow model based on the classical constitutive law (2.1). In this case, the critical Reynolds number can be written, using our notations,

$$Re_c = \frac{5}{2(1 - Bi)^2} \frac{(1 + Bi + Bi^2) \cot \theta - 3\pi Bi^2/4}{1 + 2Bi + 3Bi^2 + 3Bi^3/2} \quad (6.7)$$
Note that, unlike the criteria (6.6) and (6.5), the product $Re_c \tan \theta$ depends on the slope angle $\theta$ in the model of Fernández-Nieto et al. [26]. More generally, large differences can be observed between this third criterion and the two former. In particular, the critical Reynolds number predicted by (6.7) becomes negative if the Bingham number is larger than

$$Bi_0 = \frac{2 \left[ 1 + \sqrt{3(\pi \tan \theta - 1)} \right]}{3\pi \tan \theta - 4} \qquad (6.8)$$

This value $Bi_0$ is smaller than 1 if $\tan \theta > 4/\pi$ ($\theta$ larger than 51.9°, approximately), and the limit of $Re_c$ for $Bi \to 1$ is $-\infty$ in this case. Hence, this criterion predicts that the flow is always unstable for a sufficiently large Bingham number if the slope is larger than $4/\pi$, meaning that plasticity has a destabilizing effect for large slopes. In contrast, the limit of the critical Reynolds number for $Bi \to 1$ is always $+\infty$ in our approach and that of Balmforth and Liu [12] (the case $Bi = 1$ corresponds to a fully solid behavior in this scaling).

Note that the scaling used in our study is the same as in [26] (in particular, the slope is assumed $O(1)$). Hence, we interpret these differences between our approach and that of [26] as the consequence of expressing the shear in the pseudo-plug based on the classical expression $\tau_c \dot{\gamma} / |\dot{\gamma}|$ for the yield-stress tensor (in [26] the shear in the pseudo-plug is given solely by the plastic term). The destabilizing effect of plasticity predicted by the criterion (6.7) for large slope angles does not seem to be physical. Furthermore, the recent experimental study on the instability of viscoplastic fluids performed by Noma et al. [46] shows good agreement with the criterion obtained by Balmforth and Liu [12] in the general case of a Herschel-Bulkley fluid, even if the shear in the pseudo-plug is neglected in this criterion. In the case of a Bingham fluid, our work, based on the alternative expression of the constitutive law, confirms that the slight strain rate in the pseudo-plug has a very small effect on the instability threshold, unlike in the approach of [26]. It can thus be concluded that the classical form $\tau_c \dot{\gamma} / |\dot{\gamma}|$ of the yield-stress tensor leads to unphysical predictions for the linear instability threshold, which constitutes an additional and strong argument in favor of the specific regularization of the constitutive law proposed in our study (§2.1).

## 6.3. Simulation of roll waves

As the system of equations (4.40)-(4.42) is hyperbolic, it can be solved by simple and robust classical numerical schemes. In this work, we use a Godunov-type scheme with a HLLC Riemann solver (for more details see Toro [47]). An initially uniform flow is perturbed by applying a small sinusoidal disturbance of fixed frequency at the entrance of the system. At the entrance, the values $h$, $U$ and $\varphi$ are imposed, while at the outlet a Neumann boundary condition is considered. Parameters are chosen so that $Re > Re_c$, allowing an instability to develop. Namely, we choose the following values for the dimensional parameters: $K = 20$ Pa s, $\rho = 1000$ kg m$^{-3}$, $h_0 = 0.01$ m, $\theta = 18^\circ$ and $u_0 = gh_0^2 \sin \theta / \nu = 1.5$ m s$^{-1}$. This corresponds to the dimensionless parameters: $Re = 75.71$, $Fr = 4.96$ and $\lambda = 1$.

In order to capture the influence of the yield stress on the waves, we consider two cases, namely $\tau_c = 3$ Pa and $\tau_c = 20$ Pa, corresponding to $Bi = 0.1$ and $Bi = 0.6$ respectively. Figure 4 shows the evolution of the simulated free-surface height $h$ and flow entrophy $\varphi$. Typical roll waves, with a discontinuous shock at the front, develop for both values of $Bi$. It is observed in Figure 4a that the amplitude and the wavelength of these roll waves decreases as the Bingham number is increased. The roll waves are also associated to marked variations of entrophy, whose values become smaller as the Bingham number increases (Figure 4b).

Figure 5 shows close-ups on the shape of a roll wave. It is seen that the maximum amplitude of the wave is not reached by the shock, but that the free-surface height continues to grow upstream of the shock (Figure 5a). This feature can be attributed to the presence of entrophy in the model, and contrasts with the predictions of two-equations models, for which the peak of the wave is reached exactly at the shock (Balmforth and Liu [12]). A similar behaviour was reported for Newtonian fluids by [45], and shown to be in good agreement with experimental data. It can also be noted that the entrophy is strongly influenced by the shock, with a marked peak at the front of the wave (Figure 5b).

## 6.4. Interpretation of the entrophy variations

In this section, we relate the entrophy variations observed within the roll waves (see Figure 5) to the shape of the reconstructed velocity profile. More precisely, we analyze the deviations of the entrophy from its equilibrium value $\varphi_{eq}$ given by Eq. (4.21):

$$\varphi_{eq} = \frac{U^2}{5h^2} \left(1 - \frac{\tau_c}{\rho g h \sin \theta}\right) \left(1 + \frac{5\tau_c}{4\rho g h \sin \theta}\right) \left(1 + \frac{\tau_c}{2\rho g h \sin \theta}\right)^{-2} \qquad (6.9)$$
<figure><img src="image_6.png" /><figcaption>Figure 4: Development of roll waves as a result of uniform flow instability for two values of Bi: evolution of (a) depth h and (b) entrophy φ as functions of the distance from the system entrance (Re = 75.71, Fr = 4.96, λ = 1, and θ = 18°).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 5: Shape of the roll waves: evolution depth h and entrophy φ (Re = 75.71, Fr = 4.96, Bi = 0.1, λ = 1 and θ = 18°).</figcaption></figure>

As a reference, let us start with the Newtonian case (Bi = 0). Alekseenko et al. [48] measured the velocity profiles in a wavy Newtonian film falling on a vertical wall. These authors distinguished four types regions in the wave. In region I, the "velocity profile is described by the self-similar parabolic law". In region II, "the velocity profile is less filled as compared to the parabolic one", while in region III, the velocity profile is "more filled". In region IV, no velocity profile could be determined due to scatter in the experimental points. From the back to the front of the wave, the authors first reported a region I, then a region II, a region I again, a region III located near the maximum depth, and finally a region IV at the front. Similar observations were later made by Denner et al. [49] from experimental data and DNS. As demonstrated by Richard et al. [28], these different zones can be directly related to the deviation of the entrophy from its equilibrium value. In the Newtonian case, φ<sub>eq</sub> = U<sup>2</sup>/5h<sup>2</sup>, and the longitudinal velocity u at order O(ε) is given by (see Eq. (5.19)):

$$u = U f(z/h; 0) + \frac{7}{270} \left( \frac{g \sin \theta}{\nu} \right)^3 \frac{h^4}{\varphi^2} \left[ \varphi - \frac{U^2}{5h^2} \right] \psi(z/h; 0). \quad (6.10)$$

The first term in this expression corresponds to the classical parabolic profile obtained, e.g., in a steady uniform flow, while the second term is zero for φ = φ<sub>eq</sub>. Hence, for φ ≈ φ<sub>eq</sub>, a parabolic velocity profile is recovered, corresponding to region I. For φ > φ<sub>eq</sub>, the velocity profile is "less filled" and corresponds to region II, while for φ < φ<sub>eq</sub>, the velocity profile is "more filled" and corresponds to region III. Comparisons between φ and φ<sub>eq</sub> in roll waves simulated with our model for Bi = 0 are shown in Figure 6, and the corresponding reconstructed velocity profiles are presented in Figure 7. It is observed that the succession of the zones along the wave, as delineated from the deviation between φ and φ<sub>eq</sub> (Figure 6a), and the corresponding differences in the shape of the velocity profiles, are in good agreement with the experimental observations of Alekseenko et al. [48]. Note nevertheless that region III ends before the peak of the wave, which was not the case in the experiments. This discrepancy is caused by the strong variation of entrophy at the shock and, as shown by Richard et al. [28], can be alleviated by adding a diffusion term to the model. Indeed, the
<figure><img src="image_8.png" /><figcaption>**Figure 6:** (a) Different regions within a roll wave defined from the deviation between the entrophy $\varphi$ and its equilibrium value $\varphi_{eq}$: Newtonian case ($Re = 75.71$, $Fr = 4.96$, $Bi = 0$, $\lambda = 1$ and $\theta = 18^\circ$). Values of $\varphi_{eq}$ are computed from the local depth $h$ and average velocity $U$ along the wave (see text). (b) Corresponding free-surface profile.</figcaption></figure>

shock and associated entrophy discontinuity are suppressed by the additional diffusion, which leads to an extension of region III up to the peak of the wave.

<figure><img src="image_9.png" /><figcaption>**Figure 7:** Typical longitudinal velocity profiles reconstructed in regions I, II and III in the Newtonian case (see Figure 6). The inset shows close-ups near the free surface.</figcaption></figure>

Turning to the case of a Bingham fluid, the equilibrium entrophy is obtained from Eq. (4.21), $\varphi_{eq} = \alpha_2(\xi)U^2/5h^2$, and the longitudinal velocity at order $O(\epsilon)$ is given by Eq. (5.19). The first term in (5.19) corresponds to the classical parabolic velocity profile overlaid by an unsheared plug zone (see also section 3.1). However, unlike in the Newtonian case, the deviation of the entrophy $\varphi$ from the equilibrium value $\varphi_{eq}$ does not fully define the type of velocity profile here. Indeed, (5.19) includes a second relaxation term corresponding to the deviation of the averaged velocity $U$ from the equilibrium value $U_{eq} = \alpha_1(\xi)(\rho g h^2 \sin \theta - \tau_c h)/3\rho v$. Nevertheless, we observed that the regions where $U < U_{eq}$ (resp. $U > U_{eq}$) approximately correspond to regions where $\varphi < \varphi_{eq}$ and (resp. $\varphi > \varphi_{eq}$). Again, comparisons between $\varphi$ and $\varphi_{eq}$ along a simulated roll wave and corresponding reconstructed velocity profiles are shown in Figures 8 and 9. Globally, similar trends as for Newtonian fluids are recovered. In zones for which $\varphi \approx \varphi_{eq}$, a parabolic velocity profile with an unsheared plug zone is observed (regions I). In zones for which $\varphi > \varphi_{eq}$, the velocity profile is "less filled" than the equilibrium profile (region II), while for $\varphi > \varphi_{eq}$ the velocity profile is more filled (region III). Note also that in region III, the longitudinal velocity profile displays a slight negative shearing in the pseudo-plug. The succession of the regions along the waves, namely I-II-I-III, is also found to be similar as for the Newtonian case (Figure 8a). Here again, the addition of a diffusion term to the model would likely decrease the peak of entrophy near the front of the wave, and increase the extent of region III up to the front of the wave.
<figure><img src="image_10.png" /><figcaption>**Figure 8:** (a) Different zones within a roll wave defined from the deviation between the entrophy $\varphi$ and its equilibrium value $\varphi_{eq}$: Bingham case ($Re = 75.71$, $Fr = 4.96$, $Bi = 0.1$, $\lambda = 1$, $h_p = 0.1$ and $\theta = 18^\circ$). Values of $\varphi_{eq}$ are computed from the local depth $h$ and average velocity $U$ along the wave (see text). (b) Corresponding free-surface profile.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9:** Typical longitudinal velocity profiles reconstructed in regions I, II and III in the Newtonian case (see Figure 8). The inset shows close-ups near the free surface.</figcaption></figure>

# 7. Conclusion

In this paper, a three-equation shallow-flow model for a Bingham fluid propagating down an inclined plane is consistently derived from the governing equations. The derivation of the model is based on a new asymptotic solution describing the flow composed of a sheared layer at the base and a pseudo-plug zone, in which the strain-rate is of order $O(\epsilon)$, close to the free surface. In contrast to previous approaches, the expansion is constructed with a regularized rheology based on a perfectly rigid behavior below the yield point. A yield-stress tensor is introduced, which is here completely determined by a few additional assumptions. Namely, the norm of the yield-stress tensor is assumed to be equal to the yield stress, its trace is assumed to be zero, and we assume that there are no normal stress differences in simple shear flows. Furthermore, the existence of a normal stress difference at order zero in the pseudo-plug layer, and all terms originating from this effect, are attributed to the yield-stress tensor rather than to the viscous part of the stress tensor. With this specific regularization, unlike in the classical approach, the yield-stress tensor and the strain-rate tensor are not assumed to be aligned, although they are still aligned in the particular case of a simple shear flow.

As a consequence, shearing in the pseudo-plug is related to the contribution of viscous stress terms, which allows us to eliminate the divergence of the strain rate at the fake yield surface, and to obtain smooth longitudinal velocity profiles at order $O(\epsilon)$. There is no need for a transition layer or for matched asymptotic expansions, which considerably simplifies the problem. In addition, this asymptotic solution accounts for inertial terms at order $O(\epsilon)$ in the pseudo-plug, which was not the case in previous studies.

The final model includes the depth-averaged mass conservation equation, the depth-averaged momentum balance equation, and a depth-averaged energy balance equation obtained from the work-energy theorem. The variables of the model are the fluid depth, the average velocity and the entrophy, which is related to the deviation of the velocity with respect to its average value and represents the internal shearing of the flow. The velocity field within the flow can be
reconstructed directly from the variables of the model, the benefit of which is the absence of the derivatives of the free-surface in the corresponding formulas.

The derived three-equation model can be written in conservative form and has the same mathematical structure as Euler equations for compressible fluids with relaxation terms. As a result, the model can be solved by classical and robust numerical schemes with relatively low computational cost. In contrast, the two-equation model of Fernández-Nieto et al. [26], which was the only other consistent shallow-water model for Bingham fluids to date, has a more complex mathematical structure and does not admit the energy balance. Moreover, this former model is derived on the base of non-smooth asymptotic velocity profiles, which precludes an accurate reconstruction of the velocity field at order $O(\epsilon)$.

Several applications of the derived model are presented. The linear stability analysis of equilibrium flows demonstrates the stabilizing effect of plasticity. Furthermore, the instability threshold is close to the result of Balmforth and Liu [12] obtained from the generalized Orr-Sommerfeld equations for viscoplastic fluids with no shear in the pseudo-plug. Importantly, this analysis supports the specific regularization of the constitutive law proposed in our study, as the instability thresholds derived with the classical approach [26] lead, on the contrary, to a nonphysical destabilizing effect of plasticity at large slopes. The model is then solved numerically to simulate the roll waves appearing above the instability threshold. It is shown that the amplitude and the wavelength of these rolls waves decrease as the yield stress grows. The variations of enstrophy along the waves is also analyzed, demonstrating that deviations of this quantity from its equilibrium value characterize the type of shearing within the flow. Namely, if the enstrophy is larger than its equilibrium value, the shearing is positive in the pseudo-plug; if the enstrophy is equal to its equilibrium value, the shearing is almost zero in the pseudo-plug; and if the enstrophy is smaller than its equilibrium, the shearing is negative in the pseudo-plug. The true physical relevance of this negative shearing should however be analyzed further, as its magnitude remains relatively small in the presented simulations and artifacts related to the shallow-flow approximation could play a role. To clarify this issue, further developments shall consider the inclusion of additional diffusive terms to the model, to smooth out the shock at the front of the waves. Consideration of non-hydrostatic effects might also be helpful to enrich the physics captured by the model.

Future works will also concentrate to extending the three-equation approach to the general case of a Herschel-Bulkley fluid and to three-dimensional flows, which will allow us to make direct comparisons with experimental data. In addition, the derivation of the model highlighted the important role played by the normal stress components in such viscoplastic free-surface flows. The assumption made regarding the small magnitude of normal stresses in the sheared layer may be questioned in the light of recent experimental studies [42, 43, 50] showing that normal stresses might actually be as large as shear stress in sheared viscoplastic flows. However, a properly validated 3D constitutive law accounting for these effects, which would also be necessary to fully specify the yield-stress tensor, is still lacking at the moment.

## Acknowledgements

This work has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 955605.

## Appendix A. Derivation of the energy equation

Starting from the dimensionless form of the work-energy theorem (4.28) and dropping terms of order $O(\epsilon^2)$ or smaller, one obtains

$$ \frac{\partial}{\partial t} \left( \frac{u^2}{2} \right) + \frac{\partial}{\partial x} \left[ u \left( \frac{u^2}{2} - \frac{x \tan \theta}{\epsilon Fr^2} + \frac{z}{Fr^2} \right) + \frac{pu}{Fr^2} - \frac{Bi}{Re} \tau_{xx} u \right] \\ + \frac{\partial}{\partial z} \left[ w \left( \frac{u^2}{2} - \frac{x \tan \theta}{\epsilon Fr^2} + \frac{z}{Fr^2} \right) + \frac{pw}{Fr^2} - \frac{Bi}{\epsilon Re} \tau_{xz} u + \frac{Bi}{Re} \tau_{xx} w \right] \\ = - \frac{2Bi}{Re} \tau_{xx}^Y \frac{\partial u}{\partial x} - \frac{Bi}{\epsilon Re} \left( \tau_{xz}^Y + \frac{1}{Bi} \frac{\partial u}{\partial z} \right) \frac{\partial u}{\partial z} \quad (A.1) $$

Introducing the leading-order representation for the pressure (3.11), (3.17) and the normal stress (3.10) and averaging equation (A.1) over the depth leads to:
$$ \frac{\partial}{\partial t} \left( \frac{h \langle u^2 \rangle}{2} + \frac{h^2}{2Fr^2} \right) + \frac{\partial}{\partial x} \left( \frac{h \langle u^3 \rangle}{2} + \frac{h^2 U}{Fr^2} - 2 \frac{Bi}{Re} \int_{h-h_p}^{h} \tau_{xx}^{Y(0)} u dz \right) = \frac{\lambda h U}{\varepsilon Re} - \frac{2Bi}{Re} \int_{h-h_p}^{h} \tau_{xx}^{Y(0)} \frac{\partial u}{\partial x} dz - \frac{Bi}{\varepsilon Re} \int_{0}^{h} \left( \tau_{xz}^{Y} + \frac{1}{Bi} \frac{\partial u}{\partial z} \right) \frac{\partial u}{\partial z} dz \quad (A.2) $$

Note that the integral terms involving the normal component of the yield stress tensor $\tau_{xx}^{Y(0)}$ in (A.2) are equal at leading order:

$$ \frac{\partial}{\partial x} \int_{h-h_p}^{h} \tau_{xx}^{Y(0)} u dz = \int_{h-h_p}^{h} \tau_{xx}^{Y(0)} \frac{\partial u}{\partial x} dz + O(\varepsilon) = \delta \frac{\lambda \pi}{4} h_p (h-h_p) \frac{\partial h}{\partial x}, \quad (A.3) $$

and therefore cancel each other. Using also definition (4.12) for the entrophy $\varphi$, equation (A.2) can be rewritten as

$$ \frac{\partial}{\partial t} \left( \frac{hU^2}{2} + \frac{h^3\varphi}{2} + \frac{h^2}{2Fr^2} \right) + \frac{\partial}{\partial x} \left( \frac{h\langle u^3 \rangle}{2} + \frac{h^2U}{Fr^2} \right) = \frac{\lambda h U}{\varepsilon Re} - \frac{Bi}{\varepsilon Re} \int_{0}^{h} \left( \tau_{xz}^{Y} + \frac{1}{Bi} \frac{\partial u}{\partial z} \right) \frac{\partial u}{\partial z} dz - \frac{\partial}{\partial x} \left( \frac{h\langle u^3 \rangle}{2} \right) \quad (A.4) $$

After calculations, the integral term in (A.4) at leading order can be expressed as

$$ \frac{Bi}{\varepsilon Re} \int_{0}^{h} \left( \tau_{xz}^{Y(0)} + \frac{1}{Bi} \frac{\partial u^{(0)}}{\partial z} \right) \frac{\partial u^{(0)}}{\partial z} dz = \frac{\lambda h U^{(0)}}{\varepsilon Re}, \quad (A.5) $$

while at order 1 we obtain

$$ \frac{Bi}{Re} \int_{0}^{h} \left( \tau_{xz}^{Y(0)} + \frac{1}{Bi} \frac{\partial u^{(0)}}{\partial z} \right) \frac{\partial u^{(1)}}{\partial z} dz = \frac{\lambda h U^{(1)}}{Re}. \quad (A.6) $$

Hence, we can write

$$ \frac{Bi}{\varepsilon Re} \int_{0}^{h} \left( \tau_{xz}^{Y} + \frac{1}{Bi} \frac{\partial u}{\partial z} \right) \frac{\partial u}{\partial z} dz = \frac{\lambda h U^{(0)}}{\varepsilon Re} + 2 \frac{\lambda h U^{(1)}}{Re} - \frac{Bi}{Re} \int_{0}^{h} \tau_{xz}^{Y(0)} \frac{\partial u^{(1)}}{\partial z} dz + O(\varepsilon), \quad (A.7) $$

which leads to the averaged energy equation expressed in (4.29).

## Appendix B. Newtonian case

For a Newtonian fluid ($\tau_c = 0$), the system of equations (4.40)-(4.42) reduces to

$$ \frac{\partial h}{\partial t} + \frac{\partial hU}{\partial x} = 0 \quad (B.1) $$

$$ \frac{\partial}{\partial t} (hU) + \frac{\partial}{\partial x} (hU^2 + \Pi) = \left[ gh \sin \theta - 3\nu \frac{U}{h} \right] + \frac{7 g \bar{h} \sin \theta}{6 \varphi} \left[ \varphi - \frac{U^2}{5h^2} \right] \quad (B.2) $$
$$ \frac{\partial}{\partial t}(he) + \frac{\partial}{\partial x}\left[hU\left(e + \frac{\Pi}{h}\right)\right] = U\left[g\sin\theta - 3v\frac{U}{h}\right] + \frac{7}{18}\frac{Ugh\sin\theta}{\varphi}\left[\varphi - \frac{U^2}{5h^2}\right]. \quad (B.3) $$

The difference with the Newtonian model derived by Richard et al. [28] is that the entrophy relaxes here toward $U^2/5h^2$, while it relaxes toward the term $g \sin \theta h^2/45v^2$ in this former model. The two expressions are equivalent except for modelling rest states ($U = 0$). In particular, the expression used by Richard et al. [28] is at the origin of nonphysical sources of momentum and energy at rest, which is not the case with the present model.

## Appendix C. Derivation of the linear instability criterion

We linearize the system of equations (4.3), (4.22) and (4.35) by considering small sinusoidal perturbations around the base flow (4.2) and (4.14): $h = 1 + h'$, $U = U^{(0)} + U'$ and $\varphi = \varphi^{(0)} + \varphi'$, with $[h', U', \varphi']^T = [A_1, A_2, A_3]^T \exp[ik(x - ct)]$. The three linearized equations can be written as:

$$ \begin{bmatrix} U^{(0)} - c & 1 & 0 \\ L_h^m - R_h^m & L_U^m - R_U^m & L_\varphi^m - R_\varphi^m \\ L_h^e - R_h^e & L_U^e - R_U^e & L_\varphi^e - R_\varphi^e \end{bmatrix} \begin{bmatrix} h' \\ U' \\ \varphi' \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 0 \end{bmatrix} \quad (C.1) $$

The coefficients $L_j^m$ ($j = h, U, \varphi$) come from the left-hand side of the momentum equation (4.22), and are given by:

$$ L_h^m = i\epsilon k \left( 3\varphi^{(0)} + \frac{1}{Fr^2} \right), \quad (C.2) $$

$$ L_u^m = i\epsilon k (U^{(0)} - c), \quad (C.3) $$

$$ L_\varphi^m = i\epsilon k. \quad (C.4) $$

The coefficients $L_j^e$ ($j = h, U, \varphi$) come from the left-hand side of the energy equation (4.35), and are given by:

$$ L_h^e = \frac{i\epsilon k}{2} \left[ (U^{(0)})^2 (U^{(0)} - c) + 3\varphi^{(0)} (3U^{(0)} - c) + \frac{2}{Fr^2} (2U^{(0)} - c) \right], \quad (C.5) $$

$$ L_u^m = \frac{i\epsilon k}{2} \left[ 3(U^{(0)})^2 + 3\varphi^{(0)} + \frac{2}{Fr^2} - 2U^{(0)}c \right], \quad (C.6) $$

$$ L_\varphi^e = \frac{i\epsilon k}{2} (3U^{(0)} - c). \quad (C.7) $$

The coefficients $R_j^m$ ($j = h, U, \varphi$) come from the right-hand side of the momentum equation (4.22), and are given by:

$$ R_h^m = \frac{2\lambda}{Re} \left(1 + \frac{1}{4}h_p + \frac{1}{4}h_p^2\right) \left(1 + \frac{1}{2}h_p\right)^{-1} \left[ \alpha_1(h_p) + \frac{7}{360} \frac{\lambda^2}{\varphi^{(0)}} \beta_1(h_p) \right] \\ + \frac{7}{135} \frac{\lambda^3}{\varphi^{(0)}Re} \left(1 + \frac{3}{8}h_p - \frac{39}{16}h_p^2 - \frac{5}{8}h_p^3\right) \left(1 - 2h_p + h_p^2\right)^2 \left(1 + \frac{1}{2}h_p\right)^{-1} \beta_2(h_p), \quad (C.8) $$

$$ R_U^m = -\frac{3}{Re} (1-h_p)^{-1} (1+\frac{1}{2}h_p)^{-1} \left[ \alpha_1(h_p) + \frac{7}{360} \frac{\lambda^2}{\varphi^{(0)}} \beta_1(h_p) \right] \\ - \frac{7}{45} \frac{\lambda^2}{\varphi^{(0)}Re} (1-h_p)^3 (1+\frac{5}{4}h_p) (1+\frac{1}{2}h_p)^{-1} \beta_2(h_p), \quad (C.9) $$

$$ R_\varphi^m = \frac{7}{6} \frac{\lambda}{Re} \frac{\beta_2(h_p)}{\varphi^{(0)}}, \quad (C.10) $$
Lastly, the coefficients $R_j^e$ ($j = h, U, \varphi$) come from the right-hand side of the energy equation (4.35), and are given by:

$$R_h^e = \frac{2\lambda U^{(0)}}{Re} \left(1 + \frac{1}{4}h_p + \frac{1}{4}h_p^2\right) \left(1 + \frac{1}{2}h_p\right)^{-1} \left[\alpha_1(h_p) + \frac{7}{1080} \frac{\lambda^2}{\varphi^{(0)}} \beta_1(h_p)r(h_p)\right] \\ + \frac{7}{405} \frac{\lambda^3}{Re} \frac{U^{(0)}}{\varphi^{(0)}} \left(1 + \frac{3}{8}h_p - \frac{39}{16}h_p^2 - \frac{5}{8}h_p^3\right) \left(1 - 2h_p + h_p^2\right)^2 \left(1 + \frac{1}{2}h_p\right)^{-1} \beta_2(h_p)r(h_p) \quad (C.11)$$

$$R_U^m = -\frac{3U^{(0)}}{Re} (1-h_p)^{-1} (1+\frac{1}{2}h_p)^{-1} \left[ \alpha_1(h_p) + \frac{7}{1080} \frac{\lambda^2}{\varphi^{(0)}} \beta_1(h_p)r(h_p) \right] \\ - \frac{\lambda^2}{Re} \frac{7}{135} \frac{U^{(0)}}{\varphi^{(0)}} (1-h_p)^3 (1+\frac{5}{4}h_p) (1+\frac{1}{2}h_p)^{-1} \beta_2(h_p)r(h_p) \quad (C.12)$$

$$R_{\varphi}^{m} = \frac{7}{18} \frac{\lambda}{Re} \frac{U^{(0)}}{\varphi^{(0)}} \beta_{2}(h_{p})r(h_{p}) \quad (C.13)$$

The dispersion relation (6.4) is obtained by equating the determinant of the system (C.1) to zero to have a non-trivial solution:

$$\left[ i R_{\varphi}^{e} + \varepsilon k \left( c - 3U^{(0)} \right) / 2 \right] \left[ i R_{h}^{m} + i \left( c - U^{(0)} \right) R_{U}^{m} + \varepsilon k \left( c^{2} - 1/F r^{2} \right) - \varepsilon k \left( 2c - U^{(0)} \right) U^{(0)} - 3\varepsilon k \varphi_{0} \right] \\ - \left[ i R_{\varphi}^{m} - \varepsilon k \right] \left[ i R_{h}^{e} + i \left( c - U^{(0)} \right) R_{U}^{e} - \varepsilon k \left( 2c - U^{(0)} \right) \left( U^{(0)} \right)^{2} + \varepsilon k \left( c^{2} - 1/F r^{2} - 3\varphi^{(0)} \right) U^{(0)} \right] = 0 \quad (C.14)$$

In the long-wave limit we can write: $c = c_0 + \varepsilon kc_1 + O(\varepsilon^2)$. Substituting this expansion into relation (C.14), we obtain after calculations:

$$c_0 = \lambda(1 - h_p), \quad (C.15)$$

$$c_1 = \frac{iRe}{3} \left[ \frac{2\lambda^2}{5} (1 - h_p)^2 \left( 1 + h_p + h_p^2 - \frac{h_p^3}{4} - \frac{h_p^4}{4} \right) - \frac{1}{Fr^2} \right], \quad (C.16)$$

which corresponds to (6.4).

## References

[1] N. J. Balmforth, I. A. Frigaard, G. Ovarlez, Yielding to stress: recent developments in viscoplastic fluid mechanics, Annual review of fluid mechanics 46 (2014) 121–146.

[2] C. Ancey, Plasticity and geophysical flows: a review, Journal of Non-Newtonian Fluid Mechanics 142 (2007) 4–35.

[3] I. Frigaard, Simple yield stress fluids, Current Opinion in Colloid & Interface Science 43 (2019) 80–93.

[4] I. Frigaard, C. Nouar, On the usage of viscosity regularisation methods for visco-plastic fluid flow computation, Journal of Non-Newtonian Fluid Mechanics 127 (2005) 1–26.

[5] R. Glowinski, A. Wachs, On the numerical simulation of viscoplastic fluid flow, in: Handbook of numerical analysis, volume 16, Elsevier, 2011, pp. 483–717.

[6] Y. Dimakopoulos, M. Pavlidis, J. Tsamopoulos, Steady bubble rise in herschel–bulkley fluids and comparison of predictions via the augmented lagrangian method with those via the papanastasiou model, Journal of Non-Newtonian Fluid Mechanics 200 (2013) 34–51.

[7] P. Saramito, A. Wachs, Progress in numerical simulation of yield stress fluid flows, Rheologica Acta 56 (2017) 211–230.

[8] Y. Liu, N. Balmforth, S. Hormozi, Viscoplastic surges down an incline, Journal of Non-Newtonian Fluid Mechanics 268 (2019) 1–11.

[9] G. Whitham, Linear and nonlinear waves, New York, Wiley-Interscience, 1974.

[10] S. Kalliadasis, C. Ruyer-Quil, B. Scheid, M. G. Velarde, Falling liquid films, volume 176, Springer Science & Business Media, 2011.

[11] C. Ruyer-Quil, P. Manneville, Improved modeling of flows down inclined planes, The European Physical Journal B-Condensed Matter and Complex Systems 15 (2000) 357–369.

[12] N. Balmforth, J. Liu, Roll waves in mud, Journal of Fluid Mechanics 519 (2004) 33–54.

[13] B. Hunt, Newtonian fluid mechanics treatment of debris flows and avalanches, Journal of Hydraulic Engineering 120 (1994) 1350–1363.

[14] C. Ruyer-Quil, P. Manneville, Modeling film flows down inclined planes, The European Physical Journal B-Condensed Matter and Complex Systems 6 (1998) 277–292.
[15] X. Huang, M. H. Garcia, A herschel-bulkley model for mud flow down a slope, Journal of fluid mechanics 374 (1998) 305-333.

[16] N. J. Balmforth, R. V. Craster, A. C. Rust, R. Sassi, Viscoplastic flow over an inclined surface, Journal of Non-Newtonian Fluid Mechanics 139 (2006) 103-127.

[17] C. Ancey, S. Cochard, The dam-break problem for herschel-bulkley viscoplastic fluids down steep flumes, Journal of Non-Newtonian Fluid Mechanics 158 (2009) 18-35.

[18] N. Bernabeu, P. Saramito, C. Smutek, Numerical modeling of non-newtonian viscoplastic flows: Part ii. viscoplastic fluids and general tridimensional topographies, International Journal of Numerical Analysis and Modeling 11 (2014) 213-228.

[19] D. Benney, Long waves on liquid films, Journal of mathematics and physics 45 (1966) 150-155.

[20] A. Pumir, P. Manneville, Y. Pomeau, On solitary waves running down an inclined plane, Journal of Fluid Mechanics 135 (1983) 27-50.

[21] T. Ooshida, Surface equation of falling film flows with moderate reynolds number and large but finite weber number, Physics of fluids 11 (1999) 3247-3269.

[22] R. Usha, B. Uma, Modeling of stationary waves on a thin viscous film down an inclined plane at high reynolds numbers and moderate weber numbers using energy integral method, Physics of Fluids 16 (2004) 2679-2696.

[23] P. Noble, J.-P. Vila, Thin power-law film flow down an inclined plane: consistent shallow-water models and stability under large-scale perturbations, Journal of Fluid Mechanics 735 (2013) 29-60.

[24] M. Boutounet, J. Monnier, J.-P. Vila, Multi-regime shallow free surface laminar flow models for quasi-newtonian fluids, European Journal of Mechanics-B/Fluids 55 (2016) 182-206.

[25] A. Chesnokov, Formation and evolution of roll waves in a shallow free surface flow of a power-law fluid down an inclined plane, Wave Motion 106 (2021) 102799.

[26] E. D. Fernández-Nieto, P. Noble, J.-P. Vila, Shallow water equations for non-newtonian fluids, Journal of Non-Newtonian Fluid Mechanics 165 (2010) 712-732.

[27] G. L. Richard, M. Gisclon, C. Ruyer-Quil, J.-P. Vila, Optimization of consistent two-equation models for thin film flows, European Journal of Mechanics-B/Fluids 76 (2019) 7-25.

[28] G. Richard, C. Ruyer-Quil, J.-P. Vila, A three-equation model for thin films down an inclined plane, Journal of Fluid Mechanics 804 (2016) 162-200.

[29] I. Walton, S. Bittleston, The axial flow of a bingham plastic in a narrow eccentric annulus, Journal of Fluid Mechanics 222 (1991) 39-60.

[30] P. Coussot, Steady, laminar, flow of concentrated mud suspensions in open channel, Journal of Hydraulic Research 32 (1994) 535-559.

[31] G. Chambon, A. Ghemmour, D. Laigle, Gravity-driven surges of a viscoplastic fluid: an experimental study, Journal of Non-Newtonian Fluid Mechanics 158 (2009) 54-62.

[32] N. Balmforth, R. Craster, A consistent thin-layer theory for bingham plastics, Journal of Non-Newtonian Fluid Mechanics 84 (1999) 65-81.

[33] G. Chambon, P. Freydier, M. Naaim, J.-P. Vila, Asymptotic expansion of the velocity field within the front of viscoplastic surges: comparison with experiments, Journal of Fluid Mechanics 884 (2020) A43.

[34] K. v. Hohenemser, W. Prager, Über die ansätze der mechanik isotroper kontinua, ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik 12 (1932) 216-226.

[35] J. G. Oldroyd, A rational formulation of the equations of plastic flow for a bingham solid, in: Mathematical Proceedings of the Cambridge Philosophical Society, volume 43, Cambridge University Press, 1947, pp. 100-105.

[36] P. Saramito, A new constitutive equation for elastoviscoplastic fluid flows, Journal of Non-Newtonian Fluid Mechanics 145 (2007) 1-14.

[37] C. J. Dimitriou, G. H. McKinley, A canonical framework for modeling elasto-viscoplasticity in complex fluids, Journal of Non-Newtonian Fluid Mechanics 265 (2019) 116-132.

[38] P. Coussot, Bingham's heritage, Rheologica Acta 56 (2017) 163-176.

[39] S. Varchanis, S. J. Haward, C. C. Hopkins, A. Syrakos, A. Q. Shen, Y. Dimakopoulos, J. Tsamopoulos, Transition between solid and liquid state of yield-stress fluids under purely extensional deformations, Proceedings of the National Academy of Sciences 117 (2020) 12611-12617.

[40] A. Geffrault, H. Bessaies-Bey, N. Roussel, P. Coussot, Extensional gravity-rheometry (EGR) for yield stress fluids, Journal of Rheology 65 (2021) 887-901.

[41] C. J. Dimitriou, R. H. Ewoldt, G. H. McKinley, Describing and prescribing the constitutive response of yield stress fluids using large amplitude oscillatory shear stress (LAOStress), Journal of Rheology 27 (2013) 27-70.

[42] R. L. Thompson, L. U. Sica, P. R. de Souza Mendes, The yield stress tensor, Journal of Non-Newtonian Fluid Mechanics 261 (2018) 211-219.

[43] H. De Cagny, M. Fazilati, M. Habibi, M. M. Denn, D. Bonn, The yield normal stress, Journal of Rheology 63 (2019) 285-290.

[44] V. Teshukov, Gas-dynamic analogy for vortex free-boundary flows, Journal of Applied Mechanics and Technical Physics 48 (2007) 303-309.

[45] G. L. Richard, S. L. Gavrilyuk, A new model of roll waves: comparison with brock's experiments, Journal of Fluid Mechanics 698 (2012) 374-405.

[46] D. M. Noma, S. Dagois-Bohy, S. Millet, V. Botton, D. Henry, H. B. Hadid, Primary instability of a visco-plastic film down an inclined plane: experimental study, Journal of Fluid Mechanics 922 (2021) R2.

[47] E. F. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics, Springer Berlin, Heidelberg, 1999.

[48] S. Alekseenko, V. Y. Nakoryakov, B. Pokusaev, Wave formation on a vertical falling liquid film, AIChE Journal 31 (1985) 1446-1460.

[49] F. Denner, A. Charogiannis, M. Pradas, C. N. Markides, B. G. Van Wachem, S. Kalliadasis, Solitary waves on falling liquid films in the inertia-dominated regime, Journal of Fluid Mechanics 837 (2018) 491-519.

[50] J.-M. Piau, Carbopol gels: Elastoviscoplastic and slippery glasses made of individual swollen sponges: Meso- and macroscopic properties, constitutive equations and scaling laws, Journal of Non-Newtonian Fluid Mechanics 144 (2007) 1-29.