domain: spi math

text:
Danila Denisenko
email: danila.denisenko@inrae.fr

Gaël Loïc

Guillaume Chambon

Richard Denisenko

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

Keywords:

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 removing the assumption of alignment between the yield-stress tensor and the strain-rate tensor, but keeping the traceless property and the equality between the norm of the yield-stress tensor and the yield stress. 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 enstrophy, 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 enstrophy 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.

Introduction

Viscoplastic materials behave like solid bodies when 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) [START_REF] Balmforth | Yielding to stress: recent developments in viscoplastic fluid mechanics[END_REF]. 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 [START_REF] Ancey | Plasticity and geophysical flows: a review[END_REF]. 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 behavior in the solid-like regime [START_REF] Frigaard | Simple yield stress fluids[END_REF].

The mathematical modeling of idealized viscoplastic fluids generally relies on Hershel-Bulkley or Bingham constitutive laws. Combined with 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 behavior in the unyielded zones by a highly-viscous flow [START_REF] Frigaard | On the use of viscosity regularisation methods for visco-plastic fluid flow computation[END_REF]. 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 [START_REF] Glowinski | On the numerical simulation of viscoplastic fluid flow[END_REF][START_REF] Dimakopoulos | Steady bubble rise in Herschel-Bulkley fluids and comparison of predictions via the augmented Lagrangian method with those via the Papanastasiou model[END_REF][START_REF] Saramito | Progress in numerical simulation of yield stress fluid flows[END_REF][START_REF] Liu | Viscoplastic surges down an incline[END_REF]. 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 [START_REF] Whitham | Linear and Nonlinear Waves[END_REF][START_REF] Kalliadasis | Falling Liquid Films[END_REF]. 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 fields of interests with respect to flow aspect ratio 𝜀 = ℎ 0 ∕𝑙 0 , where ℎ 0 and 𝑙 0 denote 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 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 𝑛 if the leading terms in the model equations are of 𝑂 [START_REF] Balmforth | Yielding to stress: recent developments in viscoplastic fluid mechanics[END_REF], and if all terms vanish except for a remainder of 𝑂(𝜀 𝑛+1 ) after having inserted the asymptotic expansions obtained above into the model equations. Inconsistent reduced models lead to inaccurate predictions of, e.g., instability thresholds [START_REF] Balmforth | Roll waves in mud[END_REF][START_REF] Ruyer-Quil | Improved modeling of flows down inclined planes[END_REF].

Depth-averaged models can be formulated as systems of one, two or three equations. One-equation models usually have the simplest structure and are obtained from mass conservation by enslaving the fluid velocity to the fluid height [START_REF] Hunt | Newtonian fluid mechanics treatment of debris flows and avalanches[END_REF][START_REF] Ruyer-Quil | Modeling film flows down inclined planes[END_REF][START_REF] Huang | A Herschel-Bulkley model for mud flow down a slope[END_REF][START_REF] Balmforth | Viscoplastic flow over an inclined surface[END_REF][START_REF] Ancey | The dam-break problem for Herschel-Bulkley viscoplastic fluids down steep flumes[END_REF][START_REF] Bernabeu | Numerical modeling of non-Newtonian viscoplastic flows: Part II. Viscoplastic fluids and general tridimensional topographies[END_REF]. However, consistent one-equation models generally produce diverging or inaccurate solutions when the instability threshold for uniform flows is exceeded [START_REF] Benney | Long waves on liquid films[END_REF][START_REF] Pumir | On solitary waves running down an inclined plane[END_REF][START_REF] Ooshida | Surface equation of falling film flows with moderate Reynolds number and large but finite Weber number[END_REF]. 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 [START_REF] Ruyer-Quil | Improved modeling of flows down inclined planes[END_REF][START_REF] Ruyer-Quil | Modeling film flows down inclined planes[END_REF][START_REF] Usha | 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[END_REF][START_REF] Noble | Thin power-law film flow down an inclined plane: Consistent shallow-water models and stability under largescale perturbations[END_REF][START_REF] Boutounet | Multi-regime shallow free surface laminar flow models for quasi-Newtonian fluids[END_REF][START_REF] Chesnokov | Formation and evolution of roll waves in a shallow free surface flow of a power-law fluid down an inclined plane[END_REF]. Although relatively rare, a few studies also considered the case of viscoplastic fluids [START_REF] Balmforth | Roll waves in mud[END_REF][START_REF] Fernandez-Nieto | Shallow Water equations for non-Newtonian fluids[END_REF]. 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 [START_REF] Richard | Vila Optimization of consistent two-equation models for thin film flows[END_REF] 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 [START_REF] Richard | A three-equation model for thin films down an inclined plane[END_REF] 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 ensure 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 powerlaw fluids, for viscoplastic fluids difficulties arise from the possible coexistence of yielded and unyielded regions within the flows. 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 [START_REF] Coussot | Steady, laminar, flow of concentrated mud suspensions in open channel[END_REF][START_REF] Chambon | Gravity-driven surges of a viscoplastic fluid: an experimental study[END_REF]. At the next order of approximation, Balmforth and Craster [START_REF] Balmforth | A consistent thin-layer theory for Bingham plastics[END_REF] 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 𝑂(𝜀). These authors derived a first-order correction for the longitudinal velocity profile in inertia-less limit. Later, Chambon et al [START_REF] Chambon | Asymptotic expansion of the velocity field within the front of viscoplastic surges: comparison with experiments[END_REF] constructed the full expressions for the longitudinal velocity up to the first order. However, the obtained asymptotic solution shows two main drawbacks: [START_REF] Balmforth | Yielding to stress: recent developments in viscoplastic fluid mechanics[END_REF] 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; [START_REF] Ancey | Plasticity and geophysical flows: a review[END_REF] 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 [START_REF] Balmforth | A consistent thin-layer theory for Bingham plastics[END_REF] and Fernandez-Nieto et al [START_REF] Fernandez-Nieto | Shallow Water equations for non-Newtonian fluids[END_REF] proposed to avoid the divergence of the strain rate by introducing a transition layer between the pseudo-plug and the sheared layer. Fernandez-Nieto et al [START_REF] Fernandez-Nieto | Shallow Water equations for non-Newtonian fluids[END_REF] also derived a two-equation model, which remains the only consistent shallow-flow model for idealized viscoplastic fluids to date. However, the extra terms arising from the transition layer are complicated and of 𝑂(𝜀 4∕3 ) order, such that the expansions providing the smooth transition at the fake yield surface were not considered in deriving this model.

The goal of this paper is to generalize the three-equation approach of Richard et al [START_REF] Richard | A three-equation model for thin films down an inclined plane[END_REF] for Bingham fluid flows propagating down an inclined plane in order to derive a consistent depth-averaged model with well-posed mathematical structure. We construct a new asymptotic solution up to the first order in 𝜀, based on a reformulated version of the tensorial constitutive law. This allows us to eliminate the issue with diverging the 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 enstrophy 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 on 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 [START_REF] Balmforth | Roll waves in mud[END_REF] from generalized Orr-Sommerfeld equations. Another important advantage of the three-equation approach is also highlighted, namely the full velocity field can be consistently reconstructed directly from the variables (flow height, averaged velocity and enstrophy) of the model, without any derivatives. In particular, we analyze the physical relevance of the enstrophy in terms of shapes of the velocity profiles within a roll wave.

In §2, we formulate the equations for the fluid flow. In §3, we construct the new shallow-flow asymptotic expansion up to 𝑂(𝜀) 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 enstrophy variations predicted by the model.

Formulation of the problem

We consider a two-dimensional flow of a viscoplastic fluid propagating down an inclined plane under gravity 𝒈 (Figure 1). The angle of the slope with respect to the horizontal is 𝜃. The directions 𝑂𝑥 and 𝑂𝑧 are parallel and normal to the plane, respectively. The corresponding components of the velocity field 𝒗 are denoted by 𝑢 and 𝑤, and the components of the strain-rate tensor γ are defined as: ̇𝛾𝑥𝑥 = 2𝜕𝑢∕𝜕𝑥, ̇𝛾𝑥𝑧 = 𝜕𝑢∕𝜕𝑧 + 𝜕𝑤∕𝜕𝑥, ̇𝛾𝑧𝑧 = 2𝜕𝑤∕𝜕𝑧. The fluid depth is denoted by ℎ(𝑥, 𝑡). Lastly, the fluid is assumed to be incompressible (tr γ = 0) with a density 𝜌.

Constitutive law

The fluid is assumed to obey the Bingham constitutive law. The relation between the extra-stress tensor 𝝉 = 𝝈 + 𝑝𝑰 (with 𝝈 the total stress and 𝑝 the pressure) and the strain-rate tensor γ can be expressed as follows:

𝜏 𝑖𝑗 = 𝜏 𝑌 𝑖𝑗 + 𝐾 ̇𝛾𝑖𝑗 |𝝉| > 𝜏 𝑐 , (2.1
)

̇𝛾𝑖𝑗 = 0 |𝝉| ⩽ 𝜏 𝑐 , ( 2.2)
where the rheological parameters 𝜏 𝑐 and 𝐾 correspond to the yield stress and the Bingham viscosity of the material, respectively. The tensor norm is defined as |𝑻 | = (0.5 𝑻 ∶ 𝑻 ) 0.5 for any second-order tensor 𝑻 , where the colon denotes the double dot product. Hence, the yielding of the fluid is governed by a von Mises criterion (|𝝉| = 𝜏 𝑐 ). In (2.1), the tensor 𝝉 𝑌 corresponds to the yield-stress (or plastic) contribution to the stress. Accordingly, the second term 𝜏 𝜈 𝑖𝑗 = 𝐾 ̇𝛾𝑖𝑗 corresponds to the viscous contribution.

In most studies on viscoplastic fluids [START_REF] Balmforth | A consistent thin-layer theory for Bingham plastics[END_REF][START_REF] Chambon | Asymptotic expansion of the velocity field within the front of viscoplastic surges: comparison with experiments[END_REF][START_REF] Piau | Flow of a yield stress fluid in a long domain. Application to flow on an inclined plane[END_REF], the yield-stress tensor 𝝉 𝑌 is expressed as follows:

𝝉 𝑌 = 𝜏 𝑐 γ | γ| . (2.3)
This expression, which assumes that 𝝉 𝑌 is aligned with the strain-rate tensor γ, was first proposed by Hohenemser and Prager [START_REF] Hohenemser | Über die Ansätze der Mechanik isotroper Kontinua[END_REF]. It corresponds to a tensorial extension of the scalar Bingham law used in simple shear, as it simply reduces to 𝜏 𝑌 𝑥𝑧 = 𝜏 𝑐 in this configuration. However, as mentioned in a number of recent studies [START_REF] Piau | Carbopol gels: elastoviscoplastic and slippery glasses made of indi-vidual swollen sponges. Meso-and macroscopic properties, constitutive equa-tions and scaling laws[END_REF][START_REF] Thompson | The yield stress tensor[END_REF][START_REF] De Cagny | The yield normal stress[END_REF], this formulation fails to capture the normal stress components that develop in viscoplastic fluids at yielding. Furthermore, as will be explained later (see section 3), we argue that this expression also leads to a singularity in the expression for the strain rate at the interface between yielded and unyielded regions. For the moment, we thus keep the yield-stress tensor 𝝉 𝑌 unspecified, and only consider the two following conditions for |𝝉| > 𝜏 𝑐 :

|𝝉 𝑌 | = 𝜏 𝑐 , tr 𝝉 𝑌 = 0. (2.4)
The first condition is required to have continuity of the stress at the interface between yielded and unyielded regions (von Mises criterion). The second condition is assumed in relation to fluid incompressibility (i.e. any isotropic part of the total stress tensor is assumed to contribute to the pressure). Note that these two conditions are obviously met in the case of the classical formulation (2.3).

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 continuity equation writes

𝜕𝑢 𝜕𝑥 + 𝜕𝑤 𝜕𝑧 = 0. (2.5)
The Cauchy momentum equations in the 𝑂𝑥 and 𝑂𝑧 directions are

𝜌 ( 𝜕𝑢 𝜕𝑡 + 𝑢 𝜕𝑢 𝜕𝑥 + 𝑤 𝜕𝑢 𝜕𝑧 ) = - 𝜕𝑝 𝜕𝑥 + 𝜌𝑔 sin 𝜃 + 𝜕𝜏 𝑥𝑧 𝜕𝑧 + 𝜕𝜏 𝑥𝑥 𝜕𝑥 , ( 2.6
)

