G Chambon

Perrine Freydier

M Naaim

J.-P Vila

Asymptotic expansion of velocity field within the front of viscoplastic surges: comparison with experiments

This paper investigates the internal dynamics (velocity profiles and shear rate) of freesurface surges made of viscoplastic fluids. Compared to fluids without a yield stress, additional complexity arises from the possible coexistence of sheared and unsheared (or pseudo-plug) zones in the flow. Expanding on the thin-layer approach of Fernandez-Nieto et al. (2010), we derive formal asymptotic expansions of the velocity field and discharge up to O(ǫ) with respect to flow aspect ratio ǫ. Detailed comparisons between these theoretical predictions and experimental data reveal that, while the leading-order approximation (equivalent to a lubrication model) satisfactorily accounts for the global dynamics of the flow, considering O(ǫ) correction terms is required to capture the evolution of velocity and shear rate close to the tip. Notably, these correction terms are responsible for the vanishing of the unsheared layer in the tip region, a feature clearly observed in the experiments. Differences between the leading-order and O(ǫ) models appear to be enhanced by the viscoplastic character of the fluid. In particular, O(ǫ) correction terms related to the existence of O(1) plastic normal stresses in the pseudo-plug layer, play a critical role. This study provides important insights for future development of consistent shallow-water models adapted to viscoplastic materials.

Introduction

Since the pioneering work of Saint-Venant (1871), the modelling of free-surface flows via thin-layer approaches has been the subject of countless studies in hydraulics and fluid mechanics (e.g., [START_REF] Whitham | Linear and Nonlinear Waves[END_REF][START_REF] Chanson | The Hydraulics of Open Channel Flow: An Introduction[END_REF][START_REF] Kalliadasis | Falling Liquid Films[END_REF]. In these models, the characteristic flow thickness H is assumed to be small compared to the characteristic flow length L, i.e. flow aspect ratio ǫ = H/L ≪ 1. This includes lubrication models, in which inertia terms are neglected, as well as more general shallow-water (or Saint-Venantlike) models. Lubrication models typically reduce to a single advection-diffusion equation for the evolution of flow height h [START_REF] Benney | Long waves on liquid films[END_REF][START_REF] Hunt | Newtonian fluid mechanics treatment of debris flows and avalanches[END_REF][START_REF] Huppert | Gravity currents: a personal perspective[END_REF]; Ancey et al. 2009). Shallow-water models are generally formulated as systems of two coupled equations (in 2D) for flow height h and discharge q or depth-averaged velocity ū = q/h [START_REF] Savage | The motion of a finite mass of granular material down a rough incline[END_REF][START_REF] Ancey | Existence and features of similarity solutions for non-Boussinesq gravity currents[END_REF][START_REF] Iverson | Mechanics of debris flows and rock avalanches[END_REF][START_REF] Noble | Thin power-law film flow down an inclined plane: Consistent shallow-water models and stability under large-scale perturbations[END_REF], although recently more robust three-equation models, including an additional enstrophy variable, have also been proposed [START_REF] Richard | A three-equation model for thin films down an inclined plane[END_REF]. In these models, fluid rheology is accounted for through closure terms that implicitly depend on the distribution of longitudinal velocity u(y) (with y the cross-stream coordinate) across the flowing layer [START_REF] Ng | Roll waves on a shallow layer of mud modelled as a power-law fluid[END_REF][START_REF] Ruyer-Quil | Improved modeling of flows down inclined planes[END_REF][START_REF] Hogg | The effects of hydraulic resistance on dam-break and other shallow inertial flows[END_REF]. Most of the above-mentioned references considered the case of Newtonian or generalised Newtonian (e.g., power-law) fluids, either in laminar or turbulent regimes. In this work, we focus on a more complex class of fluids, namely viscoplastic fluids characterised by a yield stress [START_REF] Balmforth | Yielding to stress: Recent developments in viscoplastic fluid mechanics[END_REF][START_REF] Coussot | Yield stress fluid flows: A review of experimental data[END_REF]. These materials are involved in numerous industrial [START_REF] Roussel | Fyfty-cent rheometer" for yield stress measurements : From slump to spreading flow[END_REF][START_REF] Frigaard | Bingham's model in the oil and gas industry[END_REF]) and geophysical applications [START_REF] Laigle | Numerical modeling of mudflows[END_REF][START_REF] Ancey | Plasticity and geophysical flows: A review[END_REF], for which accurate thin-layer models would be needed. Yet, as will be discussed, the viscoplastic rheology poses specific modelling issues related to the possible coexistence of unyielded (solid-like) and yielded (fluid-like) zones within the flows. Hence, attempts at developing lubrication [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]Ancey & Cochard 2009;[START_REF] Bernabeu | Numerical modeling of non-Newtonian viscoplastic flows: Part ii. Viscoplastic fluids and general tridimensional topographies[END_REF][START_REF] Liu | Two-dimensional viscoplastic dambreaks[END_REF] or shallow-water [START_REF] Coussot | Steady, laminar, flow of concentrated mud suspensions in open channel[END_REF][START_REF] Liu | Roll waves on a layer of a muddy fluid flowing down a gentle slope -A Bingham model[END_REF][START_REF] Piau | Flow of a yield stress fluid in a long domain. Application to flow on an inclined plane[END_REF][START_REF] Balmforth | Roll waves in mud[END_REF]) models for viscoplastic materials remained relatively rare.