𝜌 ( 𝜕𝑤 𝜕𝑡 + 𝑢 𝜕𝑤 𝜕𝑥 + 𝑤 𝜕𝑤 𝜕𝑧 ) = - 𝜕𝑝 𝜕𝑧 -𝜌𝑔 cos 𝜃 + 𝜕𝜏 𝑥𝑧 𝜕𝑥 + 𝜕𝜏 𝑧𝑧 𝜕𝑧 . (2.7)
According to the constitutive equation (2.1) and the conditions on the yield-stress tensor (2.4), the stress components, for |𝝉| > 𝜏 𝑐 , can be written as

𝜏 𝑥𝑥 = -𝜏 𝑧𝑧 = 𝜏 𝑌 𝑥𝑥 + 𝜏 𝜈 𝑥𝑥 , 𝜏 𝜈 𝑥𝑥 = 2𝐾 𝜕𝑢 𝜕𝑥 , ( 2.8
)

𝜏 𝑥𝑧 = 𝜏 𝑌 𝑥𝑧 + 𝜏 𝜈 𝑥𝑧 , 𝜏 𝜈 𝑥𝑧 = 𝐾 ( 𝜕𝑢 𝜕𝑧 + 𝜕𝑤 𝜕𝑥 ) , ( 2.9)
with the relation for the norm of the yield stress tensor:

( 𝜏 𝑌 𝑥𝑥 ) 2 + ( 𝜏 𝑌 𝑥𝑧 ) 2 = 𝜏 2 𝑐 .
(2.10)

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

𝑢 | 𝑧=0 = 𝑤 | 𝑧=0 = 0. (2.11)
At the free surface 𝑧 = ℎ(𝑥), the following kinematic boundary condition holds:

𝜕ℎ 𝜕𝑡 + 𝑢 | 𝑧=ℎ(𝑥) 𝜕ℎ 𝜕𝑥 = 𝑤 | 𝑧=ℎ(𝑥) .
(2.12)

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 write

[ 1 - ( 𝜕ℎ 𝜕𝑥 ) 2 ] 𝜏 𝑥𝑧| 𝑧=ℎ(𝑥) = 2 𝜕ℎ 𝜕𝑥 𝜏 𝑥𝑥| 𝑧=ℎ(𝑥) , (2.13) [ 1 - ( 𝜕ℎ 𝜕𝑥 ) 2 ] 𝑝 | 𝑧=ℎ(𝑥) = - [ 1 + ( 𝜕ℎ 𝜕𝑥 ) 2 ]
𝜏 𝑥𝑥| 𝑧=ℎ(𝑥) .

(2.14)

Shallow-flow scaling

Let us define ℎ 0 the characteristic depth of the flow in the 𝑂𝑧 direction, and 𝑢 0 the characteristic velocity in the 𝑂𝑥 direction. The characteristic length in the 𝑂𝑥 direction is denoted by 𝑙 0 . The shallow-flow hypothesis corresponds to assuming that the aspect ratio 𝜀 = ℎ 0 ∕𝑙 0 is small. The main dimensionless groups of this problem are the Reynolds number 𝑅𝑒, the Froude number 𝐹 𝑟 and the Bingham number 𝐵𝑖, which are defined as