Formally, thin-layer models can be derived from the primitive Cauchy equations through asymptotic expansions with respect to flow aspect ratio ǫ. Lubrication models proceed from a leading-order expansion [START_REF] Benney | Long waves on liquid films[END_REF][START_REF] Ruyer-Quil | Modeling film flows down inclined planes[END_REF]. At the next level of approximation, derivation of consistent shallow-water models requires expansion of the velocity profile u(y) up to at least O(ǫ) terms [START_REF] Ruyer-Quil | Improved modeling of flows down inclined planes[END_REF][START_REF] Noble | Thin power-law film flow down an inclined plane: Consistent shallow-water models and stability under large-scale perturbations[END_REF]). Yet, most shallow-water models used in practical applications, either for Newtonian of non-Newtonian fluids, are based on simpler, ad-hoc assumptions. Typically, closure terms are postulated empirically, or the velocity profile is assumed identical to that obtained in a steady uniform regime (see, e.g., [START_REF] Liu | Roll waves on a layer of a muddy fluid flowing down a gentle slope -A Bingham model[END_REF][START_REF] Coussot | Steady, laminar, flow of concentrated mud suspensions in open channel[END_REF], for viscoplastic fluids). Such ad-hoc shallow-water models generally fail to predict the correct instability threshold for the flows (Ruyer-Quil & Manneville 2000; [START_REF] Kalliadasis | Falling Liquid Films[END_REF]. From an operational point of view, these models nevertheless proved capable of producing relatively accurate predictions in terms of, e.g., front propagation velocity or free-surface shape [START_REF] Hogg | The effects of hydraulic resistance on dam-break and other shallow inertial flows[END_REF][START_REF] Ancey | Viscoplastic dambreak waves: Review of simple computational approaches and comparison with experiments[END_REF][START_REF] Saingier | On the front shape of an inertial granular flow down a rough incline[END_REF][START_REF] Kirstetter | Modeling rain-driven overland flow: Empirical versus analytical friction terms in the shallow water approximation[END_REF].

Expansions of the velocity field up to O(ǫ), and the use thereof in the formulation of consistent shallow-water models, have been worked out for Newtonian and powerlaw fluids. The handling of the two types of constitutive laws is largely similar, despite a singularity in the apparent viscosity at the free surface in the power-law case. [START_REF] Ruyer-Quil | Improved modeling of flows down inclined planes[END_REF] and [START_REF] Ruyer-Quil | Wavy regime of a power-law film flow[END_REF] used a weighted-residual method to express deviations from the steady uniform profile in terms of polynomial test functions (see also [START_REF] Amaouche | Film flow for power-law fluids: modeling and linear stability[END_REF]. Alternatively, [START_REF] Fernandez-Nieto | Shallow Water equations for non-Newtonian fluids[END_REF] and [START_REF] Noble | Thin power-law film flow down an inclined plane: Consistent shallow-water models and stability under large-scale perturbations[END_REF] performed a formal asymptotic expansion of the primitive equations, in a regime where these equations reduce to steady uniform flow at leading-order. At O(ǫ), the velocity profile u(y) then explicitly depends on flow height gradient ∂ x h (with x the streamwise coordinate), and expresses as a sum of the steady uniform profile and of inertia correction terms. Shallow-water models derived from either of these two approaches were shown to correctly capture the flow instability threshold. More generally, [START_REF] Boutounet | Multi-regime shallow free surface laminar flow models for quasi-Newtonian fluids[END_REF] explained that a whole family of consistent shallow-water models, which are all equivalent at O(ǫ), can be constructed. For practical applications, additional model selection criteria would thus be needed.

Far fewer studies focused on expansions of the velocity field in viscoplastic fluids described by Bingham or Herschel-Bulkley constitutive laws. In a steady uniform regime, the longitudinal velocity profile u(y) is composed of a yielded layer at the base, overlain by an unyielded, plug layer close to the free-surface (e.g., [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]. [START_REF] Balmforth | A consistent thin-layer theory for Bingham plastics[END_REF] showed that, in thin non-uniform flows, this plug layer becomes a pseudo-plug in which the strain rate is O(ǫ). These authors worked out an asymptotic expansion of the velocity profile in the pseudo-plug up to O(ǫ) in the inertialess limit. Recently, [START_REF] Liu | Viscoplastic surges down an incline[END_REF] investigated the transition between a true plug and a pseudo-plug. [START_REF] Balmforth | Roll waves in mud[END_REF] used the weighted-residual approach to derive an O(ǫ) velocity approximation in the sheared layer, but did not consider the correction terms in the pseudo-plug. [START_REF] Fernandez-Nieto | Shallow Water equations for non-Newtonian fluids[END_REF] applied their formal asymptotic approach to the case of a Bingham fluid and showed that, besides the inertia correction terms already present for fluids without a yield stress, the O(ǫ) velocity approximation also involves correction terms due to the existence of O(1) normal stresses in the pseudoplug. These latter authors also proposed the -to date -only consistent shallow-water model for a viscoplastic rheology. However, variations of pseudo-plug thickness along the flow are not considered in this model. Lubrication-like models, on the contrary, predict a dependence of pseudo-plug thickness on flow height gradient ∂ x h [START_REF] Balmforth | A consistent thin-layer theory for Bingham plastics[END_REF]Ancey & Cochard 2009;[START_REF] Freydier | Experimental characterization of velocity fields within the front of viscoplastic surges down an incline[END_REF], and one would like O(ǫ) models to recover this feature when inertia terms are neglected.

In this context, the objectives of the present study are twofold. First, expanding on the approach of [START_REF] Fernandez-Nieto | Shallow Water equations for non-Newtonian fluids[END_REF], we derive a complete O(ǫ) expansion of longitudinal velocity profile u(y) in a Herschel-Bulkley thin free-surface flow, accounting for possible variations in pseudo-plug thickness. As will be shown, the expression obtained depends on gradients of both flow height and pseudo-plug thickness. The prospect of formulating consistent shallow-water models based on this asymptotic expansion of velocity, is left for future work. Instead, we assess the differences between the leadingorder and O(ǫ) velocity approximations in the specific case of travelling-wave solutions, for which the evolution of free surface can be computed without additional developments. The second objective is then to compare these asymptotic velocity predictions to highresolution experimental data recently obtained by our group [START_REF] Freydier | Experimental characterization of velocity fields within the front of viscoplastic surges down an incline[END_REF]. Our previous works showed in particular that, in viscoplastic free-surface surges, longitudinal velocity profiles measured close to the tip clearly diverge from the predictions of a leading-order lubrication model. Our goal here is therefore to investigate whether an O(ǫ) approximation provides a better match to the data.

Formal asymptotic expansion of the velocity field is derived in section 2. Section 3 examines the influence of the different correction terms involved in the O(ǫ) approximation, as a function of the relevant non-dimensional parameters. Experimental results are briefly recalled, and compared to the theoretical predictions, in section 4. Finally, sections 5 and 6 discuss and summarise the main outcomes of the study. Additional results obtained with Newtonian fluids are presented in Appendix B.

Thin-layer asymptotic expansion

Let us consider a free-surface viscoplastic surge propagating down an inclined plane forming an angle θ to the horizontal (figure 1). The flow is driven by gravity g. The x-axis (hereafter called longitudinal ) coincides with the base plane in the streamwise direction, while the y-axis (hereafter called vertical ) is orthogonal to the bed, pointing upwards. Flow is assumed to be 2D, i.e. independent of the spanwise direction. Fluid density is denoted ρ, local flow thickness h(x, t), and fluid velocity u(x, y, t).

The viscoplastic fluid is assumed to obey Herschel-Bulkley constitutive law: where τ ij = σ ij + pδ ij is the extra-stress tensor (with σ ij the total stress and p the pressure) and γij = ∂ j u i + ∂ j u i is the strain-rate tensor. The tensorial norm is defined as

   τ ij = τ c | γ| + K| γ| n-1 γij if |τ | > τ c , γij = 0 if |τ | τ c , (
|T | = (0.5 T ij T ij ) 0.
5 for any second-order tensor T . Rheological parameters τ c , K, and n correspond to the yield stress, the consistency, and the power index of the material, respectively. We assume n > 0, and shall consider more specifically the two cases n = 1 (Bingham fluid) and n ≈ 0.4 (shear-thinning value representative of most real viscoplastic fluids). In the following, for |τ | > τ c , we split the extra-stress tensor into its viscous and plastic contributions:

τ ij = τ v ij + τ p ij , with τ v ij = K| γ| n-1 γij and τ p ij = (τ c /| γ|) γij .

Primitive equations

Let us denote u = u x and v = u y the longitudinal and vertical velocity components, respectively. The primitive equations governing free-surface flows of viscoplastic fluids are derived in numerous previous references (e.g., [START_REF] Piau | Flow of a yield stress fluid in a long domain. Application to flow on an inclined plane[END_REF]Ancey & Cochard 2009;[START_REF] Fernandez-Nieto | Shallow Water equations for non-Newtonian fluids[END_REF]. Here, we recall them directly in dimensionless form. Hence let us define the Reynolds, Froude and Bingham numbers as:

Re = ρU 2-n H n K , Fr = U √ gH cos θ , Bi = τ c K H U n , (2.2)
where H, L and U represent a characteristic thickness, a characteristic length, and a characteristic longitudinal velocity of the flow, respectively. Pressure p is scaled by ρgH cos θ, and shear rate | γ| is scaled by U/H. We choose to scale the viscous contribution to stresses by K(U/H) n , and the plastic contribution by τ c . Lastly, the stress norm |τ | is scaled by K(U/H) n , such that dimensionless yield condition writes |τ | = Bi . Dimensionless mass and momentum conservation equations then read:

∂ x u + ∂ y v = 0 (2.3) ǫRe(∂ t u + u∂ x u + v∂ y u) = - ǫRe Fr 2 ∂ x p + λ + ∂ y τ v xy + Bi ∂ y τ p xy + ǫ(∂ x τ v xx + Bi ∂ x τ p xx ) (2.4) ǫ 3 Re(∂ t v + u∂ x v + v∂ y v) = - ǫRe Fr 2 (1 + ∂ y p) + ǫ 2 (∂ x τ v xy + Bi ∂ x τ p xy ) + ǫ(∂ y τ v yy + Bi ∂ y τ p yy ) (2.5)
where we recall that ǫ = H/L is flow aspect ratio, and the driving parameter λ is given by:

λ = ρgH n+1 sin θ KU n = Re
Fr 2 tan θ.

(2.6)

The dimensionless expressions of viscous and plastic stresses write, for |τ | > Bi :

τ v xx = -τ v yy = 2ǫ| γ| n-1 ∂ x u (2.7) τ v xy = | γ| n-1 (∂ y u + ǫ 2 ∂ x v), (2.8)
and

τ p xx = -τ p yy = 2ǫ | γ| ∂ x u
(2.9)

τ p xy = 1 | γ| (∂ y u + ǫ 2 ∂ x v).
(2.10)

The equality τ xx = -τ yy ensues from mass conservation (2.3). Lastly, the dimensionless stress norm expresses as: (2.11) and the dimensionless shear rate:

|τ | = (τ v xy + Bi τ p xy ) 2 + (ǫτ v xx + Bi τ p xx ) 2 ,
| γ| = (∂ y u + ǫ 2 ∂ x v) 2 + 4ǫ 2 (∂ x u) 2 .
(2.12)

The above balance equations are to be complemented with appropriate boundary conditions. At the base we consider a classical no-penetration and no-slip condition: u |y=0 = v |y=0 = 0. At the free surface, influence of surface tension is neglected and a stress-free condition is assumed: (τ -pI) |y=h • n = 0, with n = (-∂ x h, 1) and I the identity tensor. Rearranging terms, we obtain:

(τ v xy + Bi τ p xy ) |y=h = 2ǫ∂ x h 1 -ǫ 2 (∂ x h) 2 (ǫτ v xx + Bi τ p xx ) |y=h (2.13) p |y=h = - Fr 2 Re 1 + ǫ 2 (∂ x h) 2 1 -ǫ 2 (∂ x h) 2 (ǫτ v xx + Bi τ p xx ) |y=h .
(2.14)

Finally, the following kinematic condition holds at the free surface:

∂ t h + u |y=h ∂ x h = v |y=h .
(2.15)

For completeness, we also recall the depth-averaged mass-conservation equation:

∂ t h + ∂ x q = 0. (2.16)
where q = h 0 u dy is flow discharge. This equation is readily obtained from integration of local mass conservation (2.3) and boundary conditions.

Flow regime

Following the approach of [START_REF] Fernandez-Nieto | Shallow Water equations for non-Newtonian fluids[END_REF], we aim to derive approximations of the velocity field in the flow through formal asymptotic expansions of the previous primitive equations with respect to ǫ. The following assumptions are introduced: (i) ǫ ≪ 1 (thin-layer approximation). (ii) The factor ǫRe in front of inertia terms in (2.4) is considered O(ǫ). (iii) As the study concerns flows in which plastic effects are significant, Bingham number Bi is considered O(1). (iv) Finally, the factor ǫRe/Fr 2 in front of pressure gradient terms in (2.4)-(2.5) is also considered O(1). This latter assumption constitutes the main difference with the work of [START_REF] Fernandez-Nieto | Shallow Water equations for non-Newtonian fluids[END_REF], in which ǫRe/Fr 2 was considered O(ǫ). Accordingly, the flow regime investigated in the present study corresponds, at leading order, to a balance between pressure gradient, gravity and shear stress effects, as in classical lubrication models. Note that the assumption ǫRe/Fr 2 = O(1) implicitly requires the slope angle θ to be small, so that λ = ǫRe/Fr 2 (tan θ/ǫ) can also be O(1).

We assume the existence of regular expansions of the form: f = f (0) + ǫf (1) + ..., for all variables of interest, namely longitudinal and vertical velocities u and v, pressure p, and stress components τ ij . As will be shown, these expansions can effectively be formally worked out up to O(ǫ). Let us note, however, that the residues, whose complete evaluation lies outside the scope of this paper, are not necessarily O(ǫ 2 ) and may involve terms depending on non-integer powers ǫ a (with a > 1) [START_REF] Fernandez-Nieto | Shallow Water equations for non-Newtonian fluids[END_REF].

Leading-order expansion

At O(1) with respect to ǫ, integration of equations (2.4)-(2.5) leads to linear profiles for the shear stress and pressure:

τ v,(0) xy + Biτ p,(0) xy = Λ(x, t) (h -y),
(2.17)

p (0) = (h -y), (2.18)
where we introduced

Λ(x, t) = λ - ǫRe Fr 2 ∂ x h = λ (1 -ǫ cot θ ∂ x h). (2.19)
Moreover, Herschel-Bulkley law reduces to, for |τ | > Bi :

τ v,(0) xy + Bi τ p,(0) xy = sgn(∂ y u (0) )(∂ y u (0) ) n + Bi , (2.20) while τ p,(0) xx = τ v, (0) xx
= 0. From these leading-order expressions, we expect the existence of an unsheared layer, in which | γ| (0) = |∂ y u (0) | = 0, for y h -h p . The thickness of this layer is given by:

h p = Bi Λ .
(2.21)

We may assume, without loss of generality, that ∂ y u (0) 0. Integration of (2.17) then leads to the classical lubrication velocity profile: (2.22) where ζ = 1 -y/(h -h p ) is a reduced vertical coordinate, and

u (0) = u (0) pl 1 -ζ (n+1)/n if y < h -h p u (0) pl if y h -h p ,
u (0) pl = n n+1 Λ 1/n (h -h p ) (n+1)/n . (2.23)
It shall be emphasised that h p , Λ, and thus u

(0)
pl , are all functions of x and t through the term ∂ x h in (2.19). If ǫRe/Fr 2 = 0, expressions (2.22)-(2.23) reduce to the velocity profile obtained in a steady uniform regime, with a constant plug thickness h p = Bi /λ (e.g., [START_REF] Chambon | Experimental investigation of viscoplastic free-surface flows in steady uniform regime[END_REF].

The leading-order approximation of discharge q immediately proceeds from (2.22):

q (0) = n 2n+1 Λ 1/n (h -h p ) (n+1)/n h + n n+1 h p .
(2.24)

In the following (see §2.5), leading-order expressions of vertical velocity v (0) and derivatives ∂ x,t u (0) will also be required. pl appears incompatible with a strictly zero shear rate | γ| for y h -h p , as would be expected in a "true" unyielded plug. As observed by [START_REF] Balmforth | A consistent thin-layer theory for Bingham plastics[END_REF], this apparent paradox is overcome by noting that, since ∂ y u (0) = 0 in this layer, the expansion of longitudinal velocity takes the form:

u = u (0) pl (x, t) + ǫu (1)
pl (x, y, t) + .... Accordingly, the leading-order expression of shear rate | γ| becomes: (2.26) such that the expansion of Herschel-Bulkley law also needs to be revisited. In particular, unlike for y < h -h p , the plastic normal stress τ p xx involves an O(1) contribution:

| γ| = ǫ ∂ y u (1) pl 2 + 4 ∂ x u (0) pl 2 + ...,
τ p,(0) xx = 2∂ x u (0) pl ∂ y u (1) pl 2 + 4 ∂ x u (0) pl 2 , (2.27)
while the leading-order expression of shear stress becomes:

τ v,(0) xy + Bi τ p,(0) xy = Bi τ p,(0) xy = Bi ∂ y u (1) pl ∂ y u (1) pl 2 + 4 ∂ x u (0) pl 2 .
(2.28)

It ensues that the stress norm expresses as:

|τ | = Bi + O(ǫ n ), (2.29)
where the residue O(ǫ n ) is necessarily positive. Hence, plastic normal stresses actually make up for the small values of shear stress, and are sufficient for the yield criterion to be slightly exceeded. The unsheared layer predicted by the leading-order velocity expansion, should thus rather be considered as a pseudo-plug, in which the shear rate is O(ǫ). Accordingly, the interface at y = h -h p should be regarded as a fake yield surface.

Combining equations (2.28) and (2.17) leads to:

∂ y u (1) pl = 2(h -y)/h p 1 -((h -y)/h p ) 2 ∂ x u (0) pl , (2.30)
where (2.21) was used to express pseudo-plug thickness h p . From (2.26), the shear rate in the pseudo-plug then expresses as:

| γ| = ǫ 2 1 -((h -y)/h p ) 2 ∂ x u (0) pl + ... (2.31)
Finally, leading-order plastic normal stress in the pseudo-plug writes:

τ p,(0) xx = 1 - h -y h p 2 sgn(∂ x u (0) pl ).
(2.32)

Note that this expression for τ

p,(0) xx requires ∂ x u (0) pl = 0. If ∂ x u (0)
pl = 0, we are either in a true plug (| γ| = 0), or the shear rate in the pseudo-plug vanishes at O(ǫ) and a leading term of higher order needs to be considered in the expansion of | γ| [START_REF] Fernandez-Nieto | Shallow Water equations for non-Newtonian fluids[END_REF]. This latter case will not be considered in this study.

O(ǫ) expansion

At O(ǫ), velocity expansions in the pseudo-plug (y h -h p ) and in the sheared layer (y < h -h p ) need to be treated separately.

In the pseudo-plug (y

h -h p )
As seen in §2.4, the O(ǫ) velocity correction in the pseudo-plug directly derives from the leading-order approximation of shear stress. Integrating (2.30) between y = h -h p and y > h -h p , and accounting for the expression of ∂ x u (0) pl , leads to:

u (1) pl = Bi U N pl (x, y, t) + u + (x, t) (2.33)
where

U N pl (x, y, t) = 2Λ (1-n)/n (h -h p ) 1/n 1 - h -y h p 2 × ∂ x h + 1 n+1 (h + nh p ) ∂ x Λ Λ sgn ∂ x u (0) pl , (2.34)
and u + (x, t) is a matching term (see §2.5.3).

In the sheared layer (y < h -h p )

Expansion at O(ǫ) of momentum balance equations (2.4) and (2.5) leads to:

∂ y τ v,(1) xy + Biτ p,(1) xy = ǫRe Fr 2 ∂ x p (1) -Bi ∂ x τ p,(0) xx + Re(∂ t u (0) + u (0) ∂ x u (0) + v (0) ∂ y u (0) ) (2.35) ǫRe Fr 2 ∂ y p (1) = -Bi ∂ y τ p,(0) xx .
(2.36)

Accounting for boundary conditions at the free surface, integration of (2.36) provides:

ǫRe Fr 2 p (1) = -Bi τ p,(0) xx , (2.37)
for all values of y. Accordingly, p (1) = 0 for y < h -h p . Equation (2.37) also shows that the two terms (ǫRe/Fr 2 )∂ x p (1) and -Bi ∂ x τ p,(0) xx , involved in (2.35), are equal everywhere in the flow. Integration of (2.35) between y < h -h p and y = h then leads to the expression of the O(ǫ) shear-stress correction in the sheared layer:

τ v,(1) xy + Bi τ p,(1) xy = -Bi π 2 h p Λ sgn ∂ x u (0) pl ∂ x Λ + Re I(x, y, t), (2.38) with I(x, y, t) = - h y ∂ t u (0) + u (0) ∂ x u (0) + v (0) ∂ y u (0) dy.
(2.39)

The first term on the r.h.s. of (2.38) arises from the calculation of

2Bi h h-hp ∂ x τ p,(0) xx dy + τ v,(1) xy + Bi τ p,(1) xy |y=h
, where the integral of

∂ x τ p,(0) xx is restricted to the range [h -h p , h], since τ p,(0) xx = 0 for y < h -h p ,
and is evaluated from (2.32).

Expansion of Herschel-Bulkley law at O(ǫ) for y < h -h p writes:

τ v,(0) xy + Bi τ p,(0) xy + ǫ(τ v,(1) xy + Bi τ p,(1) xy ) = ∂ y u (0) + ǫ∂ y u (1) n + Bi , (2.40) since | γ| (0) + ǫ| γ| (1) = |∂ y u (0) + ǫ∂ y u (1) |.
We assumed, again, that ∂ y u (0) + ǫ∂ y u (1) 0.

A Taylor expansion then leads to:

∂ y u (1) = 1 n ∂ y u (0) 1-n τ v,(1) xy + Bi τ p, (1) xy .
(2.41)

Finally, inserting equation (2.38) into (2.41) and further integrating, provides the O(ǫ) longitudinal-velocity correction for y < h -h p :

u (1) = Bi U N (x, ζ, t) + Re U I (x, ζ, t) (2.42)
where

U N (x, ζ, t) = - π 2 Λ (1-2n)/n h p (h -h p ) 1/n (1 -ζ 1/n ) sgn ∂ x u (0) pl ∂ x Λ (2.43) U I (x, ζ, t) = 1 n Λ (1-n)/n (h -h p ) 1/n 1 ζ ζ ′(1-n)/n I(ζ ′ )dζ ′ (2.44)
and we recall that ζ = 1 -y/(h -h p ). Explicit calculations of the functions I and U I are lengthy but straightforward. The final expression of U I is given in (A 2), Appendix A. Note that in this expression, following Fernandez-Nieto et al. ( 2010), all factors proportional to ∂ t h have been replaced through the use of differential identity (A 1), which derives from mass conservation equation.

Matching and final expressions

The term u + in (2.33) is determined by patching the solutions obtained in the sheared layer and the pseudo-plug, i.e. enforcing continuity of u (1) at the fake yield surface y = h -h p (or ζ = 0):

u + = Bi U N (x, 0, t) + Re U I (x, 0, t).
(2.45)

Finally, the O(ǫ) longitudinal-velocity correction can thus be written as follows:

u (1) = Bi U N (x, ζ, t) + Re U I (x, ζ, t) if y < h -h p Bi [U N pl (x, y, t) + U N (x, 0, t)] + Re U I (x, 0, t) if y h -h p , (2.46)
where the functions U N pl , U N , and U I are given by (2.34), (2.43), and (A 2), respectively. Terms in factor of Bi in (2.46) represent O(ǫ) corrections arising from the existence of O(1) plastic normal stresses in the pseudo-plug. These correction terms, denoted henceforth u

(1)

N , vanish in the absence of a yield stress, i.e. for Bi = 0. It can also be noted that the normal-stress correction term affecting the sheared layer, namely Bi U N , is proportional to ∂ x Λ and thus to gradient of pseudo-plug thickness ∂ x h p . The term Re U I (ζ) in (2.46), denoted henceforth u

(1) I , represents the O(ǫ) correction associated to inertia effects. This term is present even for fluids without a yield stress, although its expression is considerably simplified when h p = 0 (see Appendix B for the Newtonian case).

From (2.46), the O(ǫ) correction on flow discharge q is readily obtained:

q (1) = Bi Q N (x, t) + Re Q I (x, t), (2.47)
where the normal-stress correction term Q N writes: (2.48) and expression of the inertia correction term Q I is given in (A 5), Appendix A. Interestingly, it is observed that the factors proportional to ∂ x Λ in U N and U N pl cancel out in the calculation of Q N , such that this latter function only involve a factor proportional to ∂ x h.

Q N (x, t) = π 2 Λ (1-n)/n h p (h -h p ) 1/n sgn ∂ x u (0) pl ∂ x h,
Finally, let us recall that the above expansions at O(ǫ) implicitly require that ∂ x u (0) pl = 0 everywhere in the flow. Accordingly, the sign factor sgn(∂ x u (0) pl ) involved in the normalstress correction terms U N pl , U N and Q N , was assumed to remain constant along the flow.

Influence of O(ǫ) correction terms for travelling-wave solutions

This section aims to compare the leading-order and O(ǫ) asymptotic approximations established above, and to investigate the relative influences of the O(ǫ) normal-stress and inertia corrections terms. We explore configurations for which Bi and Re numbers vary between 0.5 -1.5 and 0.1 -1.5, respectively, consistently with the ranges achieved in the experiments (see §4.2). A slope angle θ = 15 • is considered, and results corresponding to power indices n = 1 and n = 0.4 are compared.

Expressions of u (0) and u (1) depend on flow thickness h, on its derivative ∂ x h, and, for u (1) , on the second derivatives ∂ xx h and ∂ xt h (through ∂ x Λ and ∂ t Λ). Quantitative evaluation of the asymptotic velocity profiles thus requires specifying a free-surface evolution. In practice, values of h can be obtained either from experimental measurements or from a theoretical solution. We take here the latter approach, focusing on specific travelling-wave solutions of the form:

h(x, t) = h tw (x f ), (3.1)
with x f = u b t -x, and u b the wave propagation speed. For such solutions, mass conservation equation (2.16) implies that

d dx f (u b h -q) = 0.
If, furthermore, u b h -q = 0 is assumed at surge tip, then the relation

q = u b h (3.2)
holds everywhere along the surge. These particular solutions are thus characterised by a constant depth-averaged velocity ū = q/h = u b . Combining (3.2) with expressions (2.24) and (2.47) for the discharge q, closed-form ordinary differential equations (ODEs) for the free-surface shape h tw (x f ) can then be obtained (see below).

As will be shown, far from the tip, these travelling-wave solutions asymptotically converge to a steady uniform flow. To simplify expressions, we shall henceforth consider the following choice of characteristic scales: (i) Characteristic velocity U is taken equal to the travelling-wave propagation speed u b . (ii) Both characteristic thickness H and characteristic length L are taken equal to the flow thickness H N achieved in steady uniform regime for the average velocity u b . Note that choosing identical vertical and longitudinal length scales formally amounts to considering ǫ = 1. Obviously, we remain however in the frame of a thin-layer assumption. With this particular choice of characteristic scales, and since u b and H N are related through the steady uniform condition, numbers λ and Bi are linked as follows:

λ = 2n + 1 n n 1 - Bi λ -(n+1) 1 + n n + 1 Bi λ -n .
(3.3) 3.1. Free-surface shape

In the following, two different free-surface models are considered. For computing leading-order velocity approximations, we define model M 1 based on the leading-order approximation of discharge q:

M 1 : q (0) h tw , dh tw dx f = h tw , (3.4)
where the expression of q (0) as a function of h tw and dh tw /dx f is provided by (2.24). Equation (3.4) constitutes a first-order ODE for the free-surface shape h tw (x f ). The boundary condition h tw (x f = 0) = 0 is enforced, such that x f represents the distance from surge tip. In practice, (3.4) is integrated numerically using python function odeint (see also [START_REF] Freydier | Experimental characterization of velocity fields within the front of viscoplastic surges down an incline[END_REF]. Due to the singularity of dh tw /dx f at surge tip (see below), numerical integration is started from a small value x f,0 , with the corresponding value h tw (x f,0 ) being obtained through a finite-difference approximation.

For computing velocity expansions at order O(ǫ), we consider model M 1+ǫ based on the O(ǫ) approximation of discharge:

M 1+ǫ : q (0) h tw , dh tw dx f + q (1) h tw , dh tw dx f , d 2 h tw dx f 2 = h tw , (3.5)
where the expression of q (1) is provided by (2.47). Equation (3.5) constitutes a secondorder ODE for h tw (x f ). It is solved by finite-difference discretisation on a regular mesh, using python function newton krylov. The same boundary condition as for model M 1 is considered at x f = 0, and we impose h tw ≈ 1 sufficiently far from the tip (in practice, the value provided by model M 1 is taken). The iterations of the relaxation algorithm are initiated with the solution of model M 1 .

Figure 2 shows comparisons between the free surfaces obtained with models M 1 and M 1+ǫ for different values of numbers Bi and Re. In all cases, the front is steep and perpendicular to channel bottom at the tip, namely |∂ x h| → ∞ as x f → 0. Values of the derivative |∂ x h| exceed 0.1 for x f < 1.5 -2, typically, a position below which the validity of the thin-layer assumption may thus become questionable (see §5.2). Note however that values of |∂ x h| in excess of 1 are reached only very close to the tip, for x f < 0.1 typically. Away from the tip, the free surface monotonically increases towards the asymptotic value of 1 corresponding to a steady uniform regime. Flow height can be considered as essentially uniform for x f > 10, typically.

Globally, it appears that the free surfaces predicted by models M 1 and M 1+ǫ remain very close, regardless of the parameter values considered (figure 2). Note that model M 1 only depends on Bingham number Bi , while model M 1+ǫ depends on both Bi and Re. For n = 0.4 (figure 2a), and when Bi is small, the free surfaces obtained with the two models are hardly distinguishable. When Bi increases, small deviations between the different solutions become visible. It is observed that the free surfaces computed with model M 1+ǫ systematically lie below those obtained with model M 1 . Moreover, for a given value of Bi , the differences between the M 1 and the M 1+ǫ solutions are larger when Re is small. Similar observations are made for n = 1 (figure 2b).

Let us note that, far from the tip (x f → ∞), expression (2.23) for the leading-order velocity in the pseudo-plug can be approximated as:

u (0) pl ≈ n n+1 λ 1/n h - Bi λ (n+1)/n . Accordingly, sgn(∂ x u (0) pl ) = sgn(∂ x h) = -1 for x f → ∞. Since, furthermore, the factor sgn(∂ x u (0)
pl ) is assumed to be constant (see §2.5.3), it is taken as uniformly equal to -1 along the flow for computing solutions of model M 1+ǫ .

Velocity profiles and free-surface velocity

Approximations of the velocity profile u(y) at leading-order and at O(ǫ) are computed from expressions (2.22) and (2.46). When limiting to leading order, free surfaces obtained from model M 1 are considered. Conversely, for computing velocity approximations at O(ǫ), free surfaces coming from model M 1+ǫ are taken into account. Considering these two, slightly-different, free-surface models ensures that all computed velocity profiles effectively comply with the condition ū = 1. In the following, for convenience, the terminology M 1 and M 1+ǫ is also used for referring to the velocity profiles computed with either one of the two models. Again, the assumption sgn(∂ x u x f = 2, slight differences between the predictions of the two models begin to appear for large values of Bi (figure 3a). These differences between the models are then magnified when moving closer to the tip, as shown in figure 3b for x f = 0.4. Also clearly visible are the velocity gradients present in the pseudo-plug with model M 1+ǫ . Furthermore, the corresponding velocity profiles present a nonphysical kink at y = h -h p . As discussed below (see §5.1), this singularity arises from the patching of the O(ǫ) correction terms at the fake yield surface. Globally, the interplay of numbers Bi and Re on the velocity profiles is rather complex. Depending on the values of these parameters, and on the vertical position y, velocities predicted by model M 1+ǫ can either be larger or smaller than those predicted by model M 1 . Figure 3c shows the typical relative magnitudes of the O(ǫ) normal-stress and inertia correction terms u N is positive and, as expected, increases with Bingham number Bi . The inertia term u

(1) I is generally negative, and its absolute value increases with Reynolds number Re. More surprisingly, the inertia correction also appears to significantly increase (in absolute value) with Bi . The counterbalancing effects between the two correction terms explain that, for a given value of Bi , the net velocity correction tends to decrease when Re increases (clearly visible in figure 3b). The net correction, however, systematically increases with Bi (as seen in figures 3a and 3b).

Figure 4 presents the evolution of free-surface velocity u s along the surge. For model

M 1 , u s = u (0) pl , while for model M 1+ǫ , u s = u (0) pl + u (1)
pl | y=h . Computed predictions have been truncated at x f = 0.1, i.e. at the position where |∂ x h| typically reaches a value of 1 (see figure 2). It was checked that, with model M 1+ǫ , the asymptotic expansion

s < u (0)
s , for x f > 0.1 in all but one of the displayed cases. With n = 0.4, Bi = 1.5 and Re = 0.1, the expansion ceases to be well-ordered at x f = 0.25, typically (red star in figure 4a). Globally, it is observed that, for the two models and all considered parameter values, surface velocity u s systematically increases towards the tip. Note that the increase of u s is not incompatible with the requirement that sgn(∂ x u (0) pl ) = -1 in model M 1+ǫ (see §3.1). Again, the counterbalancing effects of the O(ǫ) normal-stress and inertia correction terms are clearly apparent for n = 0.4 (figure 4a): The normal-stress correction tends to increase the value of u s , while the inertia term generally acts in the opposite direction. The net correction, however, remains systematically positive: surface velocities obtained with model M 1+ǫ are larger than those obtained with model M 1 , with differences becoming visible for x f < 2 -6, typically. For similar values of Bi and Re, the influence of the O(ǫ) corrections, and notably of the inertia term, appear to be less marked in the case n = 1 (figure 4b).

Pseudo-plug shear rate and effective unsheared layer

At O(ǫ), a strictly positive shear rate | γ| is predicted in the pseudo-plug (see §2.4). Integrating (2.31), the following expression is obtained for the O(ǫ) vertically-averaged shear rate | γ| pl obtained from model M 1+ǫ significantly depart from 0 for x f < 2 -6, depending on the considered values of Bi and Re. Again, note that the plots are truncated at x f = 0.1. The increase in | γ|

(1) pl towards the tip is primarily driven by the normal-stress correction, and is more pronounced for larger Bi . Conversely, for n = 0.4, the increase of | γ|

(1) pl gets generally confined closer to the tip when Re increases (counterbalancing effect, see figure 5a). For n = 1, the net O(ǫ) correction is less pronounced, and the effect of Re changes sign in the close vicinity of the tip (figure 5b). Figure 6 displays a typical example of iso-contours of the O(ǫ) shear rate | γ| (1) in the pseudo-plug. Recall that the pseudo-plug corresponds to an unsheared layer for the leading-order approximation. Similarly, one can define an effective unsheared layer at order O(ǫ), as the zone located above the particular | γ| (1) -contour corresponding to a chosen shear-rate threshold. For sufficiently large values of x f , the boundary of this effective unsheared layer coincides with the fake yield surface y = h -h p , regardless of the chosen threshold. When approaching the tip, however, the thickness of the effective unsheared layer decreases much faster than h p , and abruptly collapses at a certain value of x f . The precise position at which this collapse occurs, moves closer to the tip as the considered shear-rate threshold increases (figure 6).

Comparison with experimental results

In this section, the leading-order and O(ǫ) velocity approximations are cross-compared to experimental results. The experimental dataset consists of high-resolution measurements of the internal velocity field within laboratory viscoplastic surges [START_REF] Freydier | Experimental characterization of velocity fields within the front of viscoplastic surges down an incline[END_REF]. In parallel, free-surface shapes were also recorded. Experimental setup and methods were detailed in previous papers [START_REF] Chambon | Experimental investigation of viscoplastic free-surface flows in steady uniform regime[END_REF][START_REF] Freydier | Experimental characterization of velocity fields within the front of viscoplastic surges down an incline[END_REF]. Only the elements essential for the present comparison are recalled below.

Review of experimental methods

The experimental setup consists of a 3-m-long inclined channel, whose bottom is constituted by an upward-moving conveyor belt (figure 7a). Owing to this mobile bottom, gravity-driven free-surface surges that remain stationary in the laboratory reference frame, can be generated. When a given volume of fluid is placed in the channel, flow height spontaneously adapts to the fixed bottom velocity u b [START_REF] Chambon | Gravity-driven surges of a viscoplastic fluid: An experimental study[END_REF][START_REF] Chambon | Experimental investigation of viscoplastic free-surface flows in steady uniform regime[END_REF]. The studied surges have a typical length of 1 m and a typical height of 10-30 mm, thereby complying with thin-layer assumption (see figure 8). Side-wall effects remain limited in this setup, such that transverse velocities are negligible and the flow can be regarded as effectively 2D along the centerline of the surges (except in the close vicinity of the tip: see [START_REF] Freydier | Experimental characterization of velocity fields within the front of viscoplastic surges down an incline[END_REF].

Viscoplastic experiments were performed with Carbopol (EDT 2623). Even if it can display more complex rheological trends in transient flows, this polymeric micro-gel was used as a model viscoplastic fluid in numerous previous studies (e.g., [START_REF] Piau | Carbopol gels: Elastoviscoplastic and slippery glasses made of individual swollen sponges. Meso-and macroscopic properties, constitutive equations and scaling laws[END_REF][START_REF] Bonn | Yield stress materials in soft condensed matter[END_REF]. As shown in figure 7b, the steady-state flow curve of the material is well represented by Herschel-Bulkley law. The yield stress τ c , consistency K, and power index n of the samples were determined using a laboratory rheometer (Bohlin CVOR) equipped with roughened parallel-plate tools to avoid slip. Values of n systematically lie around 0.4 (see table 1). Following our previous works [START_REF] Chambon | Experimental investigation of viscoplastic free-surface flows in steady uniform regime[END_REF][START_REF] Freydier | Experimental characterization of velocity fields within the front of viscoplastic surges down an incline[END_REF], measured values of τ c and K were increased by 11% and 19%, respectively, to ensure quantitative agreement between predictions and measurements of flow height H N in the uniform zone of the surges, far from the tip.

A few benchmark experiments were also performed with a Newtonian fluid (glucose syrup). In this case, the viscosity η of the samples was measured with a cone-plate geometry (figure 7b). Since glucose viscosity strongly depends on temperature [START_REF] Schellart | Rheology and density of glucose syrup and honey: Determining their suitability for usage in analogue and fluid dynamic models of geological processes[END_REF], particular care was taken to perform experiments as rapidly as possible, and avoid temperature variations, when working with this material. Results obtained with Newtonian surges are summarised in Appendix B.

Flows are monitored by illuminating the centerplane of the surges with a laser sheet, and acquiring images from the side with a 500 Hz camera. The transparent fluids are seeded with 10 µm tracers. We reconstruct velocity fields (see figure 8) through a spacetime analysis technique that takes full advantage of the stationary character of the surges [START_REF] Freydier | Experimental characterization of velocity fields within the front of viscoplastic surges down an incline[END_REF]. This technique achieves higher spatial resolutions than more classical correlation methods (PIV), with typical values in our experiments of 1.24 mm in the longitudinal direction and 0.62 mm in the vertical direction In addition, local flow height and free-surface shape are also extracted from the images, with an uncertainty in the range ±0.2 -0.9 mm.

An example of measured longitudinal velocity field is shown in figure 8. Far from the tip, the velocity profiles present the typical 2-layer shape expected in steady uniform regime, with an unsheared plug close to the free surface and a sheared layer underneath. When approaching the tip, velocity profiles begin to display shear throughout the entire depth of the flow. Accordingly, the unsheared zone seemingly disappears for x f /H N < 1 -2, typically [START_REF] Freydier | Experimental characterization of velocity fields within the front of viscoplastic surges down an incline[END_REF]. Experimental measurements also feature an increase in surface velocity when approaching the tip (figure 8).

Strategy of experimental-theoretical comparisons

For this study, results coming from fourteen different experiments with Carbopol, and two different experiments with glucose, were aggregated (see table 1). These experiments were performed with different bottom velocities u b (in the range 25 -150 mm.s -1 ), different slope angles θ (11.9 • and 15.3 • ), and different batches of Carbopol and glucose syrup. The yield stress τ c of Carbopol samples varied between 5 and 12 Pa. Accordingly, values of Bi and Re achieved in the viscoplastic experiments spanned the ranges 0.4 -1.1 and 0.04-1.6, respectively. Note that, in practice, numbers Bi and Re tend to be inversely correlated in the experiments, such that high (resp. low) values of Bi are generally associated to low (resp. high) values of Re (table 1). For the Newtonian experiments, the viscosity of glucose samples varied between 5 and 8 Pa.s, leading to values of Re in the range 0.4 -1.4. Overall, the values of ǫRe/Fr 2 are comprised between 0.1 and In the reference frame attached to channel bottom, the stationary surges generated in our setup correspond to travelling waves propagating at belt velocity u b . Direct comparisons with theoretical predictions of models M 1 and M 1+ǫ obtained in §3, can thus be performed. To that aim, experimental measurements of longitudinal velocity u are shifted by the value of u b . Due to an experimental artefact, the depth-averaged velocities ū obtained by integration of the shifted velocity profiles may differ from the expected value u b by a few percents [START_REF] Freydier | Experimental characterization of velocity fields within the front of viscoplastic surges down an incline[END_REF]). To normalise experimental data consistently with theoretical predictions, measured velocities are thus rescaled by the actual, experimental values of ū (constant along the surges). Similarly, positions x f and y are rescaled by the value of uniform height H N directly measured on the surges, far from the tip.

Due to the nonphysical kink at the fake yield surface in the theoretical velocity profiles at order O(ǫ) (see §3.2), detailed comparisons with the experimental velocity profiles were not attempted in the viscoplastic case. Comparisons presented hereafter shall thus concentrate on (i) surface velocity u s , and (ii) vertically-averaged shear rate in the pseudoplug layer | γ| pl . Note that we expect the evolution of | γ| pl along the surge to be essentially unaffected by the singularity at y = h -h p , since this quantity only depends on the longitudinal derivative

∂ x u (0) pl [see (3.6)].

Free-surface velocity

Figure 9 shows comparisons between experimental surface velocities and predictions of models M 1 and M 1+ǫ . Note that raw u s -measurements generally present a slight systematic discrepancy with theoretical predictions far from the tip (see insets). This discrepancy, comprised within experimental errors, is attributed to uncertainties on the rheological parameters and, thus, on experimental Bi values (see figure 4, in which the dependence of u s -predictions upon Bi is clearly observed). To better compare the modelled and experimental trends in the vicinity of the tip, datapoints in the main panel of figure 9 have been slightly shifted, so as to match theoretical predictions in the uniform zone.

As explained in §3, differences between predictions of models M 1 and M 1+ǫ get more pronounced when Bingham number Bi increases. Hence, let us first focus on an experiment characterised by a relatively large value of Bi and a small value of Re (figure 9a). In this case, it is clearly observed that experimental values of u s increase faster than predictions of model M 1 when approaching the tip, with discrepancies becoming apparent for x f < 2 typically. On the contrary, experimental measurements are in good agreement with model M 1+ǫ , even for the lowest values of x f at which measurements are available. The magnitudes of the O(ǫ) normal-stress and inertia correction terms are also represented in the figure. As expected, for this experiment, the O(ǫ) correction is dominated by the normal-stress term, while the inertia term remains much smaller and close to zero.

The better agreement provided by model M 1+ǫ close to the tip is also visible, though to a lesser extent, in figure 9b, which corresponds to an experiment with a lower value of Bi and a larger value of Re. Here, the magnitude of the inertia correction term is larger than in figure 9a, such that the net difference between models M 1 and M 1+ǫ is less pronounced. Experimental data could be viewed as agreeing with either one of the two models, although the predictions of model M 1+ǫ again appear to better capture the increasing trend close to the tip.

Finally, figure 9c shows the case of a large Re number, for which the inertia and normalstress correction terms are almost of the same order. Here, predictions of the M 1 and M 1+ǫ models are essentially indistinguishable. Experimental data, however, deviate from these predictions close to the tip. Note that measurements seem to follow the evolution of the normal-stress correction term, suggesting that the inertia correction might be overestimated in this case.

To get a synthetic view of all experiments performed, the following error ratio was computed: denote the surface velocities predicted by models M 1 and M 1+ǫ , respectively. Figure 10 shows values of r E obtained at x f = 0.7 as a function of Bi and Re numbers. Similar trends are obtained for other values of x f , as long as positions sufficiently close to the tip are considered. First, it is observed that r E displays values lower than 1, indicating that the match provided by model M 1+ǫ systematically better than that obtained with model M 1 . Moreover, the larger the values of Bi (or, correlatively, the lower the values of Re), the closer to 0 the values of r E . This reflects an excellent agreement between experimental data and predictions of model M 1+ǫ for large Bi values. On the contrary, for low values of Bi (or large values of Re), the ratio r E remains close to 1, meaning that the improvement brought by the O(ǫ) correction terms is relatively marginal in this case. In this limit of low Bi values, experimental data show significant discrepancies with both M 1 and M 1+ǫ models.

r E = |u exp s -u M1+ǫ s | |u exp s -u M1 s | , ( 4

Pseudo-plug shear rate

Figure 11 shows examples of shear rate | γ| computed by direct differentiation of the measured velocity fields (accounting for vertical and longitudinal derivative terms). To improve clarity, values of | γ| are normalised by the typical level of experimental noise γc (variable among the experiments), and the colour scale is thresholded for | γ|/ γc < 1. Hence, dark-blue/black areas correspond to the apparent unsheared layer that can be determined from the experiments. Far from the tip, this apparent unsheared layer is thicker than the predicted pseudo-plug, due to the thresholding at γc . In the tip region, on the contrary, it is observed that the experimental unsheared layer thins much faster than the pseudo-plug, and disappears for x f < 1 -2 typically. Interestingly, this feature appears to be consistent with the abrupt collapse of the theoretical effective unsheared layer obtained from model M 1+ǫ (see §3.3). If this theoretical effective unsheared layer is defined on the base of the iso-value of shear rate corresponding to measurement noise γc , both the overall shape of the layer and the location of the collapse appear in reasonable agreement with experimental data (figure 11).

More quantitatively, figure 12 compares the evolution of the vertically-averaged shear rate | γ| pl measured in the pseudo-plug with predictions of model M 1+ǫ . As for their theoretical counterparts (see §3.3), experimental values of | γ| pl are computed by averaging all | γ|-measurements obtained within the boundaries of the theoretical pseudo-plug for a given position x f . In agreement with the model, a clear increase of | γ| pl in the vicinity of the tip is visible in the data. Measurements appear to be nicely captured by the theoretical predictions for small to moderate values of Re (figure 12a and12b). Note that the slightly less compelling agreement observed for the largest value of Bi (figure 12a), can be explained by a larger noise level for this experiment. Conversely, as for surface velocities, predictions of model M 1+ǫ appear to underestimate the data for large values of Re (figure 12c).

Discussion

To the best of our knowledge, this paper presents the first complete asymptotic expansion, up to order O(ǫ), of longitudinal velocity u and flow discharge q in a thin free-surface viscoplastic surge. Our analysis expands on the previous work of [START_REF] Fernandez-Nieto | Shallow Water equations for non-Newtonian fluids[END_REF] in two respects: (i) the general case of a Herschel-Bulkley constitutive law (with n > 0) is treated , and (ii) we consider a flow regime in which ǫRe/Fr 2 = O(1), whereas this term was regarded as O(ǫ) by [START_REF] Fernandez-Nieto | Shallow Water equations for non-Newtonian fluids[END_REF]. This latter assumption provides more flexibility to the leading-order solution, which coincides with a lubrication approximation characterised by a variable pseudo-plug thickness h p depending on ∂ x h. Such a feature is required to model, e.g., surge fronts, in which pseudo-plug thickness cannot be regarded as constant. This flow regime also matches that typically achieved in the experiments used for comparisons in §4. As a counterpart, the assumption ǫRe/Fr 2 = O(1) leads to the appearance of factors proportional to ∂ x h p and ∂ t h p in the O(ǫ) correction terms. It is straightforward to check that the the expressions of correction terms derived in §2, effectively reduce to those obtained by [START_REF] Fernandez-Nieto | Shallow Water equations for non-Newtonian fluids[END_REF] when setting ǫRe/Fr 2 = 0 (and thus replacing Λ by λ and cancelling out all terms in ∂ x,t Λ).

Influence of O(ǫ) correction terms

Two sets of O(ǫ) correction terms, with different physical origins, affect the asymptotic expansions of u and q. Correction terms related to inertia effects also arise in the case of Newtonian or power-law fluids [START_REF] Fernandez-Nieto | Shallow Water equations for non-Newtonian fluids[END_REF][START_REF] Noble | Thin power-law film flow down an inclined plane: Consistent shallow-water models and stability under large-scale perturbations[END_REF][START_REF] Boutounet | Multi-regime shallow free surface laminar flow models for quasi-Newtonian fluids[END_REF]). Normal-stress correction terms, on the contrary, are specific to the viscoplastic character of the fluid, since they are directly related to the existence of the pseudo-plug layer. Note that in the flow regime considered here, and unlike in the work of [START_REF] Fernandez-Nieto | Shallow Water equations for non-Newtonian fluids[END_REF], these normal-stress corrections affect the velocity profile in both the pseudo-plug and the sheared layer underneath.

The relative influence of these O(ǫ) correction terms was investigated in the specific case of travelling-wave solutions, focusing on the tip region. The free-surface shape was shown to be only weakly sensitive to O(ǫ) corrections, regardless of the considered Bi and Re values (see figure 2). Interestingly, on the contrary, a significant effect of both inertia and normal-stress correction terms was observed on the velocity profiles, all the more for small values of n. This markedly contrasts with the case of Newtonian fluids (see Appendix B) or power-law fluids (not shown), for which leading-order and O(ǫ) velocity profiles remain essentially indiscernible, even close to the tip. In fact, an increase of Bingham number Bi was shown to result in an increased influence of both inertia and normal-stress correction terms (see figure 3). Even though they arise from different physical origins, these correction terms are strongly coupled, and both increase with the viscoplastic character of the fluid. It was also observed that the two correction terms generally act in opposite directions, such that their effects tend to counterbalance for sufficiently large values of Bi and Re. As a consequence, the net difference between predictions of leading-order and O(ǫ) models, varies in a complex and non-monotonic manner as a function of Bi and Re.

A shortcoming of the present study lies in the non-differentiability across the fake yield surface of the velocity profiles predicted at O(ǫ) (see figure 3). This singularity arises from the normal-stress correction terms in the pseudo-plug, more precisely from the factor [1 -((h -y)/h p ) 2 ] 1/2 in (2.34). Note that this factor is independent of n, such that the same singularity is also present in past studies considering Bingham fluids [START_REF] Balmforth | A consistent thin-layer theory for Bingham plastics[END_REF][START_REF] Fernandez-Nieto | Shallow Water equations for non-Newtonian fluids[END_REF]. As discussed in these previous studies, construction of more regular velocity profiles would require a proper asymptotic matching, through the introduction of a transition layer, between the sheared layer and the pseudo-plug. The complete formalisation of this asymptotic matching for general values of n, appears a promising prospect to further extend the results presented in this paper, and shall constitute the subject of future work.

Agreement with experimental data

Regardless of the number of terms retained in the asymptotics, the validity of the thin-layer assumption is clearly expected to break down in the vicinity of surge tip, where |∂ x h| → ∞. Note that, in absence of additional regularising conditions, even the primitive Cauchy equations become ill-posed at moving contact lines (e.g., [START_REF] Dussan | The moving contact line: the slip boundary condition[END_REF]. Nevertheless, since |∂ x h| typically remains less than 1 for x f > 0.1, well-ordered asymptotic expansions, and meaningful comparisons with experimental data, could be obtained relatively close to the tip. It ought to be mentioned here that, as shown in figure 13, a slight but systematic discrepancy is observed between measured and G. Chambon, P. Freydier, M. Naaim and J.-P. Vila 1) and predictions of models M1 and M1+ǫ. For the sake of clarity, only the two model predictions that differ most from one another, are shown. theoretical free-surface shapes for x f < 2, i.e. when |∂ x h| > 0.1, typically (see also [START_REF] Freydier | Experimental characterization of velocity fields within the front of viscoplastic surges down an incline[END_REF]. Arguably, this discrepancy could be interpreted as the signature of an early break-down of the thin-layer assumption. Let us however note that, in the Newtonian case, a much better agreement is recovered between measured and predicted free-surface shapes, even in the zone where |∂ x h| > 0.1 (see Appendix B). As discussed in [START_REF] Freydier | Experimental characterization of velocity fields within the front of viscoplastic surges down an incline[END_REF], this discrepancy on free-surface shapes observed with Carbopol, might also be due to more complex rheological trends of the material that are not accounted for by Herschel-Bulkley law.

Theoretical velocity profiles used for comparisons with experimental data were constructed from the modelled free-surface shapes. Using instead measured values of h, would have resulted in predicted velocity profiles characterised by non-constant depthaverages ū in the tip region, and thus more difficult to compare with experimental data. As already reported by [START_REF] Freydier | Experimental characterization of velocity fields within the front of viscoplastic surges down an incline[END_REF], for viscoplastic fluids, the leading-order approximation appears to be in excellent agreement with longitudinal velocity measurements for x f > 1 -2, typically. Closer to the tip, clear discrepancies are visible, notably a rapid collapse of the apparent unsheared layer identified in the experiments. In contrast, in the Newtonian case, the agreement between experimental velocity measurements and the leading-order approximation remains excellent all along the surge, including close to the tip (Appendix B). Note that similar comparisons between lubrication models and experiments were attempted by [START_REF] Andreini | Internal dynamics of Newtonian and viscoplastic fluid avalanches down a sloping bed[END_REF], although with less clear conclusions due to somewhat noisier data regarding the evolution of the unsheared zone in the tip, and strong discrepancies in average velocity.

The main result of the present study is the clearly improved match with experimental data obtained, for viscoplastic fluids, when O(ǫ) correction terms are included in the expansion. Importantly, the O(ǫ) approximation captures the abrupt collapse of the unsheared layer in the tip (see figure 11). The increase in shear rate and surface velocity are also quantitatively reproduced by the O(ǫ) model, and this down to values of x f as small as 0.1 -0.2 (see figures 12 and 9). Note that this good agreement regarding velocity and shear-rate data, tends to confirm that discrepancies observed on free-surface shapes for x f < 1 -2 (figure 13) are not imputable to a break-down of the asymptotics. The O(ǫ) model shows nevertheless a tendency to underestimate both velocity and shearrate measurements for Re > 1 -1.5, typically, as if the counterbalancing effect between inertia and normal-stress correction terms was overly pronounced in this case. Overall, the normal-stress correction appears to play a critical role in the improvement brought by the O(ǫ) approximation. Conversely, we suspect the magnitude of the inertia correction to be overestimated in our approach. It is however likely that the predictive capabilities of the O(ǫ) approximation would also be further improved with a more accurate free-surface model.

Conclusions

This paper presents consistent asymptotic expansions, up to O(ǫ) in flow aspect ratio ǫ, of the velocity field in free-surface surges of viscoplastic fluids. At leading-order, the expansion coincides with a lubrication approximation. The influence of O(ǫ) correction terms was investigated in the case of travelling-wave solutions, and detailed comparisons with the unique, high-resolution, experimental dataset of [START_REF] Freydier | Experimental characterization of velocity fields within the front of viscoplastic surges down an incline[END_REF] were drawn. The main conclusions of the study can be summarised as follows:

(i) While leading-order and O(ǫ) approximations are essentially indistinguishable for Newtonian and power-law fluids, the influence of O(ǫ) correction terms is enhanced by the viscoplastic character of the fluid, notably for small values of power index n. Viscoplastic flows thus provide interesting case studies, in which the differences between leading-order and O(ǫ) asymptotic models are potentially amenable to experimental measurements.

(ii) For viscoplastic fluids, inclusion of the O(ǫ) correction terms in the asymptotic analysis, qualitatively alters the description of the flow in the tip region. While the leading-order model predicts the existence of an unsheared pseudo-plug zone all along the surge, O(ǫ) correction terms entail a significant increase of the predicted shear rate, and an effective elimination of the unsheared zone, in the tip region.

(iii) For viscoplastic fluids, remarkable agreement is observed between the O(ǫ) model and experimental measurements of surface velocity and shear rate, and this even in close vicinity of surge tip, down to x f /H N ≈ 0.1. On the contrary, the leading-order model fails to reproduce evolution of surface velocity and shear rate for x f /H N < 1 -2, typically. For Newtonian fluids, both leading-order and O(ǫ) models nicely capture experimental data for all x f /H N > 0.1.

(iv) A more regular matching between the sheared layer and the pseudo-plug would be needed to construct differentiable velocity profiles above leading order in the viscoplastic case. This improvement shall allow us to achieve direct comparisons with measured velocity profiles.

To conclude, let us recall that the obtained O(ǫ) expansions of longitudinal velocity u and discharge q, constitute pre-requisites for the formulation of consistent closure terms in shallow-water models (see Introduction). Beyond the argument of mathematical consistency, the results obtained in this study also suggest that, for viscoplastic fluids, inclusion of O(ǫ) terms might have a noticeable impact on the predictions of these models in terms of, e.g., velocity values. Proper derivation of consistent shallow-water models for viscoplastic fluids, and comparisons thereof with experimental data, shall be the subject of upcoming studies. Note however that the single-equation model M 1+ǫ obtained in §3.1, can already be regarded as a specialised version of a consistent shallow-water model adapted to the case of travelling-wave solutions (for which q is enslaved to h). Let us also mention that alternative expressions of the O(ǫ) correction terms, which are all equivalent up to order O(ǫ), could be proposed. Accordingly, a whole family of shallow-water models can actually be constructed (see also [START_REF] Boutounet | Multi-regime shallow free surface laminar flow models for quasi-Newtonian fluids[END_REF]. In this work, differential identity (A 1) was used to retain only factors in ∂ x h and ∂ x,t Λ in the final expressions of velocity profiles at O(ǫ). The exact influence of these algebraic manipulations on the final predictions of shallow-water models, remains to be properly investigated. where we recall that x f = t -x with the considered choice of characteristic scales. With the boundary condition h tw (0) = 0, this equation admits an implicit analytical solution:

x f = cot θ 1 2 ln 1 + h tw 1 -h tw -h tw . (B 8)
Hence, the normalised free surface predicted by model M 1 is completely independent of fluid rheological parameters in the Newtonian case. Model M 1+ǫ is based on the O(ǫ) approximation of discharge q (0) +q (1) , with q (1) given by (B 5). The corresponding secondorder ODE takes here the following form: and Λ tw = 3 (1 + cot θ dh tw /dx f ).

M 1+ǫ :

B.3. Model predictions

Figure 14 shows that free surfaces computed with models M 1 and M 1+ǫ are essentially indistinguishable in the Newtonian case, and this regardless of the value of Re (in the range 0 -1.5). Concerning free-surface velocity u s , a constant value of u s is predicted by model M 1 , while values of u s computed with model M 1+ǫ display a decreasing trend towards the tip (figure 15). Note however that the differences between the two models remain typically tenfold smaller than in the viscoplastic case (for similar values of Re: see §3.2). As in figure 4, u s -curves in figure 15 are only represented for x f > 0.1, i.e. in the domain where |∂ x h| remains typically less than 1 (see figure 14).

B.4. Comparison with experimental data

Characteristics of the experiments performed with Newtonian samples (glucose syrup) are summarised in table 1. As shown in figure 16, experimentally-measured front shapes are in good agreement with the theoretical predictions provided by either model M 1 or M 1+ǫ . The two models tend here to slightly overestimate the data, but the differences remain significantly smaller than in the viscoplastic case (compare with figure 13).

Figure 1 .

 1 Figure 1. Free-surface viscoplastic surge down an incline: flow configuration considered in this study (see notations in text). Dashed line represents the fake yield surface.

Figure 2 .

 2 Figure 2. Free-surface shapes predicted by models M1 and M1+ǫ: evolution of flow height h and derivative ∂xh versus distance from the tip x f . (a) n = 0.4. (b) n = 1. In both plots, solutions computed for two values of Bi and two values of Re are shown (θ = 15 • ). Insets show close-ups on the slight differences among the solutions.

Figure 3 .

 3 Figure 3. (a) and (b) Vertical profiles of longitudinal velocity u computed with models M1 and M1+ǫ at two distances from the tip x f (see legend). In both plots, solutions computed for two values of Bi and two values of Re are shown (θ = 15 • , n = 0.4). (c) Vertical profiles of the O(ǫ) normal-stress (thick lines with squares) and inertia (thin lines with circles) correction terms u (1) N

  two correction terms generally display opposite signs (except for very low values of Bi ). The normal-stress term u(1)

Figure 4 .

 4 Figure 4. Evolution of free-surface velocity us predicted by models M1 and M1+ǫ as a function of distance from the tip x f . (a) n = 0.4. (b) n = 1. In both plots, solutions computed for two values of Bi and two values of Re are shown (θ = 15 • ). In plot (a), the red star marks the position where the asymptotic expansion in model M1+ǫ ceases to be well-ordered for the case Bi = 1.5, Re = 0.1 (blue dot-dashed line).

  shown in figure 5, values of | γ|(1)

Figure 5 .

 5 Figure 5. Evolution of O(ǫ) vertically-averaged shear rate in the pseudo-plug | γ| (1) pl predicted by model M1+ǫ, as a function of distance from the tip x f . (a) n = 0.4. (b) n = 1. In both plots, solutions computed for two values of Bi and two values of Re are shown (θ = 15 • ). The dotted line indicates the 0 value (i.e., M1 solution).

Figure 6 .

 6 Figure6. Iso-contours of O(ǫ) shear rate | γ|(1) in the pseudo-plug predicted by model M1+ǫ (dotted lines). The solid line represents the free-surface and the dashed line represents the fake yield surface y = h -hp (Bi = 1.5, Re = 1.5, θ = 15 • , n = 0.4).

Figure 7 .

 7 Figure 7. (a) Scheme of the conveyor-belt channel. The stationary surge is represented in white. (b) Typical steady-state flow curves of Carbopol and glucose samples used in the experiments: shear stress τ versus shear rate γ (semi-log scale). Symbols correspond to rheometrical measurements; dashed lines represent the best-fitting Herschel Bulkley and Newtonian laws (Carbopol: τc = 9.6 Pa, K = 9.0 Pa.s n , n = 0.42; Glucose: η = 5.6 Pa.s).

Figure 8 .

 8 Figure 8. Measured longitudinal velocity (u) field for experiment HB13 (Bi = 0.5, Re = 1.4). The colour scale corresponds to raw velocity values (no shifting by bed velocity u b ). Streamlines are superimposed in white. Two vertical profiles of u, obtained at positions x f = 95 mm and x f = 25 mm, respectively, are also represented. The grey area close to the tip indicates the zone where no reliable velocity measurements could be obtained.

Figure 9 .

 9 Figure 9. Free-surface velocity us as a function of distance from the tip x f : comparison between experimental data and predictions of models M1 and M1+ǫ. Error bars correspond to velocity measurement uncertainties. Light blue dotted and dot-dashed lines show contributions of the O(ǫ) normal-stress and inertia correction terms u (1) s,N and u (1) s,I , respectively (scale on right axis). Insets display actual experimental data, while main panels show data corrected by a slight vertical shift (see text). (a) Bi = 1.1, Re = 0.04 (exp. HB1). (b) Bi = 0.8 and Re = 0.2 (exp. HB4). (c) Bi = 0.5 and Re = 1.4 (exp. HB13).

Figure 10 .Figure 11 .

 1011 Figure 10. Error ratio rE (see text) computed at x f = 0.7 as a function of (a) Bingham number Bi and (b) Reynolds number Re, for the 14 experiments performed with Carbopol (see table1).

Figure 12 .

 12 Figure 12. Vertically-averaged shear rate in the pseudo-plug | γ| pl as a function of distance from the tip x f : comparison between experimental data and predictions of model M1+ǫ. The grey dotted line indicates the typical level of measurement noise γc. (a) Bi = 0.95, Re = 0.06 (exp. HB2). (b) Bi = 0.7, Re = 0.3 (exp. HB7). (c) Bi = 0.5, Re = 1.4 (exp. HB13).

Figure 13 .

 13 Figure13. Evolution of flow height h and free-surface slope ∂xh as a function of distance from the tip x f : comparison between all experimental data obtained for θ = 15.3 • (see table1) and predictions of models M1 and M1+ǫ. For the sake of clarity, only the two model predictions that differ most from one another, are shown.

Figure 14 .

 14 Figure 14. Free-surface shapes predicted by models M1 and M1+ǫ for a Newtonian fluid: evolution of flow height h and derivative ∂xh versus distance from the tip x f . Solutions of model M1+ǫ computed for three different values of Re are shown (θ = 15 • ). Inset shows a close-up on the minute differences among the solutions.

  These quantities are easily derived from (2.22) and mass conservation (2.3), and we shall not provide the full expressions in this paper (see[START_REF] Freydier | Experimental characterization of velocity fields within the front of viscoplastic surges down an incline[END_REF], for expression of v (0) ). Let us nevertheless mention that, compared to the developments of Fernandez-Nieto et al. (2010), additional complexities arise here from terms proportional to the derivatives ∂ x,t h p and ∂ x,t Λ. These derivatives can be interchanged through the following relation:

		∂ x,t h p = -	h p Λ	∂ x,t Λ.	(2.25)
		2.4. Pseudo-plug
	The dependence upon x of u	(0)

Table 1 .

 1 Characteristics

	Carbopol
	exp. # τc	K	n	θ	u b HN Bi	Re	Fr Bi /λ ǫRe/Fr 2
	HB1 13.3 9.8 0.40 15.3 22.6 14.6 1.13 0.044 0.061 0.35	0.17
	HB2 10.7 10.5 0.42 11.9 26.0 17.0 0.95 0.060 0.064 0.31	0.25
	HB3 13.3 9.8 0.40 15.3 37.0 16.2 0.95 0.098 0.095 0.32	0.18
	HB4 13.3 9.8 0.40 15.3 53.0 17.5 0.84	0.18	0.13 0.29	0.18
	HB5 9.0 6.9 0.42 15.3 46.0 12.4 0.80	0.19	0.13 0.28	0.13
	HB6 9.0 6.9 0.42 15.3 64.0 13.6 0.70	0.32	0.18 0.26	0.14
	HB7 13.3 9.8 0.40 15.3 78.0 19.0 0.74	0.34	0.18 0.27	0.19
	HB8 8.0 6.1 0.41 11.9 73.0 16.1 0.65	0.44	0.19 0.25	0.20
	HB9 13.3 9.8 0.40 15.3 128 21.6 0.62	0.76	0.28 0.24	0.21
	HB10 10.7 10.5 0.42 11.9 153 26.4 0.50	1.09	0.30 0.20	0.31
	HB11 9.0 6.9 0.42 15.3 123 15.9 0.56	0.95	0.32 0.22	0.15
	HB12 5.6 5.7 0.41 11.9 151 16.2 0.40	1.64	0.38 0.17	0.18
	HB13 8.0 6.0 0.43 15.3 148 15.5 0.50	1.37	0.39 0.20	0.14
	HB14 9.0 6.9 0.42 15.3 171 17.2 0.51	1.64	0.42 0.20	0.16
	Glucose syrup
	exp. # ρ	η	θ	u b HN Re Fr ǫRe/Fr 2
	Nt1 1383 5.6 15.3 75 18.6 0.35 0.18	0.20
	Nt2 1380 4.9 11.9 164 28.1 1.42 0.32	0.40

of the experiments performed with Carbopol and glucose syrup: fluid density ρ (kg.m -3 : only indicated for glucose, taken at 1000 kg.m -3 for Carbopol), yield stress τc (Pa; for Carbopol), consistency K (Pa.s n ; for Carbopol), power-law index n (for Carbopol), viscosity η (Pa.s; for glucose), slope angle ( • ), bottom belt velocity u b (mm.s -1 ), steady uniform height HN (mm), Bingham number Bi (for Carbopol), Reynolds number Re, Froude number Fr , plug thickness in uniform regions Bi/λ (for Carbopol), and number ǫRe/Fr 2 . 0.4. Although somewhat smaller than 1, these values are systematically larger than ǫ (calculated for a typical length of 1 m) by at least one order of magnitude, indicating that the experiments effectively comply with the flow regime considered in this study (see §2.2). Lastly, for reference, values of Fr and Bi /λ (plug thickness in uniform regions) are also indicated in table 1.

† Email address for correspondence: guillaume.chambon@irstea.fr

Acknowledgements

Financial support from AGIR PEPS project MIGRAIN (Univ. Grenoble Alpes) and CNRS (InFIniti program) is acknowledged. IRSTEA is member of Labex TEC21 (Investissements d'Avenir -Grant Agreement ANR-11-LABX-0030) and Labex Osug@2020 (Investissements d'Avenir -Grant Agreement ANR-10-LABX-0056).

Declaration of Interests. None.

Appendix A. Expressions of U I and Q I In this appendix, we provide explicit expressions for the inertia-related correction functions U I and Q I involved in the expansions of longitudinal velocity u and flow discharge q at O(ǫ) (see §2.5). As a preamble, let us note that a useful differential identity, valid at order O(1), derives from the depth-integrated mass conservation equation (2.16). Inserting the leading-order approximation of discharge q (0) given by (2.24), leads to: The function U I is obtained from (2.44) and the expressions of leading-order longitudinal and vertical velocities u (0) and v (0) . The lengthy algebraic calculations were performed with wxMaxima symbolic software. The final expression writes:

where q U ,1 (ζ, ̟ p ) and q U ,2 (ζ, ̟ p ) are polynomials in ζ and ̟ p :

It can be checked that q U ,1 (1, ̟ p ) = q U ,2 (1, ̟ p ) = 0, as imposed by the no-slip basal boundary condition.

The function Q I is calculated by a further integration of U I :

. The final expression writes:

Appendix B. Newtonian case

This appendix presents: (i) the specialised version of the asymptotic expansion derived in §2 for the case of Newtonian fluids, and (ii) comparisons with the experimental results obtained with glucose syrup (see table 1).

B.1. Velocity profile and discharge

Setting Bi = 0 (thus h p = 0) and n = 1 in (2.22), a classical parabolic profile is recovered for longitudinal velocity u at leading order:

with ζ = 1 -y/h. Accordingly, leading-order discharge q (0) writes:

At O(ǫ), only inertia correction terms are present, such that (2.46) reduces to:

with here:

4) Similarly, the O(ǫ) correction on discharge writes:

with

. Travelling-wave solutions

The same choice of characteristic scales as in §3, results here in the identity λ = 3. Like for the viscoplastic case, two free-surface models can be constructed from the relation q = h tw characterising travelling-wave solutions. Model M 1 is based on the leading-order approximation (B 2) of q. This model writes here: Similarly, as shown in figure 17, measured velocity profiles are remarkably captured by both asymptotic models. At the scale used in the figure, the velocity profiles predicted by models M 1 and M 1+ǫ are indistinguishable. The agreement with experimental data remains excellent even for the smallest values of x f investigated (figures 17c and 17e).