𝑅𝑒 = 𝜌𝑢 0 ℎ 0 𝐾 , 𝐹 𝑟 = 𝑢 0 √ 𝑔ℎ 0 cos 𝜃 , 𝐵𝑖 = 𝜏 𝑐 ℎ 0 𝐾𝑢 0 . (2.15)
These parameters, as well as the slope angle 𝜃, are assumed to be of order 𝑂 [START_REF] Balmforth | Yielding to stress: recent developments in viscoplastic fluid mechanics[END_REF] with respect to aspect ratio 𝜀. In order to reformulate the problem (2. with the condition on the norm of the dimensionless yield-stress tensor:

( 𝜏 𝑌 𝑥𝑥 ) 2 + ( 𝜏 𝑌 𝑥𝑧 ) 2 = 1. (2.23)
The no-penetration and the no-slip conditions (2.11) keep the same form

𝑢 | 𝑧=0 = 𝑤 | 𝑧=0 = 0, (2.24)
while the dimensionless kinematic and dynamic boundary conditions at the free surface (2.12)-(2.14) are rewritten as

𝜕ℎ 𝜕𝑡 + 𝑢 | 𝑧=ℎ(𝑥) 𝜕ℎ 𝜕𝑥 = 𝑤 | 𝑧=ℎ(𝑥) , (2.25)
[ 1 -𝜀 2 ( 𝜕ℎ 𝜕𝑥 ) 2 ] 𝜏 𝑥𝑧| 𝑧=ℎ(𝑥) = 2𝜀 𝜕ℎ 𝜕𝑥 𝜏 𝑥𝑥| 𝑧=ℎ(𝑥) , (2.26)
[ 1 -𝜀 2 ( 𝜕ℎ 𝜕𝑥 ) 2 ] 𝑝 | 𝑧=ℎ(𝑥) = - 𝐵𝑖𝐹 𝑟 2 𝑅𝑒 [ 1 + 𝜀 2 ( 𝜕ℎ 𝜕𝑥 ) 2 ] 𝜏 𝑥𝑥| 𝑧=ℎ(𝑥) . (2.27)
Lastly, the norms of the dimensionless stress and strain-rate tensors express as

|𝝉| = √ 𝜏 2 𝑥𝑥 + 𝜏 2 𝑥𝑧 | γ| = √ ( 𝜕𝑢 𝜕𝑧 + 𝜀 2 𝜕𝑤 𝜕𝑥 ) 2 + 4𝜀 2 ( 𝜕𝑢 𝜕𝑥 ) 2 .
(2.28)

Asymptotic expansions

Let us assume the existence of regular expansions of the form 𝑓 = 𝑓 (0) + 𝜀𝑓 (1) + … (3.1) for all variables of the problem (2.17)-(2.27), namely longitudinal and normal velocities 𝑢 and 𝑤, pressure 𝑝 and stress components 𝜏 𝑖𝑗 . As mentioned in introduction, the structure of thin viscoplastic flows generally consist of a pseudoplug, in which the strain rate vanishes at leading order, overlying a sheared layer. These two layers are separated by a fake yield surface. We shall construct the expansions in these two layer separately, using tilde notations for the pseudo-plug. In this layer, the leading-order longitudinal velocity ũ(0) is assumed to be independent of the normal coordinate 𝑧 [START_REF] Balmforth | A consistent thin-layer theory for Bingham plastics[END_REF]:

𝑢 = ũ(0) (𝑥, 𝑡) + 𝜀 ũ(1) (𝑥, 𝑧, 𝑡) + … (3.2)
Accordingly, the strain rate in the pseudo-plug writes

| γ| = 𝜀 √ ( 𝜕 ũ(1) 𝜕𝑧 ) 2 + 4 ( 𝜕 ũ(0) 𝜕𝑥 ) 2 + 𝑂(𝜀 2 ). (3.3)

Leading-order expansion

At 𝑂 [START_REF] Balmforth | Yielding to stress: recent developments in viscoplastic fluid mechanics[END_REF] with respect to 𝜀, integration of the momentum equation (2.18) with the boundary condition (2.26) gives the following shear-stress profile for both the sheared and the pseudo-plug layers:

𝜏 (0) 𝑥𝑧 = 𝜆 𝐵𝑖 (ℎ -𝑧), (3.4) τ(0) 𝑥𝑧 = 𝜆 𝐵𝑖 (ℎ -𝑧). (3.5)
In the sheared layer, the flow behaves as a viscous fluid, for which the smallness of the normal stress results from the shallow-flow assumption. We keep here this result and assume that the leading-order normal stress is zero in the sheared layer:

𝜏 (0) 𝑥𝑥 = 0. (3.6)
This leads to |𝝉| (0) = (𝜆∕𝐵𝑖)(ℎ -𝑧) in the sheared layer. The thickness of the pseudo-plug ℎ 𝑝 can then be obtained from the yielding criterion |𝝉| (0) = 1 at the fake-yield surface 𝑧 = ℎℎ 𝑝 , which gives

ℎ 𝑝 = 𝐵𝑖 𝜆 .
Further we consider the expansions for the pseudo-plug and the sheared layer separately.

In the pseudo-plug (𝑧 ⩾ ℎℎ 𝑝 )

Owing to the expansions (3.2), the constitutive law (2.21)-(2.22) reduces to

τ(0) 𝑥𝑥 = τ𝑌 (0) 𝑥𝑥 , (3.7) τ(0) 𝑥𝑧 = τ𝑌 (0) 𝑥𝑧 . (3.8)
From the condition on the norm of the yield-stress tensor (2.23), one then obtain

τ(0) 𝑥𝑥 = 𝛿 √ 1 - ( ℎ -𝑧 ℎ 𝑝 ) 2 (3.9)
with 𝛿 = sgn(𝜏 𝑥𝑥 ). 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, | τ| (0) = 1. Note that these normal stresses vanish at the fake yield surface 𝑧 = ℎℎ 𝑝 . Further, the pressure profile is obtained from integration of the momentum equation (2.19) with the dynamic condition (2.27):

p(0) = ℎ -𝑧 -𝛿 𝐵𝑖𝐹 𝑟 2 𝑅𝑒 √ 1 - ( ℎ -𝑧 ℎ 𝑝 ) 2 . (3.10)
The expression for the longitudinal velocity ũ(0) (𝑥, 𝑡) will be given later from matching with the sheared zone. The normal velocity w(0) can then be derived by integration of the continuity equation (2.17):

w(0) = -𝑧 𝜕 ũ(0) 𝜕𝑥 + 𝑤 (0) + , ( 3.11)
where, again, the term 𝑤 (0) + will be obtained from matching.

In the sheared zone (𝑧 < ℎℎ 𝑝 )

In this zone, the constitutive law expresses as

𝜏 (0) 𝑥𝑥 = 𝜏 𝑌 (0) 𝑥𝑥 , (3.12) 𝜏 (0) 𝑥𝑧 = 𝜏 𝑌 (0) 𝑥𝑧 + 1 𝐵𝑖 𝜕𝑢 (0) 𝜕𝑧 . (3.13)
The condition on the smallness of the normal stress in the sheared layer (3.6) leads to 𝜏 𝑌 (0) 𝑥𝑥 = 0. Then the von Mises criterion (2.23) implies that 𝜏 𝑌 (0) 𝑥𝑧 = 1 and integration of equation (3.13) with the no-slip condition (2.24) results in a parabolic longitudinal velocity profile:

𝑢 (0) = 𝜆𝑧 ( ℎ -ℎ 𝑝 - 𝑧 2 ) . (3.14)
Note that 𝜕𝑢 (0) ∕𝜕𝑧 = 0 at 𝑧 = ℎℎ 𝑝 , consistently with the expansion (3.2) assumed in the pseudo-plug. The normal velocity 𝑤 (0) is then found from integration of equation (2.17) with the no-slip boundary condition (2.24):

𝑤 (0) = -𝜆𝑧 2 2 𝜕ℎ 𝜕𝑥 . (3.15)
Finally, integration of the momentum equation (2.19) and stress continuity at the fake yield surface 𝑧 = ℎℎ 𝑝 provide the pressure profile:

𝑝 (0) = ℎ -𝑧. ( 3.16)
Note that, unlike in the pseudo-plug, a classical hydrostatic pressure distribution is recovered in the sheared layer.

Matching

The continuity condition ũ(0) = 𝑢 (0) at the fake yield surface 𝑧 = ℎℎ 𝑝 provides the expression of the longitudinal velocity in the pseudo-plug:

ũ(0) = 𝜆 2 (ℎ -ℎ 𝑝 ) 2 . (3.17)
It is worth noting that the leading-order longitudinal velocity profile given by Eqs. (3.14) and (3.17) is identical to the profile that would be obtained in a steady uniform flow of height ℎ, 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 𝑂(𝜀).

Finally, matching the solutions 𝑤 (0) and w(0) at 𝑧 = ℎℎ 𝑝 gives the following expression for the normal velocity in the pseudo-plug:

w(0) = -𝜆(ℎ -ℎ 𝑝 ) ( 𝑧 - ℎ -ℎ 𝑝 2 ) 𝜕ℎ 𝜕𝑥 . (3.18)

O(𝜀) expansion

Here we construct the expansions at order 𝑂(𝜀) for the shear stress 𝜏 (1) 𝑥𝑧 and the longitudinal velocity 𝑢 (1) . The expansions of the other variables will not be needed for the derivation of a depth-averaged model consistent at first order. The leading-order solution derived above does not depend on whether one considers the classical expression (2.3) for the yield-stress tensor 𝝉 𝑌 , or one only assumes the conditions (2.4). This does not remain true, however, at 𝑂(𝜀). With the classical approach, equations (3.3) and (3.8) lead to the following relation in the pseudo-plug:

τ(0) 𝑥𝑧 = τ𝑌 (0) 𝑥𝑧 = 𝜕 ũ(1) 𝜕𝑧 √ ( 𝜕 ũ(1) 𝜕𝑧 ) 2 + 4 ( 𝜕 ũ(0) 𝜕𝑥 ) 2 (3.19)
and thus:

𝜕 ũ(1) 𝜕𝑧 = 2(ℎ -𝑧)∕ℎ𝑝 √ 1 -((ℎ -𝑧)∕ℎ 𝑝 ) 2 | | | | 𝜕 ũ(0) 𝜕𝑥 | | | | . (3.20)
Expression (3.20) was used in former studies to derive the first-order velocity correction ũ(1) in the pseudo-plug [START_REF] Fernandez-Nieto | Shallow Water equations for non-Newtonian fluids[END_REF][START_REF] Balmforth | A consistent thin-layer theory for Bingham plastics[END_REF][START_REF] Chambon | Asymptotic expansion of the velocity field within the front of viscoplastic surges: comparison with experiments[END_REF]. However, it is easily seen that 𝜕 ũ( 1) ∕𝜕𝑧 diverges at the fake-yield surface 𝑧 = ℎℎ 𝑝 , which contradicts the assumption that the strain-rate should remain small in the pseudo-plug. As already suggested by Balmforth and Craster [START_REF] Balmforth | A consistent thin-layer theory for Bingham plastics[END_REF] and explored by Fernandez-Nieto et al [START_REF] Fernandez-Nieto | Shallow Water equations for non-Newtonian fluids[END_REF], this inconsistency could be alleviated by introducing a transition layer between the sheared zone and the pseudo-plug. However, this significantly complicates the process of constructing the solution. Another drawback of expressions (3.19) (3.20) is that the first-order correction 𝜕 ũ(1) ∕𝜕𝑧 is controlled solely by terms related to the yield-stress tensor. Instead, it could be expected that the slight shearing of the pseudo-plug be rather associated to viscous stresses.

As already mentioned, these drawbacks motivated us to relax the assumption that 𝝉 𝑌 is aligned with the strain-rate tensor γ, and to only consider the conditions (2.4). As will be shown, these conditions together with the assumption of small normal stresses in the sheared layer are actually sufficient to build the solution at order 𝑂(𝜀), even though the yield-stress tensor 𝝉 𝑌 is not fully specified.

In the pseudo-plug (𝑧 > ℎℎ 𝑝 )

At order 𝑂(𝜀), the momentum equation along 𝑂𝑥 (2.18) leads to

𝐵𝑖 𝑅𝑒 𝜕 τ(1) 𝑥𝑧 𝜕𝑧 = 𝜕 ũ(0) 𝜕𝑡 + ũ(0) 𝜕 ũ(0) 𝜕𝑥 + 1 𝐹 𝑟 2 𝜕 p(0) 𝜕𝑥 - 𝐵𝑖 𝑅𝑒 𝜕 τ(0) 𝑥𝑥 𝜕𝑥 = 𝜆(ℎ -ℎ 𝑝 ) 𝜕ℎ 𝜕𝑡 + 𝜆 2 2 (ℎ -ℎ 𝑝 ) 3 𝜕ℎ 𝜕𝑥 + 1 𝐹 𝑟 2 𝜕ℎ 𝜕𝑥 -𝛿 2𝐵𝑖 𝑅𝑒 𝜕 𝜕𝑥 ⎛ ⎜ ⎜ ⎝ √ 1 - ( ℎ -𝑧 ℎ 𝑝 ) 2⎞ ⎟ ⎟ ⎠ . (3.21)
Note that the last term in (3.21) can be rewritten as

-2𝛿 𝜕 𝜕𝑥 ⎛ ⎜ ⎜ ⎝ √ 1 - ( ℎ -𝑧 ℎ 𝑝 ) 2⎞ ⎟ ⎟ ⎠ = 2𝛿 𝜕 𝜕𝑧 ⎛ ⎜ ⎜ ⎝ √ 1 - ( ℎ -𝑧 ℎ 𝑝 ) 2⎞ ⎟ ⎟ ⎠ 𝜕ℎ 𝜕𝑥 , ( 3
τ(1) 𝑥𝑧 = 𝑅𝑒 𝐵𝑖 (𝑧 -ℎ) ( 𝜆(ℎ -ℎ 𝑝 ) 𝜕ℎ 𝜕𝑡 + 𝜆 2 2 (ℎ -ℎ 𝑝 ) 3 𝜕ℎ 𝜕𝑥 ) + 𝑅𝑒 𝐵𝑖 𝐹 𝑟 2 (𝑧 -ℎ) 𝜕ℎ 𝜕𝑥 + 2𝛿 √ 1 - ( ℎ -𝑧 ℎ 𝑝 ) 2 𝜕ℎ 𝜕𝑥 . (3.24)
Recall that stresses express as the sum of a yield-stress and of a viscous contributions, i.e. τ(1) 𝑥𝑧 = τ𝑌 (1) 𝑥𝑧 + τ𝜈 (1) 𝑥𝑧 ∕𝐵𝑖. By identification, the first-order correction to the yield-stress term writes

τ𝑌 (1) 𝑥𝑧 = 2𝛿 √ 1 - ( ℎ -𝑧 ℎ 𝑝 ) 2 𝜕ℎ 𝜕𝑥 , ( 3.25)
while the first-order correction to the viscous contribution is

τ𝜈 (1) 𝑥𝑧 = 𝑅𝑒(𝑧 -ℎ) ( 𝜆(ℎ -ℎ 𝑝 ) 𝜕ℎ 𝜕𝑡 + 𝜆 2 2 (ℎ -ℎ 𝑝 ) 3 𝜕ℎ 𝜕𝑥 ) + 𝑅𝑒 𝐹 𝑟 2 (𝑧 -ℎ) 𝜕ℎ 𝜕𝑥 . (3.26)
In addition, from equation (2.22) it is found 𝜕 ũ(1) ∕𝜕𝑧 = τ𝜈 (1) 𝑥𝑧 . We thus obtain the following expression for the first-order correction of the velocity:

ũ(1) = 𝑅𝑒 ( 𝑧 2 2 -ℎ𝑧 ) ( 𝜆(ℎ -ℎ 𝑝 ) 𝜕ℎ 𝜕𝑡 + 𝜆 2 2 (ℎ -ℎ 𝑝 ) 3 𝜕ℎ 𝜕𝑥 ) + 𝑅𝑒 𝐹 𝑟 2 ( 𝑧 2 2 -ℎ𝑧 ) 𝜕ℎ 𝜕𝑥 + 𝑢 + (𝑥, 𝑡), (3.27)
where 𝑢 + is a matching term to be determined later. The first term on right side of (3.27) 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.20) derived from the classical constitutive law (2.3), the velocity derivative 𝜕 ũ(1) ∕𝜕𝑧 remains here bounded everywhere in the pseudo-plug zone. Moreover, the first-order correction ũ( 1) is here controlled by the viscous stress τ𝜈 [START_REF] Balmforth | Yielding to stress: recent developments in viscoplastic fluid mechanics[END_REF] 𝑥𝑧 and not by the yield-stress term τ𝑌 (0) 𝑥𝑧 . As a consequence, the expression (3.27) for ũ(1) includes the inertial terms, which was not the case in the previous approach.

In the sheared layer (𝑧 ⩽ ℎℎ 𝑝 )

The momentum equation along 𝑂𝑥 (2.18) at 𝑂(𝜀) writes here 𝐵𝑖 𝑅𝑒 𝜕𝜏 (1)

𝑥𝑧 𝜕𝑧 = 𝜕𝑢 (0) 𝜕𝑡 + 𝑢 (0) 𝜕𝑢 (0) 𝜕𝑥 + 𝑤 (0) 𝜕𝑢 (0) 𝜕𝑧 + 1 𝐹 𝑟 2 𝜕𝑝 (0) 𝜕𝑥 = 𝜆𝑧 𝜕ℎ 𝜕𝑥 + 𝜆 2 𝑧 2 2 (ℎ -ℎ 𝑝 ) 𝜕ℎ 𝜕𝑥 + 1 𝐹 𝑟 2 𝜕ℎ 𝜕𝑥 . (3.28)
Integration of this equation yields

𝜏 (1) 𝑥𝑧 = 𝑅𝑒 𝐵𝑖 ( 𝜆𝑧 2 2 𝜕ℎ 𝜕𝑡 + 𝜆 2 𝑧 3 6 (ℎ -ℎ 𝑝 ) 𝜕ℎ 𝜕𝑥 ) + 𝑧 𝑅𝑒 𝐵𝑖𝐹 𝑟 2 𝜕ℎ 𝜕𝑥 + 𝜏 (1) 𝑥𝑧 | 𝑧=0 , ( 3.29)
where the unknown term 𝜏 (1) 𝑥𝑧 | 𝑧=0 can be found from the stress continuity condition at the fake yield surface 𝜏 (1) 𝑥𝑧 | 𝑧=ℎ-ℎ 𝑝 = τ(1) 𝑥𝑧 | 𝑧=ℎ-ℎ 𝑝 . This leads to the following expression for the shear stress correction:

𝜏 (1) 𝑥𝑧 = 𝑅𝑒𝜆 2𝐵𝑖 (𝑧 2 -ℎ 2 + ℎ 2 𝑝 ) 𝜕ℎ 𝜕𝑡 + 𝑅𝑒𝜆 2 6𝐵𝑖 (ℎ -ℎ 𝑝 ) ( 𝑧 3 -(ℎ -ℎ 𝑝 ) 2 (ℎ + 2ℎ 𝑝 ) ) 𝜕ℎ 𝜕𝑥 + 𝑅𝑒 𝐵𝑖 𝐹 𝑟 2 (𝑧 -ℎ) 𝜕ℎ 𝜕𝑥 . (3.30)
By identification, the yield-stress contribution 𝜏 𝑌 (1) 𝑥𝑧 is zero, since all terms involve the factor 1∕𝐵𝑖. Then from (2.22) and the no-slip condition (2.24), the following velocity profile is obtained:

𝑢 (1) = 𝑅𝑒𝜆 2 𝑧 ( 𝑧 2 3 -ℎ 2 + ℎ 2 𝑝 ) 𝜕ℎ 𝜕𝑡 + 𝑅𝑒𝜆 2 6 (ℎ -ℎ 𝑝 )𝑧 ( 𝑧 3 4 -(ℎ -ℎ 𝑝 ) 2 (ℎ + 2ℎ 𝑝 ) ) 𝜕ℎ 𝜕𝑥 + + 𝑅𝑒 𝐹 𝑟 2 𝑧 ( 𝑧 2 -ℎ ) 𝜕ℎ 𝜕𝑥 . (3.31)

Matching

The matching term 𝑢 + in (3.27) is obtained from the continuity condition 𝑢 (1) | 𝑧=ℎ-ℎ 𝑝 = ũ(1) | 𝑧=ℎ-ℎ 𝑝 . Finally, the correction of the velocity profile in the pseudo-plug thus expresses as

ũ(1) = 𝑅𝑒𝜆 6 (ℎ -ℎ 𝑝 ) ( 3𝑧 2 -6ℎ𝑧 + (ℎ -ℎ 𝑝 ) 2 ) 𝜕ℎ 𝜕𝑡 + 𝑅𝑒𝜆 2 8 (ℎ -ℎ 𝑝 ) 3 ( 2𝑧 2 -4ℎ𝑧 + (ℎ -ℎ 𝑝 ) 2 ) 𝜕ℎ 𝜕𝑥 + + 𝑅𝑒 𝐹 𝑟 2 𝑧 ( 𝑧 2 -ℎ ) 𝜕ℎ 𝜕𝑥 . (3.32)

Velocity profile comparisons

In this section we compare the asymptotic expansion of the longitudinal velocity 𝑢 derived by using the classical expression of the yield-stress tensor (2.3) 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 𝑂(𝜀)-solution (blue dashed curve) and for the new 𝑂(𝜀)-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 𝑢 and 𝑤, pressure 𝑝 and stress components 𝜏 𝑖𝑗 ) are the same in both cases (see Figure 2a). At order 𝑂(𝜀), 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.26), 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 𝑂(𝜀) correction is expected to introduce a slight shearing in this zone. However, for the classical approach, the derivative (3.20) diverges at 𝑧 = ℎℎ 𝑝 , 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.

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 𝑂(𝜀) 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.17)-(2. [START_REF] Benney | Long waves on liquid films[END_REF], and that the boundary conditions (2.24)-(2.27) 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 𝐴 of the flow, let us define the depth-averaged value ⟨𝐴⟩ by

⟨𝐴⟩ = 1 ℎ ℎ ∫ 0 𝐴𝑑𝑧. (4.1)
For the averaged longitudinal velocity, it is convenient to use further the special notation ⟨𝑢⟩ = 𝑈 . In particular, the expression for the leading order term 𝑈 (0) is readily obtained from the equations (3.14) and (3.17):

𝑈 (0) = 𝜆ℎ 2 3 ( 1 - ℎ 𝑝 ℎ ) 2 ( 1 + ℎ 𝑝 2ℎ ) . (4.2)

Mass conservation

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

𝜕ℎ 𝜕𝑡 + 𝜕ℎ𝑈 𝜕𝑥 = 0. (4.3)
Introducing the leading-order asymptotic solution for 𝑈 , i.e. 𝑈 = 𝑈 (0) into (4.3) results in a kinematic wave equation:

𝜕ℎ 𝜕𝑡 + 𝜆ℎ(ℎ -ℎ 𝑝 ) 𝜕ℎ 𝜕𝑥 = O(𝜀), (4.4)
which indicate that a perturbation of depth is propagated at the speed 𝑐 0 = 𝜆ℎ(ℎℎ 𝑝 ). Equation (4.4) provides a useful expansion for 𝜕ℎ∕𝜕𝑡:

𝜕ℎ 𝜕𝑡 = -𝜆ℎ(ℎ -ℎ 𝑝 ) 𝜕ℎ 𝜕𝑥 + O(𝜀). (4.5)
In particular, this expansion allows us to write the first-order correction 𝑈 (1) in the following form:

3𝑈 (1) 𝑅𝑒 ℎ = [ 2 5 𝜆 2 ℎ 3 ( 1 - ℎ 𝑝 ℎ ) 2 ( 1 + ℎ 𝑝 ℎ + ℎ 2 𝑝 ℎ 2 - ℎ 3 𝑝 4ℎ 3 - ℎ 4 𝑝 4ℎ 4 ) - 1 𝐹 𝑟 2 ] ℎ 𝜕ℎ 𝜕𝑥 . (4.6)

Momentum balance

Averaging the momentum balance equation in the 𝑂𝑥-direction (2.18), together with the no-slip condition (2.24), the kinematic boundary condition (2.25), and the dynamic boundary conditions (2.26) and (2.27), provides

𝜕 𝜕𝑡 (ℎ𝑈 ) + 𝜕 𝜕𝑥 ⎛ ⎜ ⎜ ⎜ ⎝ ℎ⟨𝑢 2 ⟩ + 1 𝐹 𝑟 2 ℎ ∫ 0 𝑝𝑑𝑧 - 𝐵𝑖 𝑅𝑒 ℎ ∫ ℎ-ℎ 𝑝 𝜏 𝑥𝑥 𝑑𝑧 ⎞ ⎟ ⎟ ⎟ ⎠ = 𝜆ℎ -𝐵𝑖𝜏 𝑥𝑧 (0) 𝜀𝑅𝑒 , ( 4.7)
where the quantity 𝜏 𝑥𝑧 (0) denotes the shear stress at the bottom 𝑧 = 0. To express the integral of the pressure in (4.7) we use the equations (3.10) and (3.16):

1 𝐹 𝑟 2 ℎ ∫ 0 𝑝𝑑𝑧 = ℎ 2 2𝐹 𝑟 2 - 𝐵𝑖 𝑅𝑒 ℎ ∫ ℎ-ℎ 𝑝 𝜏 (0) 𝑥𝑥 𝑑𝑧 + 𝑂(𝜀). ( 4.8)
Taking into account the expression for the normal stress at leading-order (3.9), the integral term on the right-hand side of (4.8) writes

ℎ ∫ ℎ-ℎ 𝑝 𝜏 (0) 𝑥𝑥 𝑑𝑧 = 𝜋ℎ 𝑝 4 . (4.9)
As a result, the integral terms involving the normal stress 𝜏 (0) 𝑥𝑥 are constant and disappear from equation (4.7) after differentiation. Further, using the leading-order representation for the shear stress (3.4), the averaged momentum equation (4.7) reduces to

𝜕 𝜕𝑡 (ℎ𝑈 ) + 𝜕 𝜕𝑥 ( ℎ⟨𝑢 2 ⟩ + ℎ 2 2𝐹 𝑟 2 ) = - 𝐵𝑖𝜏 (1) 𝑥𝑧 (0) 𝑅𝑒 + 𝑂(𝜀 2 ). (4.10)
The quantity ⟨𝑢 2 ⟩ in (4.10) can be expressed by considering the velocity as the sum of its average value 𝑈 and a deviation 𝑢 * : 𝑢(𝑥, 𝑧, 𝑡) = 𝑈 (𝑥, 𝑡) + 𝑢 * (𝑥, 𝑧, 𝑡).

(4.11)

By definition ⟨𝑢 * ⟩ = 0, so that ⟨𝑢 2 ⟩ = 𝑈 2 + ⟨𝑢 * 2 ⟩. In equation (4.10), the deviation ⟨𝑢 * 2 ⟩ could be estimated at order 0 either as a function of ℎ or as a function of 𝑈 . This would lead to a closed system of two coupled equations for the flow height ℎ and the depth-averaged velocity 𝑈 [START_REF] Noble | Thin power-law film flow down an inclined plane: Consistent shallow-water models and stability under largescale perturbations[END_REF][START_REF] Fernandez-Nieto | Shallow Water equations for non-Newtonian fluids[END_REF]. However, as shown by Richard et al [START_REF] Richard | Vila Optimization of consistent two-equation models for thin film flows[END_REF] 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 [START_REF] Teshukov | Gas-dynamics analogy for vortex free-boundary flows[END_REF] and expanded by Richard and Gavrilyuk [START_REF] Richard | A new model of roll waves: comparison with Brock's experiments[END_REF], is thus to define an independent variable related to ⟨𝑢 * 2 ⟩. 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 ℎ, its average velocity 𝑈 and the variance of its velocity. In fact, it is more convenient to use

𝜑 = ⟨𝑢 * 2 ⟩ ℎ 2 (4.12)
as the third variable of the model, since it plays the role of an entropy for the system [START_REF] Richard | A new model of roll waves: comparison with Brock's experiments[END_REF]. In the particular case of a constant vorticity, the quantity 𝜑 is proportional to the square of the vorticity [START_REF] Teshukov | Gas-dynamics analogy for vortex free-boundary flows[END_REF]. For this reason, 𝜑 is called enstrophy. The expansion of this new variable 𝜑 = 𝜑 (0) + 𝜀𝜑 (1) + … (4.13) leads to 𝜑 (0) = ⟨ ( 𝑢 * (0) ) 2 ⟩ ∕ℎ 2 and 𝜑 (1) = 2⟨𝑢 * (0) 𝑢 * (1) ⟩∕ℎ 2 . The calculation at leading order gives

𝜑 (0) = 𝜆 2 ℎ 2 45 ( 1 - ℎ 𝑝 ℎ ) 5 ( 1 + 5ℎ 𝑝 4ℎ ) . (4.14)
For the expansion at order 𝑂(𝜀), we obtain

𝜑 (1) 𝑅𝑒 = 2𝜆ℎ 2 45 [ 3 7 𝜆 2 ℎ 3 ( 1 - ℎ 𝑝 ℎ ) 2 ( 1 + 101ℎ 𝑝 48ℎ + 29ℎ 2 𝑝 12ℎ 2 + 43ℎ 3 𝑝 48ℎ 3 - 7ℎ 4 𝑝 12ℎ 4
)

- 1 𝐹 𝑟 2 ( 1 + 9ℎ 𝑝 8ℎ + 3ℎ 2 𝑝 8ℎ 2 ) ] ( 1 - ℎ 𝑝 ℎ ) 3 𝜕ℎ 𝜕𝑥 . (4.15)
Finally, with the introduction of the enstrophy 𝜑, the equation (4.10) can be written

𝜕 𝜕𝑡 (ℎ𝑈 ) + 𝜕 𝜕𝑥 ( ℎ𝑈 2 + ℎ 3 𝜑 + ℎ 2 2𝐹 𝑟 2 ) = - 𝐵𝑖𝜏 (1) 𝑥𝑧 (0) 𝑅𝑒 . (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 𝜏 (1) 𝑥𝑧 (0) given by (3.30) reduces to

𝐵𝑖𝜏 (1) 𝑥𝑧 | 𝑧=0 𝑅𝑒 = [ 1 3 𝜆 2 ℎ 3 ( 1 - ℎ 𝑝 ℎ ) 2 ( 1 + ℎ 𝑝 ℎ + ℎ 2 𝑝 ℎ 2 ) - 1 𝐹 𝑟 2 ] ℎ 𝜕ℎ 𝜕𝑥 , ( 4.17)
which can also be rewritten in terms of 𝑈 (1) :

𝐵𝑖𝜏 (1) 𝑥𝑧 | 𝑧=0 𝑅𝑒 = 3𝑈 (1) 𝑅𝑒 ℎ - 𝜆 2 ℎ 4 15
( 1 - ℎ 𝑝 ℎ ) 3 ( 1 + 2ℎ 𝑝 ℎ + 3ℎ 2 𝑝 ℎ 2 + 3ℎ 3 𝑝 2ℎ 3 ) 𝜕ℎ 𝜕𝑥 . (4.18)
Note that the quantity 𝑈 (1) can be expressed as:

𝑈 (1) = 𝑈 -𝑈 (0) 𝜀 + 𝑂(𝜀) (4.19)
which has the structure of a relaxation term in 𝑈 . To express the derivative 𝜕ℎ∕𝜕𝑥 in (4.18), let us consider the following relation obtained from the expansions for 𝑈 and 𝜑:

𝜑 - 𝑈 2 5ℎ 2 ( 1 - ℎ 𝑝 ℎ ) ( 1 + 5ℎ 𝑝 4ℎ ) ( 1 + ℎ 𝑝 2ℎ ) -2 = 2𝜀 35 𝑅𝑒𝜆𝜑ℎ 3 ( 1 - ℎ 𝑝 ℎ ) 2 ( 1 + 5ℎ 𝑝 4ℎ ) -1 ( 1 + 61ℎ 𝑝 16ℎ + 63ℎ 2 𝑝 8ℎ 2 + 63ℎ 3 𝑝 8ℎ 3 + 21ℎ 4 𝑝 16ℎ 4 ) 𝜕ℎ 𝜕𝑥 + 𝜆ℎ 𝑝 20 ℎ ( 1 - ℎ 𝑝 ℎ ) 3 ( 1 + 5ℎ 𝑝 2ℎ + ℎ 2 𝑝 2ℎ 2 ) ( 1 + ℎ 𝑝 2ℎ ) -1 [ 𝑈 -𝑈 (0) ] + 𝑂(𝜀 2 ). (4.20)
The left-hand side of (4.20) has the structure of a relaxation term in 𝜑. Indeed, using expressions (4.2) and (4.14) we obtain

𝜑 (0) - ( 𝑈 (0) ) 2 5ℎ 2 ( 1 - ℎ 𝑝 ℎ ) ( 1 + 5ℎ 𝑝 4ℎ ) ( 1 + ℎ 𝑝 2ℎ ) -2 = 0. (4.21)
As a result, the derivative 𝜕ℎ∕𝜕𝑥 can be expressed as a sum of two relaxation terms for 𝑈 and 𝜑. Finally, the integrated momentum equation (4.16) can thus be written as

𝜕 𝜕𝑡 (ℎ𝑈 ) + 𝜕 𝜕𝑥 ( ℎ𝑈 2 + ℎ 3 𝜑 + ℎ 2 2𝐹 𝑟 2 ) = 1 𝜀𝑅𝑒 [ 𝜆ℎ -𝐵𝑖 - 3𝑈 ℎ𝛼 1 (𝜉) ] [ 𝛼 1 (𝜉) + 7 360
𝜆 2 ℎ 2 𝜑 𝛽 1 (𝜉) ] + 1 𝜀𝑅𝑒 7 6 𝜆ℎ 𝜑 [ 𝜑 - 𝑈 2 5ℎ 2 𝛼 2 (𝜉) ] 𝛽 2 (𝜉) (4.22)
where 𝜉 = ℎ 𝑝 ∕ℎ and the functions 𝛼 1 , 𝛼 2 , 𝛽 1 and 𝛽 2 are defined as: 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 𝑂𝑥-axis. The second term corresponds to a relaxation for the enstrophy.

𝛼 1 (𝜉) = (1 -𝜉) ( 1 + 𝜉 2 ) (4.23) 𝛼 2 (𝜉) = (1 -𝜉) ( 1 + 5 4 𝜉 ) ( 1 + 𝜉 2 ) -2 (4.24) 𝛽 1 (𝜉) = 𝜉 (1 -𝜉) 5 ( 1 + 5 4 𝜉 ) ( 1 + 5 2 𝜉 + 1 2 𝜉 2 ) ( 1 + 2𝜉 + 3𝜉 2 + 3 2 𝜉

Kinetic energy equation

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

𝜕 𝜕𝑡 ( 1 2 𝜌𝑣 2 ) + div ( 1 2 𝜌𝑣 2 𝒗 ) = div(𝝈𝒗) -𝝈 ∶ γ + 𝜌𝒈𝒗 (4.27)
where we recall that 𝒗 denotes the velocity field and 𝝈 = -𝑝𝑰 + 𝝉. Introducing the components of the vectors and tensors involved in (4.27), we get the following expression written in dimensionless form:

𝜕 𝜕𝑡 ( 𝑢 2 2 + 𝜀 2 𝑤 2 2 ) + 𝜕 𝜕𝑥 [ 𝑢 ( 𝑢 2 2 + 𝜀 2 𝑤 2 2 - 𝑥 tan 𝜃 𝜀𝐹 𝑟 2 + 𝑧 𝐹 𝑟 2 ) + 𝑝𝑢 𝐹 𝑟 2 - 𝐵𝑖 𝑅𝑒 (𝜏 𝑥𝑥 𝑢 + 𝜀𝜏 𝑥𝑧 𝑤) ] + 𝜕 𝜕𝑧 [ 𝑤 ( 𝑢 2 2 + 𝜀 2 𝑤 2 2 - 𝑥 tan 𝜃 𝜀𝐹 𝑟 2 + 𝑧 𝐹 𝑟 2 ) + 𝑝𝑤 𝐹 𝑟 2 - 𝐵𝑖 𝜀𝑅𝑒 𝜏 𝑥𝑧 𝑢 - 𝐵𝑖 𝑅𝑒 𝜏 𝑧𝑧 𝑤 ] = - 2𝐵𝑖 𝑅𝑒 ( 𝜏 𝑝 𝑥𝑥 + 2 𝜀 𝐵𝑖 𝜕𝑢 𝜕𝑥 ) 𝜕𝑢 𝜕𝑥 - 𝐵𝑖 𝜀𝑅𝑒 ( 𝜏 𝑥𝑧 + 1 𝐵𝑖 ( 𝜕𝑢 𝜕𝑧 + 𝜀 2 𝜕𝑤 𝜕𝑧 )) ( 𝜕𝑢 𝜕𝑧 + 𝜀 2 𝜕𝑤 𝜕𝑧 ) . (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:

𝜕 𝜕𝑡 ( ℎ𝑈 2 2 + ℎ 3 𝜑 2 + ℎ 2 2𝐹 𝑟 2 ) + 𝜕 𝜕𝑥 ( ℎ⟨𝑢 3 ⟩ 2 + ℎ 2 𝑈 𝐹 𝑟 2 ) = - 𝜆ℎ𝑈 (1) 𝑅𝑒 + 𝐵𝑖 𝑅𝑒 ℎ ∫ 0 𝜏 𝑌 (0) 𝑥𝑧 𝜕𝑢 (1) 𝜕𝑧 𝑑𝑧. (4.29)
From the equation (4.11) and the definition of enstrophy 𝜑 , we have ⟨𝑢 2 ⟩ = 𝑈 2 + ℎ 2 𝜑 and

⟨𝑢 3 ⟩ = 𝑈 3 + 3ℎ 2 𝑈 𝜑 + ⟨𝑢 * 3 ⟩. (4.30)
Accordingly, (4.29) can be rewritten as:

𝜕 𝜕𝑡 ( ℎ𝑈 2 2 + ℎ 3 𝜑 2 + ℎ 2 2𝐹 𝑟 2 ) + 𝜕 𝜕𝑥 ( ℎ𝑈 3 2 + 3ℎ 3 𝑈 𝜑 2 + ℎ 2 𝑈 𝐹 𝑟 2 ) = = - 𝜆ℎ𝑈 (1) 𝑅𝑒 + 𝐵𝑖 𝑅𝑒 ℎ ∫ 0 𝜏 𝑌 (0) 𝑥𝑧 𝜕𝑢 (1) 𝜕𝑧 𝑑𝑧 - 𝜕 𝜕𝑥 ( ℎ⟨𝑢 * 3 ⟩ 2 ) . (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 water flows, Richard and Gavrilyuk [START_REF] Richard | A new model of roll waves: comparison with Brock's experiments[END_REF] considered a weakly-sheared flow assumption, which allowed them to neglect the term involving the cubic deviation ⟨𝑢 * 3

⟩ 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 ⟨𝑢 * 3 ⟩ as a fourth variable would lead to an infinite hierarchy of equations. To close the system, we thus calculate ⟨𝑢 * 3

⟩ as a function of ℎ by employing the asymptotic expansions derived above. Namely, using expression (4.30), we obtain the following expression at leading order:

⟨𝑢 * 3 ⟩ = - 2 945 𝜆 3 ℎ 6 ( 1 - ℎ 𝑝 ℎ ) 7 ( 1 + 49ℎ 𝑝 16ℎ + 35ℎ 2 𝑝 8ℎ 2 ) + 𝑂(𝜀). (4.32)
The derivative of ⟨𝑢 * 3 ⟩ involved in (4.31) has thus the form:

- 𝜕 𝜕𝑥 ( ℎ⟨𝑢 * 3 ⟩ 2 ) = 𝜆 3 ℎ 6 135 ( 1 - ℎ 𝑝 ℎ ) 6 ( 1 + 21ℎ 𝑝 8ℎ + 57ℎ 2 𝑝 16ℎ 2 + 5ℎ 3 𝑝 4ℎ 3 ) 𝜕ℎ 𝜕𝑥 + 𝑂(𝜀). (4.33)
The remaining term to be calculated in (4.31) is the integral:

𝐵𝑖 𝑅𝑒 ℎ ∫ 0 𝜏 𝑌 (0) 𝑥𝑧 𝜕𝑢 (1) 𝜕𝑧 𝑑𝑧 = 𝜆ℎ 𝑝 2 [ 5 12 𝜆 2 ℎ 5 ( 1 - ℎ 𝑝 ℎ ) 2 ( 1 + ℎ 𝑝 ℎ + 3ℎ 2 𝑝 5ℎ 2 - ℎ 3 𝑝 ℎ 3 ) -ℎ 2 𝐹 𝑟 -2 ( 1 - ℎ 2 𝑝 3ℎ 2
) ] 𝜕ℎ 𝜕𝑥 (4.34)

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

𝜕 𝜕𝑡 ( ℎ𝑈 2 2 + ℎ 3 𝜑 2 + ℎ 2 2𝐹 𝑟 2 ) + 𝜕 𝜕𝑥 ( ℎ𝑈 3 2 + 3ℎ 3 𝑈 𝜑 2 + ℎ 2 𝑈 𝐹 𝑟 2 ) = 𝑈 𝜀𝑅𝑒 [ 𝜆ℎ -𝐵𝑖 - 3𝑈 ℎ𝛼 1 (𝜉) ] [ 𝛼 1 (𝜉) + 7 1080 𝜆 2 ℎ 2 𝜑
(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 𝑂𝑥 axis and the power of the yield stress and the viscous friction forces. The second term corresponds to the relaxation for the enstrophy 𝜑.

Equation of enstrophy

From equations (4.3), (4.22) and (4.35), the following evolution equation for the enstrophy 𝜑 can be derived:

ℎ 2 2 ( 𝜕ℎ𝜑 𝜕𝑡 + 𝜕ℎ𝑈 𝜑 𝜕𝑥 ) = 𝑈 𝜀𝑅𝑒 [ 𝜆ℎ -𝐵𝑖 - 3𝑈 ℎ𝛼 1 (𝜉) ] [ 7 360 𝜆 2 ℎ 2 𝜑 𝛽 1 (𝜉) ( 𝑟(𝜉) 3 - 1
)]

+ 𝑈 𝜀𝑅𝑒 7 6 𝜆ℎ 𝜑 [ 𝜑 - 𝑈 2 5ℎ 2 𝛼 2 (𝜉) ] 𝛽 2 (𝜉) ( 𝑟(𝜉) 3 - 1
) . (4.37)

Note that, since 𝑟(𝜉)∕3 -1 < 0, the enstrophy properly relaxes toward its equilibrium value 𝛼 2 (𝜉)𝑈 2 ∕5ℎ 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 enstrophy, which plays the role of an entropy of the system, should increase according to a Rankine-Hugoniot relation (for more details see [START_REF] Richard | A new model of roll waves: comparison with Brock's experiments[END_REF]). Hence, we shall only consider the system consisting of equations (4.3), (4.22), (4.35) in what follows.

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 [START_REF] Richard | A three-equation model for thin films down an inclined plane[END_REF] in that the enstrophy relaxes toward a more physical value 𝑈 2 ∕5ℎ 2 .

Π = ℎ 3 𝜑 + 𝑔ℎ 2 cos 𝜃 2 , ( 4

Reconstruction of the velocity field

The solution of the system (4.40)-(4.42) provides us with the free-surface ℎ, the averaged velocity 𝑈 and the enstrophy 𝜑 at order 𝑂(𝜀). For a comprehensive analysis of the flows, it can also be useful to reconstruct the velocity field (𝑢, 𝑤) from these computed variables ℎ, 𝑈 and 𝜑. Let us recall that, unlike the previous studies based on the classical expression (2.3) for the yield-stress tensor, a notable benefit of our model is the possibility to consistently reconstruct velocity profiles that are smooth at order 𝑂(𝜀) (see Section 3.3).

At leading order, the longitudinal velocity 𝑢 can be expressed as

𝑢 (0) = 𝑈 (0) 𝑓 𝑠ℎ (𝑧) for 𝑧 < ℎ -ℎ 𝑝 , (5.1) ũ(0) = 𝑈 (0) 𝑓 𝑝𝑙 for 𝑧 ⩾ ℎ -ℎ 𝑝 , ( 5.2)
with

𝑓 𝑠ℎ (𝑧) = 3𝑧 ℎ ( 1 - 𝑧 2(ℎ -ℎ 𝑝 ) ) ( 1 - ℎ 𝑝 ℎ ) -1 ( 1 + ℎ 𝑝 2ℎ ) -1 , 𝑓 𝑝𝑙 = 3 2 ( 1 + ℎ 𝑝 2ℎ ) -1 . (5.3)
Note that since the quantity 𝑈 (0) is the leading-order solution for the average velocity, we have

1 ℎ ∫ ℎ-ℎ 𝑝 0 𝑓 𝑠ℎ (𝑧) d𝑧 + 1 ℎ ∫ ℎ ℎ-ℎ 𝑝 𝑓 𝑝𝑙 d𝑧 = 1. (5.4)
At order 𝑂(𝜀), the expansions for the longitudinal velocity are written (1) -𝑈 (1) 𝑓 𝑠ℎ ] + 𝑂(𝜀 2 ) for 𝑧 < ℎℎ 𝑝 , (5.5)

𝑢 = 𝑈 𝑓 𝑠ℎ + 𝜀 [ 𝑢
𝑢 = 𝑈 𝑓 𝑝𝑙 + 𝜀 [ ũ(1)
-𝑈 (1) 𝑓 𝑝𝑙 ] + 𝑂(𝜀 2 ) for 𝑧 ⩾ ℎℎ 𝑝 , (5.6) where 𝑈 𝑓 𝑠ℎ and 𝑈 𝑓 𝑝𝑙 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 ℎ, 𝑈 and 𝜑, excluding their derivatives, since the derivatives 𝜕ℎ∕𝜕𝑡 and 𝜕ℎ∕𝜕𝑥, which are involved in the first-order corrections 𝑢 1 and ũ1 (see equations (3.31)

] × × ( 7 24𝜑ℎ 2 [ 𝜆ℎ -𝐵𝑖 - 3𝑈 ℎ𝛼 1 (𝜉) ] 𝜁 1 (𝜉) + 7 18
𝜆 𝜑 2 ℎ [ 𝜑 - 𝑈 2 5ℎ 2 𝛼 2 (𝜉) ] 𝜁 2 (𝜉) ) (5.7)
with 𝜉 = ℎ 𝑝 ∕ℎ and 𝜂 = 𝑧∕ℎ. As above, the function 𝑓 , which corresponds to 𝑓 𝑠ℎ and 𝑓 𝑝𝑙 , gives the profile of order 0 and the other terms give the corrections of order 1. For the sheared layer, 𝑧 < ℎℎ 𝑝 , the functions 𝑓 , 𝜒 and 𝜓 are defined by

𝑓 (𝜂; 𝜉) = 3𝜂 (1 -𝜉) 2 ( 1 -𝜉 - 𝜂 2 ) ( 1 + 1 2 𝜉 ) -1 , (5.8) 𝜒(𝜂; 𝜉) = 𝜉𝜂 (1 -𝜉) 2 ( 1 - 3 2 𝜂 -𝜉 2 + 1 2 𝜂𝜉 2 ) ( 1 + 1 2 𝜉 ) -1 (
) (1 -𝜉) 2 ( 1 + 1 2 𝜉 ) -1 .
(5.13)

The functions 𝜁 1 and 𝜁 2 are given by Regarding the normal velocity 𝑤, the derivative 𝜕ℎ∕𝜕𝑥 in (3.15) can be expressed consistently in a similar way. Since 𝑤 is a first-order quantity compared to 𝑢, it is sufficient to express 𝑤 0 as a sum of relaxation terms to obtain an accuracy at order 1. In the end, the following expression is obtained:

𝜁 1 (𝜉) = 𝜉 (1 -𝜉) 2 ( 1 + 5 4 𝜉 ) ( 1 + 5 2 𝜉 + 1 2 𝜉 2 ) ( 1
𝑤 = 𝜆ℎ 2 𝜀𝑅𝑒 ( 7 24𝜑ℎ 2 [ 𝜆ℎ -𝐵𝑖 - 3𝑈 ℎ𝛼 1 (𝜉) ] 𝜁 1 (𝜉) + 7 18 𝜆 𝜑 2 ℎ [ 𝜑 - 𝑈 2 5ℎ 2 𝛼 2 (𝜉) ] 𝜁 2 (𝜉) ) 𝑞(𝜂; 𝜉) (5.16)
The function 𝑞 for the sheared layer, 𝑧 < ℎℎ 𝑝 , is given by

𝑞(𝜂; 𝜉) = - 𝜂 2 2 ,
(5.17)

and for the pseudo-plug, 𝑧 ⩾ ℎℎ 𝑝 , by

𝑞(𝜂; 𝜉) = -(1 -𝜉) ( 𝜂 - 1 -𝜉 2 ) .
(5.18)

In dimensional variables the reconstructed longitudinal velocity expresses as

𝑢 = 𝑈 𝑓 (𝜂; 𝜉) + [ 𝜒(𝜂; 𝜉) 𝑔ℎ 2 cos 𝜃 2 + 𝜓(𝜂; 𝜉) 𝑔 2 ℎ 5 sin 2 𝜃 15𝜈 2 ] × ( 7 24𝜈𝜑ℎ 2 [ 𝑔ℎ sin 𝜃 - 𝜏 𝑐 𝜌 - 3𝜈𝑈 ℎ𝛼 1 (𝜉) ] 𝜁 1 (𝜉) + 7 18
𝑔 sin 𝜃 𝜈𝜑 2 ℎ [ 𝜑 - 𝑈 2 5ℎ 2 𝛼 2 (𝜉) ] 𝜁 2 (𝜉) ) (5.19)
with 𝜉 = 𝜏 𝑐 ∕(𝜌𝑔ℎ sin 𝜃). The dimensional normal velocity expresses as

𝑤 = 𝑔ℎ 2 sin 𝜃 𝜈 ( 7 24𝜑ℎ 2 [ 𝑔ℎ sin 𝜃 - 𝜏 𝑐 𝜌 - 3𝜈𝑈 ℎ𝛼 1 ( ξ) ] 𝜁 1 (𝜉) + 7 18
𝑔 sin 𝜃 𝜈𝜑 2 ℎ [ 𝜑 - 𝑈 2 5ℎ 2 𝛼 2 (𝜉) ] 𝜁 2 (𝜉)
) 𝑞(𝑧∕ℎ; 𝜉). (5.20)

Applications

Characteristic velocities

The left-hand side of the system (4.40)-(4.42) is identical as in Richard and Gavrilyuk [START_REF] Richard | A new model of roll waves: comparison with Brock's experiments[END_REF]. The right-hand side involves only relaxation terms. Accordingly, the model is hyperbolic as proven by Teshukov [START_REF] Teshukov | Gas-dynamics analogy for vortex free-boundary flows[END_REF]. The system can be rewritten in the matrix form: where 𝑽 = [ℎ, 𝑈 , 𝜑] 𝑇 , 𝑺 is the matrix of the source terms and the matrix 𝑨 is given by

𝜕𝑽 𝜕𝑡 + 𝑨 𝜕𝑽 𝜕𝑥 = 𝑺, ( 6.1)
𝑨 = ⎡ ⎢ ⎢ ⎣ 𝑈 ℎ 0 𝑔 cos 𝜃 + 3ℎ𝜑 𝑈 ℎ 2 0 0 𝑈 ⎤ ⎥ ⎥ ⎦ (6.2)
The three characteristic velocities of the model are given by the eigenvalues of 𝑨, namely:

𝑈 , 𝑈 - √ 𝑔ℎ cos 𝜃 + 3ℎ 2 𝜑, 𝑈 + √ 𝑔ℎ cos 𝜃 + 3ℎ 2 𝜑 ( 6.3)
As clear from these expressions, the shearing effect contributes to the characteristic velocities through the term 3ℎ 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 yieldstress 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 ℎ = 1 + ℎ ′ , 𝑈 = 𝑈 (0) + 𝑈 ′ and 𝜑 = 𝜑 (0) + 𝜑 ′ , where ℎ ′ , 𝑈 ′ and 𝜑 ′ are small sinusoidal perturbations. Namely, we take the perturbations of the form

[ℎ ′ , 𝑈 ′ , 𝜑 ′ ] 𝑇 = [𝐴 1 , 𝐴 2 , 𝐴 3 ] 𝑇 exp[i𝑘(𝑥 -𝑐𝑡)],
where 𝑘 is the wavenumber and 𝑐 is the phase velocity. The dispersion relation is found by equating the determinant of the linearised system to zero. The details of this lengthy calculation are given in Appendix C. Up to the first order in 𝜀, the relation writes

𝑐 = 𝜆(1 -ℎ 𝑝 ) + i𝑘𝜀 𝑅𝑒 3 [ 2𝜆 2 5 (1 -ℎ 𝑝 ) 2 ( 1 + ℎ 𝑝 + ℎ 2 𝑝 - ℎ 3 𝑝 4 - ℎ 4 𝑝 4 ) - 1 𝐹 𝑟 2 ] + 𝑂(𝜀 2 ). (6.4)
The base flow is stable if Im(𝑐) < 0. Balmforth and Liu [START_REF] Balmforth | Roll waves in mud[END_REF] studied a linear stability of the base flow for the linearized Cauchy equations (known as generalized Orr-Sommerfeld equations) in the case of a Herschel-Bulkley fluid with power flow index 𝑛. In order to compare our result with one obtained by these authors in particular case of the Bingham fluid (𝑛 = 1), we choose the same characteristic velocity, namely 𝑢 0 = 𝑔ℎ 2 0 sin 𝜃∕𝜈. Note that this choice imposes 𝜆 = 1. As a result, in the long-wave limit (𝜀 → 0) stability occurs for 𝑅𝑒 < 𝑅𝑒 𝑐 , where the critical Reynolds number 𝑅𝑒 𝑐 is given by

𝑅𝑒 𝑐 = 10 cot 𝜃 (1 -𝐵𝑖) 2 ( 4 + 4𝐵𝑖 + 4𝐵𝑖 2 -𝐵𝑖 3 -𝐵𝑖 4 ) -1 . (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 [START_REF] Balmforth | Roll waves in mud[END_REF][START_REF] Fernandez-Nieto | Shallow Water equations for non-Newtonian fluids[END_REF]. Figure 3 (a) compares our condition (6.5) with the one obtained by Balmforth and Liu [START_REF] Balmforth | Roll waves in mud[END_REF] from generalized Orr-Sommerfeld equations in the particular case of the Bingham fluid. As shown, our model reproduces almost the same instability threshold. The slight difference (almost not visible in the figure) can be attributed to the different formulation of the constitutive law. Note that the recent experimental study on the instability of viscoplastic fluids made by Mounkaila Noma et al [START_REF] Mounkaila Noma | Primary instability of a visco-plastic film down an inclined plane: experimental study[END_REF] shows good agreement with the criteria obtained by Balmforth and Liu [START_REF] Balmforth | Roll waves in mud[END_REF] in the general case of a Herschel-Bulkley fluid.

0 0.5 1 1.5 2

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 [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF]). 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 ℎ, 𝑈 and 𝜑 are imposed, while at the outlet a Neumann boundary condition is considered. Parameters are chosen so that 𝑅𝑒 > 𝑅𝑒 𝑐 , allowing an instability to develop. Namely, we choose the following values for the dimensional parameters: 𝐾 = 20 Pa s, 𝜌 = 1000 kg m -3 , ℎ 0 = 0.01 m, 𝜃 = 18 • and 𝑢 0 = 𝑔ℎ 2 0 sin 𝜃∕𝜈 = 1.5 m s -1 . This corresponds to the dimensionless parameters: 𝑅𝑒 = 75.71, 𝐹 𝑟 = 4.96 and 𝜆 = 1. In order to capture the influence of the yield stress on the waves, we consider two cases, namely 𝜏 𝑐 = 3 Pa and 𝜏 𝑐 = 20 Pa, corresponding to 𝐵𝑖 = 0.1 and 𝐵𝑖 = 0.6 respectively. Figure 4 shows the evolution of the simulated free-surface height ℎ and flow enstrophy 𝜑. Typical roll waves, with a discontinuous shock at the front, develop for both values of 𝐵𝑖. 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 enstrophy, 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 enstrophy 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 [START_REF] Balmforth | Roll waves in mud[END_REF]). A similar behaviour was reported for Newtonian fluids by [START_REF] Richard | A new model of roll waves: comparison with Brock's experiments[END_REF], and shown to be in good agreement with experimental data. It can also be noted that the enstrophy is strongly influenced by the shock, with a marked peak at the front of the wave (Figure 5b).

Interpretation of enstrophy variations

In this section, we relate the enstrophy 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 enstrophy from its equilibrium value 𝜑 𝑒𝑞 given by Eq. (4.21):

𝜑 𝑒𝑞 = 𝑈 2 5ℎ 2 ( 1 - 𝜏 𝑐 𝜌𝑔ℎ sin 𝜃 ) ( 1 + 5𝜏 𝑐 4𝜌𝑔ℎ sin 𝜃 ) ( 1 + 𝜏 𝑐 2𝜌𝑔ℎ sin 𝜃 ) -2 . ( 6.6)
As a reference, let us start with the Newtonian case (𝐵𝑖 = 0). Alekseenko and Nakoryakov [START_REF] Alekseenko | Wave formation on vertical falling liquid films[END_REF] 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 [START_REF] Denner | Solitary waves on falling liquid films in the inertia-dominated regime[END_REF] from experimental data and DNS. As demonstrated by Richard et al [START_REF] Richard | A three-equation model for thin films down an inclined plane[END_REF], these different zones can be directly related to the deviation of the enstrophy from its equilibrium value. In the Newtonian case, 𝜑 𝑒𝑞 = 𝑈 2 ∕5ℎ 2 , and the longitudinal velocity 𝑢 at order 𝑂(𝜀) is given by (see Eq. (5.19)):

𝑢 = 𝑈 𝑓 (𝑧∕ℎ; 0) + 7 270 ( 𝑔 sin 𝜃 𝜈 ) 3 ℎ 4 𝜑 2 [ 𝜑 - 𝑈 2 5ℎ 2 ] 𝜓(𝑧∕ℎ; 0). (6.7)
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 𝜑 = 𝜑 𝑒𝑞 . Hence, for 𝜑 ≈ 𝜑 𝑒𝑞 , a parabolic velocity profile is recovered, corresponding to region I. For 𝜑 > 𝜑 𝑒𝑞 , the velocity profile is "less filled" and corresponds to region II, while for 𝜑 < 𝜑 𝑒𝑞 , the velocity profile is "more filled" and corresponds to region III. Comparisons between 𝜑 and 𝜑 𝑒𝑞 in roll waves simulated with our model for 𝐵𝑖 = 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 𝜑 𝑒𝑞 (Figure 6a), and the corresponding differences in the shape of the velocity profiles, are in good agreement with the experimental observations of Alekseenko and Nakoryakov [START_REF] Alekseenko | Wave formation on vertical falling liquid films[END_REF]. 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 enstrophy at the shock and, as shown by Richard et al [START_REF] Richard | A three-equation model for thin films down an inclined plane[END_REF], can be alleviated by adding a diffusion term to the model. Indeed, the shock and associated enstrophy discontinuity are suppressed by the additional diffusion, which leads to an extension of region III up to the peak of the wave. Turning to the case of a Bingham fluid, the equilibrium enstrophy is obtained from Eq. (4.21), 𝜑 𝑒𝑞 = 𝛼 2 (𝜉)𝑈 2 ∕5ℎ 2 , and the longitudinal velocity at order 𝑂(𝜀) 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 enstrophy 𝜑 from the equilibrium value 𝜑 𝑒𝑞 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 𝑈 from the equilibrium value 𝑈 𝑒𝑞 = 𝛼 1 (𝜉)(𝜌𝑔ℎ 2 sin 𝜃 -𝜏 𝑐 ℎ)∕3𝜌𝜈. Nevertheless, we observed that the regions where 𝑈 < 𝑈 𝑒𝑞 (resp. 𝑈 > 𝑈 𝑒𝑞 ) approximately correspond to regions where 𝜑 < 𝜑 𝑒𝑞 and (resp. 𝜑 > 𝜑 𝑒𝑞 ). Again, comparisons between 𝜑 and 𝜑 𝑒𝑞 along a simulated roll wave and corresponding reconstructed velocity profiles are shown in Figures 8 and9. Globally, similar trends as for Newtonian fluids are recovered. In zones for which 𝜑 ≈ 𝜑 𝑒𝑞 , a parabolic velocity profile with an unsheared plug zone is observed (regions I). In zones for which 𝜑 > 𝜑 𝑒𝑞 , the velocity profile is "less filled" than the equilibrium profile (region II), while for 𝜑 > 𝜑 𝑒𝑞 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 enstrophy near the front of the wave, and increase the extent of region III up to the front of the wave.

Conclusion

In this paper, a three-equation shallow-flow model for a Bingham fluid propagating down an inclined 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 𝑂(𝜀), close to the free surface. In contrast to previous approaches, the expansion is constructed by relaxing the classical assumption of alignment between the yield-stress tensor and the strain-rate. As a consequence, shearing in the pseudoplug 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 𝑂(𝜀). In addition, this asymptotic solution accounts for inertial terms at order 𝑂(𝜀) 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 enstrophy, 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 Nieto et al [START_REF] Fernandez-Nieto | Shallow Water equations for non-Newtonian fluids[END_REF], 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 𝑂(𝜀).

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 in good agreement with the result of Balmforth and Liu [START_REF] Balmforth | Roll waves in mud[END_REF] obtained from the generalized Orr-Sommerfeld equations for viscoplastic fluids. The model is 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 [START_REF] Piau | Carbopol gels: elastoviscoplastic and slippery glasses made of indi-vidual swollen sponges. Meso-and macroscopic properties, constitutive equa-tions and scaling laws[END_REF][START_REF] Thompson | The yield stress tensor[END_REF][START_REF] De Cagny | The yield normal stress[END_REF] 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. The difference with the Newtonian model derived by Richard et al [START_REF] Richard | A three-equation model for thin films down an inclined plane[END_REF] is that the enstrophy relaxes here toward 𝑈 2 ∕5ℎ 2 , while it relaxes toward the term 𝑔 sin 𝜃ℎ 2 ∕45𝜈 2 in this former model. The two expressions are equivalent except for modelling rest states (𝑈 = 0). In particular, the expression used by Richard et al [START_REF] Richard | A three-equation model for thin films down an inclined plane[END_REF] is at the origin of nonphysical sources of momentum and energy at rest, which is not the case with the present model.

Appendix A. Derivation of the energy equation

Appendix C. Derivation of the linear instability criterion

We linearize the system of equations

⎤ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎣ ℎ ′ 𝑈 ′ 𝜑 ′ ⎤ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎣ 0 0 0 ⎤ ⎥ ⎥ ⎦ (C.1)
The coefficients 𝐿 𝑚 𝑗 (𝑗 = ℎ, 𝑈 , 𝜑) come from the left-hand side of the momentum equation (4.22), and are given by:

𝐿 𝑚 ℎ = i𝜀𝑘 ( 3𝜑 (

Figure 1 :

 1 Figure 1: Definition sketch

. 22 )

 22 while the first dynamic boundary condition(2.26) at O(𝜀) is τ(1) 𝑥𝑧 | 𝑧=ℎ(𝑥) = 2𝛿 𝜕ℎ 𝜕𝑥 . (3.23) Integration of equation (3.21) coupled with (3.23) thus leads to the following expression for the shear stress correction:

Figure 2 :

 2 Figure 2: Velocity profiles for the leading order solution (black dashed-dot curve), for the classical 𝑂(𝜀)-solution (dashed blue curve) and for the new 𝑂(𝜀)-solution (red curve): (a) dimensionless velocity profiles obtained with ℎ = 1, 𝜕ℎ∕𝜕𝑥 = -0.05, 𝜕ℎ∕𝜕𝑡 = 0.05, 𝑅𝑒 = 1, 𝐹 𝑟 = 0.6, 𝐵𝑖 = 0.3 and 𝜃 = 20 • ; (b) velocity profiles normalized by the depth-averaged value 𝑈 , obtained with ℎ = 1, 𝜕ℎ∕𝜕𝑥 = -0.1, 𝜕ℎ∕𝜕𝑡 = 0.07, 𝑅𝑒 = 0.1, 𝐹 𝑟 = 0.18, 𝐵𝑖 = 0.3 and 𝜃 = 18 • .

Figure 3 :

 3 Figure 3: Comparison between the instability threshold derived with our new model and that obtained by Balmforth and Liu [11]: critical Reynolds number 𝑅𝑒 𝑐 as a function of the Bingham number 𝐵𝑖.

Figure 4 :

 4 Figure 4: Development of roll waves as a result of uniform flow instability for two values of 𝐵𝑖: evolution of (a) depth ℎ and (b) enstrophy 𝜑 as functions of the distance from the system entrance (𝑅𝑒 = 75.71, 𝐹 𝑟 = 4.96, 𝜆 = 1, and 𝜃 = 18 • ).

Figure 5 :

 5 Figure 5: Shape of the roll waves: evolution depth ℎ and enstrophy 𝜑 (𝑅𝑒 = 75.71, 𝐹 𝑟 = 4.96, 𝐵𝑖 = 0.1, 𝜆 = 1 and 𝜃 = 18 • ).

Figure 6 :Figure 7 :

 67 Figure 6: (a) Different regions within a roll wave defined from the deviation between the enstrophy 𝜑 and its equilibrium value 𝜑 𝑒𝑞 : Newtonian case (𝑅𝑒 = 75.71, 𝐹 𝑟 = 4.96, 𝐵𝑖 = 0, 𝜆 = 1 and 𝜃 = 18 • ). Values of 𝜑 𝑒𝑞 are computed from the local depth ℎ and average velocity 𝑈 along the wave (see text). (b) Corresponding free-surface profile.

Figure 8 :Figure 9 :

 89 Figure 8: (a) Different zones within a roll wave defined from the deviation between the enstrophy 𝜑 and its equilibrium value 𝜑 𝑒𝑞 : Bingham case (𝑅𝑒 = 75.71, 𝐹 𝑟 = 4.96, 𝐵𝑖 = 0.1, 𝜆 = 1, ℎ 𝑝 = 0.1 and 𝜃 = 18 • ). Values of 𝜑 𝑒𝑞 are computed from the local depth ℎ and average velocity 𝑈 along the wave (see text). (b) Corresponding free-surface profile.

  (4.3),(4.22) and(4.35) by considering small sinusoidal perturbations around the base flow (4.2) and (4.14): ℎ = 1 + ℎ ′ , 𝑈 = 𝑈 (0) + 𝑈 ′ and 𝜑 = 𝜑 (0) + 𝜑 ′ , with [ℎ ′ , 𝑈 ′ , 𝜑 ′ ] 𝑇 = [𝐴 1 , 𝐴 2 , 𝐴 3 ] 𝑇 exp[i𝑘(𝑥 -𝑐𝑡)].The three linearized equations can be written as:

  𝐾∕𝜌, 𝜉 = 𝜏 𝑐 ∕(𝜌𝑔ℎ sin 𝜃), and the functions 𝑟, 𝛼 𝑖 , 𝛽 𝑖 (𝑖 = 1, 2) defined in (4.23)-(4.26) and (4.36). Hence, the

																		.38)
	𝑒 =	1 2	(	𝑈 2 + ℎ 2 𝜑 + 𝑔ℎ cos 𝜃	)	.		(4.39)
	Equations (4.3), (4.22) and (4.35) can then be written as (in dimensional form):
	𝜕ℎ 𝜕𝑡		+	𝜕ℎ𝑈 𝜕𝑥	= 0,								(4.40)
	𝜕 𝜕𝑡	(ℎ𝑈 ) +	𝜕 𝜕𝑥	( ℎ𝑈 2 + Π	)	=	[ 𝑔ℎ sin 𝜃 -	𝜏 𝑐 𝜌	-3𝜈	𝑈 ℎ𝛼 1 (𝜉)	] [ 𝛼 1 (𝜉) +	7 360	𝜈 2 𝜑 (𝑔ℎ sin 𝜃) 2	𝛽 1 (𝜉)	]
																		+	7 6	𝑔ℎ sin 𝜃 𝜑	[	𝜑 -	𝑈 2 5ℎ 2	] 𝛼 2 (𝜉)	𝛽 2 (𝜉), (4.41)
	𝜕 𝜕𝑡	(ℎ𝑒) +	𝜕 𝜕𝑥	[	ℎ𝑈	( 𝑒 +	Π ℎ	)]	= 𝑈	[ 𝑔ℎ sin 𝜃 -	𝜏 𝑐 𝜌	-3𝜈	𝑈 ℎ𝛼 1 (𝜉)	] [ 𝛼 1 (𝜉) +	7 1080	𝜈 2 𝜑 (𝑔ℎ sin 𝜃) 2	𝛽 1 (𝜉)𝑟(𝜉) ]
																		+	7 18	𝑈 𝑔ℎ sin 𝜃 𝜑	[ 𝜑 -	𝑈 2 5ℎ 2	] 𝛼 2 (𝜉) 𝛽 2 (𝜉)𝑟(𝜉) (4.42)
	with 𝜈 =

left-hand side of system (4.40)-(4.42) indeed has the form of Euler equations for compressible fluids with Π, 𝑒 and 𝜑 playing the roles of pressure, energy and enstrophy, respectively.

The reduced system obtained for the Newtonian case (𝜏 𝑐 = 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

  Starting from the dimensionless form of the work-energy theorem (4.28) and dropping terms of order 𝑂(𝜀 2 ) or smaller, one obtains Note that the integral terms involving the normal component of the yield stress tensor 𝜏𝑌 (0) 𝑥𝑥 in (A.2) are equal at leading order: and therefore cancel each other. Using also definition (4.12) for the enstrophy 𝜑, equation (A.2) can be rewritten as After calculations, the integral term in (A.4) at leading order can be expressed as For a Newtonian fluid (𝜏 𝑐 = 0), the system of equations (4.40)-(4.42) reduces to

	while at order 1 we obtain
				ℎ
	𝐵𝑖 𝑅𝑒	∫	( 𝜏 𝑌 (0) 𝑥𝑧 +	1 𝐵𝑖	𝜕𝑢 (0) 𝜕𝑧	)	𝜕𝑢 (1) 𝜕𝑧	𝑑𝑧 =	𝜆ℎ𝑈 (1) 𝑅𝑒	.	(A.6)
			0
	Hence, we can write
				ℎ													ℎ
	𝐵𝑖 𝜀𝑅𝑒	∫	( 𝜏 𝑌 𝑥𝑧 +	1 𝐵𝑖	𝜕𝑢 𝜕𝑧	) 𝜕𝑢 𝜕𝑧	𝑑𝑧 =	𝜆ℎ𝑈 (0) 𝜀𝑅𝑒	+ 2	𝜆ℎ𝑈 (1) 𝑅𝑒	-	𝐵𝑖 𝑅𝑒	∫	𝜏 𝑌 (0) 𝑥𝑧	𝜕𝑢 (1) 𝜕𝑧	𝑑𝑧 + 𝑂(𝜀),	(A.7)
				0														0
	which leads to the averaged energy equation expressed in (4.29).
	Appendix B. Newtonian case
	𝜕 𝜕𝑡 + 𝜕ℎ 𝜕𝑡 𝜕 𝜕𝑡	( 𝜕 𝜕𝑧 𝑢 2 2 [ 𝑤 ) + 𝜕ℎ𝑈 + ( 𝜕𝑥 (ℎ𝑈 ) + 𝜕𝑥 𝜕 𝜕𝑥 𝑢 2 2 -[ 𝑢 𝑥 tan 𝜃 ( 𝑢 2 2 -𝜀𝐹 𝑟 2 + 𝑥 tan 𝜃 𝜀𝐹 𝑟 2 𝑧 𝐹 𝑟 2 ) = 0 𝜕 ( ℎ𝑈 2 + Π ) = [ 𝑔ℎ sin 𝜃 -3𝜈 + 𝑧 𝐹 𝑟 2 ) + 𝑝𝑤 𝐹 𝑟 2 -+ 𝐵𝑖 𝑝𝑢 𝐹 𝑟 2 𝜀𝑅𝑒 𝜏 𝑥𝑧 𝑢 + -𝐵𝑖 𝑅𝑒 𝑅𝑒 𝜏 𝑥𝑥 𝑢 𝐵𝑖 𝜏 𝑥𝑥 𝑤 ] ] = -2𝐵𝑖 𝑅𝑒 𝑈 ℎ ] + 7 6 𝑔 h sin 𝜃 𝜑 [ 𝜑 -𝑈 2 𝜏 𝑌 𝑥𝑥 ] 5ℎ 2	𝜕𝑢 𝜕𝑥	-	𝐵𝑖 𝜀𝑅𝑒	( 𝜏 𝑌 𝑥𝑧 +	1 𝐵𝑖	𝜕𝑢 𝜕𝑧	) 𝜕𝑢 𝜕𝑧	(B.1) (A.1) (B.2)
	Introducing the leading-order representation for the pressure (3.10), (3.16) and the normal stress (3.9) and averaging
	equation (A.1) over the depth leads to: 𝜕 𝜕𝑡 ( ℎ⟨𝑢 2 ⟩ 2 + ℎ 2 2𝐹 𝑟 2 ) + 𝜕 𝜕𝑥 ⎛ ⎜ ⎜ ℎ⟨𝑢 3 2 𝜕 𝜕𝑡 (ℎ𝑒) + 𝜕 𝜕𝑥 [ ℎ𝑈 ( 𝑒 + Π ℎ )] = 𝑈	⟩ [ 𝑔ℎ sin 𝜃 -3𝜈 + ℎ 2 𝑈 𝐹 𝑟 2 -2 𝑅𝑒 𝐵𝑖 𝑈 h ] ∫ + ℎ	𝜏 𝑌 (0) 𝑥𝑥 𝑢𝑑𝑧 ⎟ ⎟ ⎞ 7 18 𝑈 𝑔ℎ sin 𝜃 𝜑	[ 𝜑	𝑈 2 5ℎ 2	] .	(B.3)
																		⎜ ⎝	ℎ-ℎ 𝑝	⎟ ⎠
																		ℎ	ℎ
																		=	𝜆ℎ𝑈 𝜀𝑅𝑒	-	2𝐵𝑖 𝑅𝑒	∫	𝜏 𝑌 (0) 𝑥𝑥	𝜕𝑢 𝜕𝑥	𝑑𝑧 -	𝐵𝑖 𝜀𝑅𝑒	∫	( 𝜏 𝑌 𝑥𝑧 +	1 𝐵𝑖	𝜕𝑢 𝜕𝑧	) 𝜕𝑢 𝜕𝑧	𝑑𝑧 (A.2)
																		ℎ-ℎ 𝑝	0
				ℎ										ℎ
	𝜕 𝜕𝑥		∫	𝜏 𝑌 (0) 𝑥𝑥 𝑢𝑑𝑧 =		∫	𝜏 𝑌 (0) 𝑥𝑥	𝜕𝑢 𝜕𝑥	𝑑𝑧 + 𝑂(𝜀) = 𝛿	𝜆𝜋 4	ℎ 𝑝 (ℎ -ℎ 𝑝 )	𝜕ℎ 𝜕𝑥	,	(A.3)
	ℎ-ℎ 𝑝							ℎ-ℎ 𝑝
	𝜕 𝜕𝑡	(	ℎ𝑈 2 2	+	ℎ 3 𝜑 2		+	ℎ 2 2𝐹 𝑟 2	)	+	𝜕 𝜕𝑥	(	ℎ⟨𝑢 3 2	⟩	+	ℎ 2 𝑈 𝐹 𝑟 2	)	=
																		ℎ
																		=	𝜆ℎ𝑈 𝜀𝑅𝑒	-	𝐵𝑖 𝜀𝑅𝑒	∫	( 𝜏 𝑌 𝑥𝑧 +	1 𝐵𝑖	𝜕𝑢 𝜕𝑧	) 𝜕𝑢 𝜕𝑧	𝑑𝑧 -	𝜕 𝜕𝑥	(	ℎ⟨𝑢 ′3 2	⟩	)	(A.4)
																		0
				ℎ
	𝐵𝑖 𝜀𝑅𝑒	∫	( 𝜏 𝑌 (0) 𝑥𝑧 +	1 𝐵𝑖	𝜕𝑢 (0) 𝜕𝑧	)	𝜕𝑢 (0) 𝜕𝑧	𝑑𝑧 =	𝜆ℎ𝑈 (0) 𝜀𝑅𝑒	,	(A.5)
				0

  The coefficients 𝑅 𝑚 𝑗 (𝑗 = ℎ, 𝑈 , 𝜑) come from the right-hand side of the momentum equation(4.22), and are given by: Lastly, the coefficients 𝑅 𝑒 𝑗 (𝑗 = ℎ, 𝑈 , 𝜑) come from the right-hand side of the energy equation (4.35), and are given by: 𝛽2 (ℎ 𝑝 )𝑟(ℎ 𝑝 ) (C.11) 𝛽 2 (ℎ 𝑝 )𝑟(ℎ 𝑝 ) (C.12)The dispersion relation (6.4) is obtained by equating the determinant of the system (C.1) to zero to have a non-trivial solution:In the long-wave limit we can write: 𝑐 = 𝑐 0 + 𝜀𝑘𝑐 1 + 𝑂(𝜀 2 ). Substituting this expansion into relation (C.14), we obtain after calculations:

	𝐿 𝑚 𝑢 =	i𝜀𝑘 2	[	3	(	𝑈 (0) ) 2 + 3𝜑 (0) +	2 𝐹 𝑟 2	-2𝑈 (0) 𝑐	]	,	(C.6)
	𝐿 𝑒 𝜑 =	i𝜀𝑘 2		(	3𝑈 (0) -𝑐	)	.	(C.7)
	𝑅 𝑚 ℎ =	2𝜆 𝑅𝑒		(	1 +	1 4	ℎ 𝑝 +	1 4	ℎ 2 𝑝	) ( 1 +	1 2	ℎ 𝑝	) -1 [ 𝛼 1 (ℎ 𝑝 ) +	7 360	𝜑 (0) 𝜆 2	𝛽 1 (ℎ 𝑝 ) ]
																		+	7 135	𝜆 3 𝜑 (0) 𝑅𝑒	( 1 +	3 8	ℎ 𝑝 -	39 16	ℎ 2 𝑝 -	5 8	ℎ 3 𝑝	) ( 1 -2ℎ 𝑝 + ℎ 2 𝑝	) 2 ( 1 +	1 2	ℎ 𝑝	) -1	𝛽 2 (ℎ 𝑝 ), (C.8)
	𝑅 𝑚 𝑈 = -	3 𝑅𝑒			( 1 -ℎ 𝑝	) -1	( 1 +	1 2	ℎ 𝑝	) -1 [ 𝛼 1 (ℎ 𝑝 ) +	7 360	𝜑 (0) 𝜆 2	𝛽 1 (ℎ 𝑝 )	]
																		-	7 45	𝜆 2 𝜑 (0) 𝑅𝑒	(	1 -ℎ 𝑝	) 3	( 1 +	5 4	ℎ 𝑝	) ( 1 +	1 2	ℎ 𝑝	) -1	𝛽 2 (ℎ 𝑝 ), (C.9)
	𝑅 𝑚 𝜑 =	7 6	𝜆 𝑅𝑒			𝛽 2 (ℎ 𝑝 ) 𝜑 (0)	,	(C.10)
	𝑅 𝑒 ℎ =	2𝜆𝑈 (0) 𝑅𝑒	(	1 +	1 4	ℎ 𝑝 +	1 4	ℎ 2 𝑝	) ( 1 +	1 2	ℎ 𝑝	) -1 [ 𝛼 1 (ℎ 𝑝 ) +	7 1080	𝜑 (0) 𝜆 2	𝛽 1 (ℎ 𝑝 )𝑟(ℎ 𝑝 ) ]
															+	7 405	𝜆 3 𝑅𝑒	𝑈 (0) 𝜑 (0)	(	1 +	3 8	ℎ 𝑝 -	39 16	ℎ 2 𝑝 -	5 8	ℎ 3 𝑝	) ( 1 -2ℎ 𝑝 + ℎ 2 𝑝	) 2 ( 1 +	1 2	ℎ 𝑝	) -1
	𝑅 𝑚 𝑈 = -	3𝑈 (0) 𝑅𝑒	( 1 -ℎ 𝑝	) -1	( 1 +	1 2	ℎ 𝑝	) -1 [ 𝛼 1 (ℎ 𝑝 ) +	7 1080	𝜑 (0) 𝜆 2	𝛽 1 (ℎ 𝑝 )𝑟(ℎ 𝑝 )	]
																		-	𝜆 2 𝑅𝑒	7 135	𝑈 (0) 𝜑 (0)	(	1 -ℎ 𝑝	) 3	( 1 +	5 4	ℎ 𝑝	) ( 1 +	1 2	ℎ 𝑝	) -1
	𝑅 𝑚 𝜑 =	7 18	𝜆 𝑅𝑒	𝑈 (0) 𝜑 (0)	𝛽 2 (ℎ 𝑝 )𝑟(ℎ 𝑝 )	(C.13)
	[	i𝑅 𝑒 𝜑 + 𝜀𝑘	(	𝑐 -3𝑈 (0) )	∕2 ]	[ i𝑅 𝑚 ℎ + i	( 𝑐 -𝑈 (0) )	𝑅 𝑚 𝑈 + 𝜀𝑘	(	𝑐 2 -1∕𝐹 𝑟 2 )	-𝜀𝑘	(	2𝑐 -𝑈 (0) )	𝑈 (0) -3𝜀𝑘𝜑 0	]
			-	[	i𝑅 𝑚 𝜑 -𝜀𝑘	] [ i𝑅 𝑒 ℎ + i	(	𝑐 -𝑈 (0) )	𝑅 𝑒 𝑈 -𝜀𝑘	( 2𝑐 -𝑈 (0) ) ( 𝑈 (0) ) 2 + 𝜀𝑘	( 𝑐 2 -1∕𝐹 𝑟 2 -3𝜑 (0) )	𝑈 (0)	]	= 0
															0) +	𝐹 𝑟 2 1	)	,	(C.2) (C.14)
	𝐿 𝑚 𝑢 = i𝜀𝑘		( 𝑈 (0) -𝑐	)	,	(C.3)
	𝐿 𝑚 𝜑 = i𝜀𝑘.												(C.4)
	The coefficients 𝐿 𝑒 𝑐 1 = i𝑅𝑒 3 [ 2𝜆 2 5	(1 -ℎ 𝑝 ) 2	(	1 + ℎ 𝑝 + ℎ 2 𝑝 -	ℎ 3 𝑝 4	-	ℎ 4 𝑝 4	)	-	1 𝐹 𝑟 2
	𝐿 𝑒 ℎ =	i𝜀𝑘 2		[	( 𝑈 (0) ) 2 (	𝑈 (0) -𝑐	)	+ 3𝜑 (0) ( 3𝑈 (0) -𝑐	)	+	2 𝐹 𝑟 2	( 2𝑈 (0) -𝑐	) ]	,	(C.5)

𝑗 (𝑗 = ℎ, 𝑈 , 𝜑) come from the left-hand side of the energy equation (4.35), and are given by:

𝑐 0 = 𝜆(1 -ℎ 𝑝 ),

(C.15)

] , (C.

[START_REF] Balmforth | Viscoplastic flow over an inclined surface[END_REF]

which corresponds to (6.4).

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.