Towards a new friction model for shallow water equations through an interactive viscous layer

The derivation of shallow water models from Navier-Stokes equations is revisited yielding a class of two-layer shallow water models.An improved velocity profile is proposed, based on the superposition of an ideal fluid and a viscous layer inspired by the Interactive Boundary Layer interaction used in aeronautics. This leads to a new friction law which depends not only on velocity and depth but of the variations of velocity and thickness of boundary layer. The resulting system is an extended shallow water model consisting of three depth-integrated equations: the first two are mass and momentum conservation in which a slight correction on hydrostatic pressure has been made; the third one, known as von Karman equation, describes the evolution of the viscous layer. This coupled model is shown to be conditionally hyperbolic, and a Godunov-type finite volume scheme is also proposed. Several numerical examples are provided and compared to the"Multi-Layer Saint-Venant"model. They emphasize the ability of the model to deal with unsteady viscous effects. They illustrate also the phase-lag between friction and topography, and even recover possible reverse flows.


Introduction
Many phenomena in fluvial or maritime hydraulics involve free surface flows in shallow waters for the study e.g. of floods and tides. Shallow water equations were originally introduced by Saint-Venant in 1871 [12] in the context of channel modelling. Since then, the model has been widely extended and is used in the modelling and numerical simulation of a number of natural or man-made phenomena such as river flow [25,6], flood forecasting [7,33], pollutant transport [48,26], dam-break [1,53], tsunami [20,32,45], overland flow [16,52,13], soil erosion [8,40] and many others. The shallow water system can be derived from the incompressible Navier-Stokes equations under several hypotheses; the main one being the long wave approximation meaning that the characteristic wavelength is much larger than the water depth (see figure 1 for a sketch and definitions). Two consequences follow then: the hydrostatic pressure law, and the viscous term vanishing in the horizontal direction. Next, to proceed from Navier-Stokes to shallow water, the equations are integrated along the vertical direction. At this point, care has to be taken of the vertical velocity profile, which on the one hand has to be approximated to deal with nonlinearities of the momentum flux, but on the other hand drives the bottom boundary condition, hence the friction phenomena.
Two classical assumptions on the longitudinal velocity profile along the vertical direction lead to explicit integrations. The first one is a viscous Poiseuille-like (i.e. parabolic) profile on the whole water depth which gives rise to a linear (with respect to the depth-averaged velocity) friction term, sometimes referred to as laminar friction. The second one is a constant profile, somehow corresponding to an ideal fluid; but, by construction, there is a priori no friction term in the integrated equations. Friction has to be added afterwards using empirical laws such as Manning, Chézy, etc (see e.g. [11]). The main drawback of these classical approaches is the nonadaptability of the friction terms for large variations of velocity because the assumed profiles (parabolic or flat) do not hold.
We intend here to pay a particular attention to the fact that these empirical laws are unable to describe the fluid inertia or more precisely to predict the phase-lag between the bottom friction and a perturbation of the bed which is known as an essential mechanism for dune or ripple formation [30]. Indeed, it is well reported that for the case of a sub-critical flow on a bump, the maximum of the friction must be slightly shifted upstream of the crest to drag the particles from the troughs up to the crest. Consequently, coupling classical shallow water equations with a mass conservation equation for sediment transport (e.g. Exner equation [17] for bedload case) cannot predict the bed instability, see [10,34] for more details. We look for a more flexible model with a better understanding of how the no-slip boundary condition gives rise to the friction term in the depthintegrated equations. This will allow to recover this phase-lag phenomenon as well as boundary layer separation, a manifestation of the recirculation of the flow near bottom. This is done by using an asymptotic description of the fluid as a superposition of an ideal fluid over a viscous layer located at the bottom, with a strong interaction between the layers. The thickness of the viscous layer is quantified by a small parameterδ related to the inverse of the Reynolds number of the flow. Integrating the incompressibility equation under this consideration leads to the same mass conservation equation of the usual shallow water system. On the contrary, the integration of the momentum equation exhibits major differences. On the one hand, it turns out that the order of magnitude of the friction term is larger as expect: precisely of orderδ, while the above mentioned Poiseuille profile leads to aδ 2 order of magnitude. In the case of an ideal fluid (δ = 0), the model degenerates, of course, into the classical shallow water one. On the other hand, we introduce a new closure for the momentum flux which involves an additional pressure law of orderδ. At this stage, we obtain a system of two equations which are similar in structure with the usual shallow water system, but involving several additional unknown functions related to the viscous layer.
The next step consists therefore in a careful analysis of the viscous layer. Following a classical methodology in aerodynamics, see e.g. [49], we integrate the Prandtl equation along the vertical axis to obtain the so-called von Kármán equation. It describes the evolution of the so-called displacement thickness δ 1 (see Figure 1), which is involved in the definition of the afore mentioned corrective pressure, and can be interpreted as some physical thickness of the viscous layer. The system has to be complemented by the velocity equation of the ideal fluid, since it is involved in the von Kármán equation. We will discuss on the effects for flows over short bumps. The acceleration induced by the bump will change a lot the basic flow so that the shape velocity profile is no longer a half-Poiseuille nor a flat one. This study aims to understand this kind of flows which are not taken into account by the shallow water equations themselves. Furthermore, we will present the link between our system and the Multi Layer Saint Venant model proposed in [2].
The outline of the paper is as follows. In a first section we recall the Navier-Stokes system, and state the long wave approximation which is convenient for shallow water approximation. Next we turn to the viscous layer analysis, and derive Prandtl and von Kármán equations. Velocity profiles are also introduced. The third section is devoted to the derivation of various formulations of the Extended Shallow Water System. In Section 4 we derive some formal properties of the model, together with the numerical scheme. Finally, we evidence several properties of the model by numerical simulations.

From Navier-Stokes to shallow water equations
In this section we recall how classical models for shallow waters are obtained from Navier-Stokes equations. The first assumption is a long wave approximation, stating that indeed we deal with a thin layer of water. Next, we integrate along the vertical direction, assuming a given velocity profile on the whole water depth. h 0 , the bottom is a given function f b of characteristic length L and η is the free surface. Two families of velocity profiles are displayed for the flow over the topography, first with the usual half-Poiseuille description (dashed), and second with the flat profile with a boundary layer (plain). Note that the shear (slope of the velocity at the wall) is completely different in those two descriptions even flux and depth are the same.

Navier-Stokes equations
We consider a fluid evolving in a time-dependant domain Ω t = R × {f b (x) ≤ y ≤ η(t, x)}. This thin layer is limited below by a fixed bottom represented by a function y = f b (x) and above by the free surface described by y = η(t, x). We denote the water depth h = η − f b . In this study, the properties of the air above the free surface are completely neglected (see Figure 1).
Our starting point is to consider the dimensionless Navier-Stokes equations expressing the mass and momentum conservations of an incompressible Newtonian fluid [49]. For the sake of simplicity, we limit ourselves in this work to consider only laminar flows. Indeed, the asymptotic from the Navier-Stokes equations is clearer and the resulting description is quantitative. A similar study for turbulent flows can be made with a modified Reynolds tensor. Nondimensionalization has been made with the same characteristic length h 0 , e.g. a reference water depth, for both the abscissa and the ordinate. The dimensionless Navier-Stokes equations write where u, v are the horizontal and vertical velocities respectively and p is the pressure. We have defined the Reynolds Re and Froude F r numbers given by in which u 0 and ν being the reference velocity and kinematic viscosity respectively and the constant g stands for the acceleration due to gravity. The Reynolds number expresses the ratio between the inertia force and the viscosity; the Froude number represents the ratio between the kinetic and potential energies.
The system is complemented with the following boundary conditions: • at the free surface y = η(t, x): kinematic boundary condition: v = ∂ t η + u∂ x η continuity of the stress tensor: is the stress tensor and n = 1 is the outer unit normal to the free surface.

Long wave scaling
Up to now, no hypothesis has been taken into account for the size order of the characteristic quantities u 0 , h 0 . We have in mind applications to rivers or coastal flows for which the following conditions may be observed: • the Reynolds number is large, • the horizontal velocity has small variation along the vertical, • the vertical velocity is small compared to the horizontal velocity.
We introduce the aspect ratio where L is a characteristic wave length and h 0 a characteristic depth. Let us start by investigating the third one, which justifies the following scaling for the velocities: v = εṽ, u =ũ.
Then the mass conservation equation (1.1) enforces also a scaling for the space variables y x since Hence there are two options for the variable scaling: 1. Long wave scaling: x =x ε , y =ỹ, t =t ε .
It means that the long wave hypothesis needs a long time study. Furthermore, since and so the bottom slope needs to be small enough.
2. Thin layer scaling: This scaling restricts the study of small vertical velocity only to a thin water depth which tends to zero when ε → 0. It is the classical scaling used in the boundary layer approach.
So, the long wave scaling is compatible with the hypothesis about small vertical velocity compared to horizontal velocity. Let us see the consequences for the set of equations (1.1)-(1.3) and the boundary conditions: Taking an approximation at order O(ε) leads to: • Cancellation of the viscosity in horizontal direction. Equation (1.4) reduces to: • Simplified version of the stress tensor continuity at the free surface: • Hydrostatic pressure with (1.5) andp = 0 at the surface: This approximation implies that ∂xp does not depend onȳ. In other words, the pressure gradient is conserved over the vertical. This result for pressure at order O(ε) is already observed (see e.g. [29]). Remark 1.1. We emphasize here that the above properties, which are classical shallow water hypotheses, are brought out solely by the long wave approximation.
To summarize, the long wave approximation of the Navier-Stokes equations consists of the following system, sometimes called the RNSP equations (Reduced Navier-Stokes/Prandtl [38]). This will be our reference system for the remaining of this article. In these equations, it is more convenient to define an effective Reynolds Re h number which takes into account the aspect ratio of the model Dropping all the tildes from the variables and unknowns, the system can be written: (1.11) In the limit Re h → ∞, the case of an incompressible ideal fluid, the RNSP equations degenerate to the hydrostatic Euler system. Indeed, the equations of mass conservation (1.6), hydrostatic pressure (1.8) and boundary condition at free surface (1.9)-(1.10) are unchanged, but the momentum conservation equation is replaced by We conclude this paragraph by giving the behaviour of horizontal velocity.
Proposition 1.2. The horizontal velocity u of a smooth solution of the hydrostatic Euler system is constant along the vertical direction, so we write u := u e (t, x), and verifies the following equation: Proof. Applying a partial derivative in y on equation (1.12) and taking into account the incompressibility (1.6) together with the hydrostatic pressure law (1.8) leads to d dt (∂ y u) = 0. It means that ∂ y u remains constant along the characteristic curve x (t) = u, y (t) = v, hence ∂ y u = ∂ y u| (x0,η0) = 0 by condition (1.10) where (x 0 , η 0 ) is the foot of characteristic curve starting at the free surface y = η. Consequently the horizontal velocity u of hydrostatic Euler system is independent on y.
As it is well-known in Ideal Fluid theory, the depth-independent horizontal velocity component leads to a slip velocity at the bed. Thus, the no-slip boundary condition (1.11) is not relevant for ideal fluid and has to be replaced by a weaker one called the non-penetration condition (1.14) Recovering some connection between the ideal fluid equations and the no-slip condition is precisely the aim of the viscous layer theory, which we will present in section 2.

Classical shallow water equations
Let us recall briefly the classical way to obtain the shallow water model by vertical integration the RNSP equations over the whole water depth. Integrating the mass conservation equation (1.6), from f b to η, yields in which we have applied the kinematic condition (1.9) at the surface and the no-slip condition (1.11)-or the weaker one (1.14)-at the bottom. As f b is time-independent, we can write ∂ t η = ∂ t (η − f b ); and by using definition (1.15), the mass conservation in its integrated form reads In the same way, we now derive the depth-integrated momentum balance equation. The main point to notice at this stage is the appearance of the nonlinear momentum flux and the friction term. Indeed, by integrating the momentum equation (1.7) over the depth and using the fact that ∂ x p does not depend on y, and also the condition ∂ y u| y=η = 0 from (1.10), we get Using definition (1.15) and applying the free surface condition (1.9) and the bottom condition (1.11) or (1.14), we can rewrite the integrated momentum conservation equation in the form in which we have introduced the so-called Boussinesq coefficient β and the bottom shear stress τ b , also called friction, which are defined by Therefore, evaluating β and τ b requires the knowledge of the flow. It can be checked that β ≥ 1 since, by definition (1.18), we can write Without complementary equations, a closure relation on the velocity profile is needed in order to compute the Boussinesq coefficient and to express the friction term in function of the conservative variables (h, hU ). Let us recall two constitutive profiles which are often adopted in the context of shallow water flows.
Flat profile. This is the most classical approach in hydraulic river modelling. Scaling analysis reveals that the velocity profile is quasi-flat except within a very thin-layer close to the river bed. Based on this consideration, the velocity profile can be assumed to be flat over the whole water depth, so that β equals one. The friction term τ b is not properly defined in this context. The flow can be described using the Euler system resulting in an inviscid shallow water model -equation (1.17) without the friction term τ b . A large family of empirical friction laws exist, which express the friction as a quadratic function of U with a friction coefficient C f = O(Re −1/4 ) for a smooth bottom, see [49]. This coefficient depends on h and U as well, for instance with Chézy, Manning laws, see [11] for a bibliographical study. In summary, this kind of model consists in writing (1.20) We emphasize again that the friction term derives here from empirical considerations.
Poiseuille profile. This type of profile is inspired from a analytic solution of the RNSP equations in the case of an uniform flow on a negative constant slope. The balance between the friction and the driving force of the slope gives a self-similar parabolic solution, also known as the half-Poiseuille or Nusselt solution: This choice of profile leads to The friction is indeed linear with respect to the mean velocity, and is referred to as laminar friction.
Using such a prescribed profile in shallow water equations leads to some important restriction of the model, in particular when dealing with large variation of the velocity. For example it has been reported in [28] that a constant value for the Boussinesq coefficient is less adapted to describe the dynamics of the fluid layer close to a dry-wet transition.
We conclude this section by evidencing that both these models do not present a phase-lag for the flow over a bump. Let us study a steady linearized solution of the usual shallow water equations. Consider such a small perturbation of the bed f b that we can write in the form f b = εf 1 b ,where ε is just a small parameter and not necessarily the aspect ratio defined before. We look for the solution in the form For high Reynolds numbers, we see from relations (1.20) and (1.21) that the friction is negligible. We can consider therefore h, U solution of the following frictionless and steady state shallow water equations (1.23) Inserting (1.22) in (1.23) and identifying powers of ε leads to a cascade of equations for each terms h i , U i , i = 0, 1. Then, it should be checked that the zero-order terms h 0 and U 0 are needed constant. For first-order terms, a straightforward calculation leads to By introducing the local Froude number F r 0 , we can express the linearized solution in the form (1.24) As we can see, U is exactly in-phase with f b . As a consequence, local maxima of the friction estimated by empirical formulas (1.20) or (1.21) are always reached at f b = 0, that is at the crest of the bump. Indeed,

Viscous layer analysis
We turn now to the main step towards the model we look for. It mainly consists in dividing the fluid in two layers: • an ideal fluid layer dealing with the free surface; • a thin viscous layer with the no-slip condition at the bottom.
In the first layer addressing to ideal fluid, we take advantage of the explicit integration along the vertical. In the second one describing viscous layer, we take into account the viscosity in the vertical direction and recover some friction in the integrated equations. This section is devoted to the study of the viscous layer, and to the analysis of the interactions between the two layers.
We introduce a small parameterδ, whose magnitude will be specified below. It is related to the thickness of the viscous layer, but does not correspond to its actual physical value. We follow the classical strategy used in the boundary layer theory [46,49] except that in that caseδ → 0, whereas we keep a finite value here. The first step is to rescale again the RNSP equations with the thin layer scaling to obtain a set of the well-known Prandtl equations. The next step consists in vertical integration of these equations over the viscous layer. This leads to the so-called von Kármán equation, where extra unknowns are introduced. Finally, some suitable assumptions have to be made on the velocity profile in order to obtain a closed model.

Prandtl equations
We introduce the following change of variables, referred to as the Prandtl shift: There are several possible scalings forδ in terms of Re h : • ifδ verifies Re hδ 2 1, we recover the ideal fluid equations; • ifδ satisfies Re hδ 2 1, we obtain ∂ 2 yū = 0 that leads toū = 0 due to the continuity of the stress tensor and the no-slip condition. So we do not consider this trivial case; • the last possibility is Re hδ 2 ∼ 1, which balances the convective terms and the diffusive one. It is called dominant balance or least degeneracy principle [54], and allows to preserve as many as possible terms in the equations. This is why in what follows, we consider the scalinḡ With this choice ofδ, the viscous term appears with the same order as the other terms in the momentum equation. We obtain the Prandtl equations written in viscous layer variables: We notice that this system of equations is in the same form as the Prandtl equations obtained directly from Navier-Stokes equations with classical boundary layer scaling ( [49], ch. VII) except for the topography term in the momentum equation (2.4).
Up to now, we do not have enough boundary conditions for the viscous layer. The natural connection consists in assuming that the velocity at the "top" of the viscous layer coincides with the velocity of inviscid layer. Precisely, we impose the following matching boundary condition which is obviously compatible with the Prandtl shift (2.1) sincex = x andt = t. Notice that in classical boundary layer theory the limit is given by the asymptotic matchingū(t,x,ȳ → ∞) → u e (t, x).

Von Kármán equation
The von Kármán equation expresses the defect of velocity between the ideal fluid and the viscous layer. A classical way to obtain such an equation consists in writing the Prandtl equation, then introducing the velocity defect (u e −ū), and finally integrating it on the viscous layer (see Schlichting [49]). We introduce the following two integrated quantities, see Figure 2: Let U be the depth-averaged velocity. We define • the displacement thickness δ 1 given by

8)
• the momentum thickness δ 2 given by Physically, the displacement thickness expresses the distance by which the ground should be displaced to obtain an ideal fluid with velocity u e and the same flow rate hU (see Figure 2). In the same way, the momentum thickness accounts for the loss of momentum in the viscous layer. A simple computation from (2.8) and (2.9) leads to the following expressions for these quantities: In the limitδ → 0, we recover the classical formulae for δ 1 , δ 2 in the boundary layer scaling [49]: Proposition 2.2. The evolution of the displacement and momentum thicknesses are ruled by the so-called von Kármán equation: Proof. First we notice that the Prandtl shift leads to following relations Momentum equation (1.13) for inviscid flow can be rewritten as in which we have used (2.5) to rewrite the right-hand side. The difference between this equation and (2.4) gives Through (2.3) and (2.6), we can rearrange the termv = − ȳ 0 ∂xū. Furthermore we get Using integration by parts, the last term in the left-hand side can be rewritten as Now we integrate the resulting equation overȳ, between 0 andη, together with the matching boundary condition (2.7) to obtain the momentum integral equation The last three terms of the left-hand side are now rewritten as ∂x η 0ū (ū e −ū) dȳ. Moreover, since u e is independent ofȳ, we have the relations Coupled with equation (1.13) on the velocity u e of ideal fluid, the von Kármán equation (2.10) gives only a partial representation of the boundary layer, since it involves four additional unknowns, namely u e , δ 1 , δ 2 , and τ b . To proceed further towards an integrated model, we need to specify velocity profiles to close the von Kármán equation.

Velocity profile in the viscous layer
Through the viscous layer, the velocityū varies from 0 (at the bottom) to the ideal fluid velocity u e . Therefore we introduce a profile function ϕ as well as a scaling factor ∆(t,x), chosen in such a way that ∆ quantifies the physical thickness of the viscous layer. Following [49], we wish to havē Therefore we choose 0 < ∆ ≤η, and a profile function ϕ(ξ) such that Hence, by definition of δ 1 and δ 2 , we can write To link these variables, we introduce the shape factor H which only depends on the profile function ϕ 14) The parietal constraints can also be expressed in terms of ϕ and u e as where the parameter f 2 is known as the friction factor (see [49,38]) and f 2 H := α 1 ϕ (0). Using definitions (2.14) and (2.15), we can rewrite the Von Kármán equation (2.10) in following form At this stage, choosing a velocity profile in the viscous layer amounts to impose a closure formula on the shape factor H and the friction factor f 2 . Once this is done, the von Kármán will be given a closed form, in terms of u e and the displacement thickness δ 1 . Several shapes can be used for the profile, including turbulent ones. As far as laminar profiles are concerned, we refer to [49, ch.X] for elements of comparisons between different profiles. We shall assume that ϕ depends solely on the variable ξ according to the similarity principle [49] on velocity profile over a flat plane at zeroincidence. For the sake of clarity, let us briefly present in what follows several classical profile functions in viscous layer from which we establish some instructive laws on H and f 2 for our study.
Pohlhausen polynomial profile. This kind of approach is known as the Pohlhausen solution which consists in considering a polynomial approximation of velocity profile. The polynomial coefficients are chosen such that ϕ(ξ) verifies boundary condition (2.13). The first and also the simplest case consists in a linear profile Higher order polynomial can be derived by imposing additional conditions 0 = ϕ (1) = ϕ (1) = · · · which mean that the transition between the viscous layer and the inviscid layer must be smooth. Thus, the second and third-order profiles write We notice that although these profiles are quite different, both lead to very close values of H, f 2 . Nevertheless, they do not allow the observation of separation for decelerated flows. This difficulty can be overcome only by using a fourth-order polynomial with a free parameter Λ. The resulting profile function, named Pohlhausen4 in the following, takes the form The parameter Λ is related to the pressure gradient, and therefore to the variation of velocity of the inviscid flow. Indeed, taking the second-order derivative ϕ (0) and by evaluating Prandtl momentum equation (2.4) at y = 0 we deduce where the last equality is the stationary version of the momentum equation (1.12) in the inviscid layer. Since by definition ϕ(ξ) ≤ 1, we have Λ ≤ 12, and we emphasize that for Λ ≤ −12, the velocity profiles exhibit negative regions that correspond to reverse flow (see Figure 3, left side).
We turn now to study the shape and friction factors based on the profile under consideration. First, substituting the fourth-order polynomial into definitions (2.14) and (2.15) allows to express H, f 2 as explicit functions of Λ, but they are omitted here for the sake of compactness. Nevertheless we notice that these relations are just formal since the "physical thickness" ∆(t,x) remains unknown once the velocityū in viscous layer has not yet been solved. In practice, it is more convenient to replace ∆ by the displacement thickness δ 1 . A possible way, according to [38,37], is to introduce a new parameter Λ 1 , inspired from the definition of Λ, given by Consequently, the factors H, f 2 can also be expressed as functions of Λ 1 . We represent on Figure 3, on the right side, the functions H(Λ 1 ) and f 2 (Λ 1 ) with Λ ranging from −24 to 12 that corresponds to −6 ≤ Λ 1 ≤ 0.48. Finally, we present in Tab. 1 values of H and f 2 corresponding to special cases: Λ = 12 (limit of physical range), Λ = 0 (no pressure gradient namely Blasius solution) and Λ = −12 (incipient separation). Falkner Skan profile. Polynomial profiles, despite their simplicity, are rather artificial. Their construction is based only on some suitable boundary conditions. An alternative approach, that might be more interesting, is to use exact solutions of boundary layer equations in order to establish more physical closures. Such an approach can be done by employing the solution to Falkner-Skan equation [18]. It plays an important role to illustrate the main physical features of boundary layer phenomena. This solution describes the form of an external laminar boundary layer of a flow over a wedge. The Blasius solution for a flat plate is a particular case of this solution.
Falkner-Skan equation consists of a third-order boundary value problem whose resolution is still complicate (see e.g. [9,56,31]). We do not present here any details on the resolution of Falkner-Skan equation but focus on the construction of closure formulae and compare the obtained results with those given by Pohlhausen4. First, we solve the Falkner-Skan equation on the whole physical range of pressure gradient, corresponding to the case of accelerated, decelerated and reverse flows, to obtain all the values of the triplet (Λ 1 , H, f 2 ). Next, we find out a numerical relation between these parameters by inspiring from the approach presented for Pohlhausen4. On Figure 4, it is found that the Pohlhausen4 closure, although its purely algebraic derivation, presents a good agreement with Falkner-Skan when Λ 1 ≥ 0 corresponding to accelerated flows. In particular, exact value of Blasius solution   As we can see on figure 4, this numerical law presents a good agreement near Blasius solution (both accelerated and decelerated flows). These regions are also the most concerned cases in river hydraulic application (i.e. with small bed perturbation). Especially, when plotting f 2 as function of H, an excellent agreement is found. This is why we adopt (2.18) for our numerical study of the present model.
Finally, it is important to emphasize that both polynomial and numerical closure for H and f 2 are based on steady solutions of boundary-layer equations.

Extended shallow water model
We are now in position to obtain the extended model we are looking for. Depth-integration the mass and momentum conservation equations of RNSP system yields shallow water equations (1.16) and (1.17), as presented in Sec. 1.3. Compared to that of classical shallow water model [21], we have noticed, on the one hand, the dependence in δ 1 , δ 2 of the momentum equation, and on the other hand, the rise of parietal constraints in the right-hand side at order 1 inδ. This motivates a new closure for the momentum flux and so the system has to be coupled with the von Kármán equation presented above.

Towards the extended model
This section is devoted precisely to the coupling between depth-integrated shallow water equations and the von Kármán equation. It merely emphasizes that relation (2.9) does not give any closure for the momentum flux. This is achieved by obtaining a closure on the momentum thickness δ 2 , through the study of the viscous layer, as we did above in Section 2.3. In what follows, we show that expression (2.9) for the momentum flux is in fact the most convenient, since it takes into account the effect of the viscous layer, and we clarify clarify the role of the von Kármán equation. To this end, we start from the system of depth-integrated equations (1.16), (1.17) and ideal fluid equation (1.13). We rewrite the momentum equation (1.17) with a generic form of the flux together with relation (2.11) on parietal constraints: where J is the momentum flux for which we seek a closure. We evidence now the fact that a convenient definition of J allows to recover the von Kármán equation from this system of integrated equations.
Proof. We start from the von Kármán equation and introduce (2.8) to obtain To this equation we add (3.1), the parietal term disappears, leading to Developing the time derivative and simplifying with (1.13) we obtain Finally, we use (1.16) to eliminate the time derivative, regroup terms and get so that, up to a constant which can be taken equal to zero by considering that the flux is zero when the velocity is zero, we have which together with (2.8) gives precisely (3.2). Conversely, we consider (2.8) and use successfully the mass and momentum balance equations Noting that ∂ x (δ 1 u 2 e ) = u e ∂ x (δ 1 u e ) + δ 1 u e ∂ x u e we recover as required the von Kármán equation.
In  Denoting by δ * 1 the thickness obtained using (2.8), we have ∂ t (u e (δ 1 − δ * 1 )) − u e ∂ x (u e (δ 1 − δ * 1 )) = 0. In other words the error between these displacement thickness is constant along the characteristics of the ideal fluid. Hence, if initially δ 1 = δ * 1 then it is true for all times.
Proof. From (2.8) Now using the von Kármán equation to eliminate δ 2 , we obtain which is the desired result.
In summary, propositions 3.1 and 3.2 show that there are two equivalent possibilities to compute δ 1 in order to close the system. The first one consists in adding the algebraic relation (2.8), which gives somewhat an equation of state. The second one makes use of the von Kármán equation that we can also rewrite in the following form emphasizing the fact that δ 1 u e is advected with velocity u e . This will be useful in particular for numerical purposes.
Putting together the results of Proposition 3.1 and the closure formulae (2.14) and (2.15) on the velocity profile in the viscous layer, we can now rewrite momentum equation (1.17)-or (3.1)-of the viscous shallow water model in the form Moreover, the relations between h, U , u e and δ 1 through the displacement thickness (2.8) and (2.9) directly imply the following expressions for the momentum flux The last expression clearly emphasizes that we are able to compute a non constant Boussinesq coefficient, indeed from (1.19) we have Forδ = 0 we recover the classical shallow water system with Boussinesq coefficient equal to 1. As soon as viscosity effects arise, that isδ > 0, we have not only the friction term on the right-hand side of (3.5) but also a correction of the same order one inδ to the hydrostatic pressure. Notice that in [47] a similar correction in the flux of the shallow water system is proposed to improve the study of roll-waves. In the context of a thin viscous layer (δδ 1 /h 1) we consider here, one can observe from (3.6) that the Boussinesq coefficient is very close to unity, so that its impact on the the velocity correction is negligible. Therefore we focus on the study of the friction term.

Final equivalent ideal fluid formulations
The aim of this section is to propose extended shallow water (ESW) models involving an ideal fluid with velocity u e and to couple them with the von Kármán equation describing evolution of the displacement thickness δ 1 . The first step towards these models consists in providing two systems where we keep equation (1.16) for mass conservation and constitutive relation (2.8) on the displacement thickness. With these relations, momentum equation (3.5) can easily be reformulated as It is clear on this formulation that among the three following equations: inviscid momentum (1.13), viscous momentum (3.5) and von Kármán (3.4), once two equations are satisfied then so is the third one. From this consideration, the ESW can be expressed with two equivalent formulations: the first one describes an ideal fluid living on a viscous layer, which becomes some apparent topography; the second one represents an equivalent ideal fluid over the whole depth h. It turns out that this last formulation is the most suitable for numerical studies. Here these systems are obtained by straightforward manipulations of the equations, but it is noteworthy that they can be derived as well from the Euler system with appropriate boundary conditions, as we shall see below. (3.7) One can see that the two first equations of the system represent a shallow flow, of thickness H, over a modified topography, namely f b +δδ 1 . This has to be related to the so-called "apparent topography" formulation, where for numerical purposes the friction term is rewritten as the derivative of some function, see [4,5]. Here this derivative arises in a natural way, together with additional time derivatives in the mass and momentum equations. As mentioned, the apparent topography formulation can be obtained as well from integration of the Euler system over the vertical, but from the modified topography f b +δδ 1 to the free surface η. The key point is that the boundary condition on the interface between the viscous layer and the ideal fluid is a non-penetration As for the classical shallow water model, it is straightforward to obtain from the first two equations of system (3.7) the following energy balance equation: In contrast with the classical shallow water system, the "dissipation" of energy here is driven by the dynamical behaviour in time of the displacement thickness δ 1 . If ∂ t (δδ 1 ) ≥ 0, as for the Blasius-Stokes solution presented below in Section 5.2, we have actually dissipation of energy by the bottom friction, thus some stability of the solutions. If this not the case, the problem of energy dissipation is open.
Interactive Boundary Layer formulation. This model describes an ideal fluid on the whole water depth h. Multiplying (3.4) byδ and substituting it into (3.5) leads to a model being very similar to the usual shallow water one. The resulting system reads (3.10) Once the conservative variables h, hu e and δ 1 u e are solved, the averaged velocity U , if needed, can be recovered using relation (2.8). An noticeable feature of this formulation is that the friction is no longer explicitly present in the momentum equation, it is replaced by an advection term on the "momentum" δ 1 u e . This term actually represents the momentum exchange between these two layers. System (3.10) enjoys a more conservative structure than the previous formulation, which makes it more suitable for numerical discretization. This formulation of the model can also be obtained by integrating the Euler system over the whole water depth, but the non-penetration condition (1.14) has to be replaced by a slightly modified one, called transpiration condition, that writes v e | y=f b = u e f b +δ∂ x (δ 1 u e ). (3.11) which is in fact a formulation of the "Interactive Boundary Layer" (IBL) or "Viscous Inviscid Interaction" in aerodynamics, see [36]. Comparing with (1.14), one can see that the transpiration condition is a correction of order one inδ of the non-penetration condition due to the development of viscous layer. We shall refer to (3.10) as the IBL formulation in the following. For the sake of completeness, we briefly recall how condition (3.11) is derived. We wish to estimate the vertical velocity v e of the ideal fluid at the bottom y = f b . We first notice that, since ∂ y u e = 0, integrating incompressibility equation (1.6) of ideal fluid over the whole water depth yields Putting things together leads to the required transpiration boundary condition (3.11).

Numerical analysis of the IBL formulation
We propose here a numerical implementation of the IBL formulation (3.10). This model was chosen because of its relative simplicity compared to the apparent topography and to the four equations models. Also, it can be related to the previous implementation in the context of rigid pipes, see [38]. In order to design a finite volume solver, we first rewrite this system in vector form as where the conservative variable W , the flux F (W ), the convective term B(W ) and the source term τ (W ) are defined by The system is numerically solved by a splitting method. First we solve the so-called convective part together with the topography source term − hf b F r 2 , in a well-balanced way. Next, we solve the friction part by a semi-implicit method.

Eigenvalue analysis
For numerical purpose and stability analysis, we first study the hyperbolicity of the convective part (4.1) of the model. To this end, we rewrite (4.1) in quasi-linear form ∂ t W + A(W )∂ x W = 0. The system is said to be hyperbolic if the convective matrix A(W ) is R-diagonalizable, and strictly hyperbolic if the eigenvalues are distinct. Estimating these eigenvalues is also important to design an explicit finite volume scheme for the system. As eigenvalues are invariant by changing variables [23], it is therefore more convenient to make the variable change W → Y (W ) := (h, u e , δ 1 u e ) t and study eigenvalues of the corresponding convective matrix where we have denoted the partial derivatives Their explicit expressions can be computed as well once a closure formula for the shape factor is provided, since the shape factor H is a function of Λ 1 , so H depends only on u e and δ 1 u e .
The characteristic polynomial ofÃ(Y ) reads The polynomial P SW has always 3 roots, denoted λ 0 1,2,3 , which represent the wave speeds of the system with no coupling between shallow water equations and the Von Kármán one. Indeed, the first two roots λ 0 1,2 express the propagation velocities of ideal fluid while the last one λ 0 3 approximate that of the viscous layer: From this it follows that, forδ small enough, the characteristic polynomial admits 3 real eigenvalues and the system turns out to be hyperbolic. More precisely, denoting λ − < λ + the roots of P SW (λ), the system is hyperbolic when dδ lies between P SW (λ ± ), see figure 5, that is It is straightforward to verify that for the case of closure (2.18) based on Falkner-Skan solutions, the third wave speed reads 4) and in particular when Λ 1 = 0, i.e. the Blasius solution, the viscous layer propagates downstream at velocity λ 0 3 = u e /H 0.39u e (see figure 6-left). Even within the hyperbolic regime, the eigenstructure is given only implicitly, due to the form of nonlinear coupling between the ideal fluid and viscous layer. As a result, when dealing with numerical methods requiring characteristic field decomposition, e.g. Roe type method, the eigenvalues need to be computed by numerical root finding. In the absence of analytic expressions for the eigenvectors, building desirable properties for such a scheme may be more difficult. We propose hereafter a HLL type scheme [27] taking advantage that the numerical flux arises directly from the governing equations (4.1) and only an estimation of lowest and the largest wave speeds λ L,R is required. Such a wave speeds estimation can be done by using accurate Nickalls's bounds [43], which writes in which we have used (4.4) provided from the Falkner-Skan closure. On figure 6 we compare the shallow water wave speeds λ 0 1,2 , the viscous layer one λ 0 3 and the bounds λ L,R given by (4.5). Rescaling these velocities by u e , the result can be considered as functions of the local Froude number F r 0 := F r u 2 e /h and the pressure parameter Λ 1 . The left figure shows the dependence on F r 0 in the case of Blasius solution, i.e. Λ 1 = 0, both for subcritical and supercritial regimes. We find that the estimation λ R for the largest wave speed is very accurate. On the right figure, we display the dependence on Λ 1 for a fixed value of F r 0 , e.g. by considering the case of critical flow so λ 0 1 = 0. As we can see, the velocity λ 0 3 can be negative-the wave associated to viscous layer propagate upstream-for large reverse flow.

A Godunov-type finite volume scheme
Let us recall some basic notations of finite volume discretization. We introduce a space step ∆x and a time step ∆t, both assumed to be constant for simplicity. The computational domain is discretized by a sequence of points x j+1/2 := j∆x for j ∈ Z. We define W 0 j a piecewise constant approximation of initial condition on each control volume C j :=]x j−1/2 , x j+1/2 [. The time step ∆t for a mesh size ∆x has to atisfiy the well-known CFL condition ∆t ≤ ∆x 2|λ max | , (4.6) where λ max is the largest eigenvalue expressing the fastest wave speed of the system. This condition ensures that information of each Riemann problem at a cell's interface does not cross more than one cell.
Convection step. Assume that the solution W n j at time t n is known. Godunov-type schemes compute the solution to (4.1) at the next time level t n+1 := t n + ∆t by building first an approximate solution W ∆ (x, t) of the Riemann problem at each interface x j+1/2 with initial data {W n j } j∈Z , and next averaging W ∆ (x, ∆t) on each control volume to obtain a piecewise constant solution W n+1/2 j . Given initial data (W L , W R ) of local Riemann problem, we adopt a simple Riemann solver W ∆ (x, t) composed by three discontinuity waves propagating with velocities λ L ≤ 0, λ 0 = 0, λ R ≥ 0 and two intermediate states W * L and W * R in the star region, see figure 7. In the following λ L,R are given by formula (4.5), and the zero velocity λ 0 corresponds to the stationary contact discontinuity associated with the topography. The third eigenvalue λ 3 was not used here because the analytical expression of the Riemann invariant is unknown.
Under CFL condition (4.6), the first-order three-points finite volume scheme writes where the left-and right-numerical fluxes F L,R j+1/2 := F L,R (W n j , W n j+1 ) are given by (4.8) Therefore, designing such a scheme consists in determinating the intermediate states W * L,R in the star region. Figure 7 -A three-waves approximate Riemann problem.
According to [27] the approximate solver W ∆ (x, t) must be consistent with the exact solution W R (x, t) in the sense that Applying to the conservation law (4.1), this integral consistency condition provides that the intermediate states satisfy the following relations (4.10) It is well-known that the non-conservative products arising in the source terms may not make sense as distributions. It is possible to give a rigorous definition using Vol'pert's calculus on BV functions [55]. This choice leads to following approximations 12) which preserve at least the lake-at-rest equilibrium state, that is u = 0, Consistency conditions (4.9)-(4.11) have to be complemented by three additional relations in order to solve completly the intermediate states. The Riemann invariants associated to the stationary contact wave can be used to provide these missing relations (u e ) * 2 Note that the first two equations are the analytical Riemann invariants of the system applied for intermediate states while the third one related to δ 1 u e is just an approximation. In fact, the last analytical Riemann invariant is not known explicitly due to the form of nonlinear coupling between the ideal fluid and viscous layer. From the numerical point of view, this consists in approximating (δ 1 u e ) * L,R in the star region by only one averaged value (δ 1 u e ) * , as in the classical HLL scheme [27]. Moreover, this choice leads to (hu e ) * L = (hu e ) * R := q * and, together with (4.14), we find again the well-known Bernoulli relations of shallow water equations Plugging (4.13), (4.15) into (4.10) and (4.11) allows us to solve the discharges q * and (δ 1 u e ) * in the star region. Next, intermediate water depths h * L,R are obtained from (4.9) and (4.16). Similarly to the case of classical shallow water model, it can be shown that the present scheme is accurate and well-balanced. We refer to [24] for more technical details and discussions.
Friction step. Once the convection step (4.7) is done, the next step is to take into account the friction source term of von Kármán equation to modify the displacement thickness. This consists in solving equation (4.2) with initial data W n+1/2 to get the solution W n+1 at the next time step. To this end, we use a simple semi-implicit scheme which writes This discretization leads to a second-order equation for δ n+1

1
. Under condition (δ n+1/2 1 ) 2 + 4(f 2 H) n ∆t ≥ 0, which implies an additional restriction on ∆t only in the case of reverse flow In practice when using the Falkner-Skan closure, condition (4.17) is not restrictive, the time step is rather controlled by the CFL condition (4.6) of the convection step.

Multi Layer formulation
We conclude this section by a short presentation of the so-called multi-layer Saint-Venant model proposed in [2], which we shall use as a reference for comparison in the next section. The authors consider a superposition of shallow water systems each one interacting with its neighbours: where for each layer α = 1, . . . , N , h α = α h denotes the layer thickness, α > 0 being a given constant, N α=1 α = 1, and u α (t, x) the averaged velocity in layer α. The source terms m α+1/2 , τ α+1/2 stand for the momentum exchange and the friction between layers α and α + 1 respectively.
As proposed in [2], we solve this multilayer shallow water (MLSW) system (4.19) using a first-order finite volume scheme in which the numerical flux is built by a kinetic formulation. The friction term between the layers τ α+1/2 is discretized in an implicit way by We notice that the third expression is due to the no-stress condition at the surface while the second one, expressing the bottom friction, is based on the no-slip condition and on a first-order expansion of the velocity. Using this model imposes some constraints on the vertical discretization: it requires a very thin layers near the bottom in order to accurately compute the friction while the velocity of two adjacent layers (in the viscous region) must not be too different in order to preserve the hyperbolicity of the model. This later condition implies that relative thickness α of the layers can be varied but only gradually. For the present study, a discretization such as α = z α − z α−1 with z α := e 10α/N − 1 e 10 − 1 for 0 ≤ α ≤ N := 100 seems to give satisfactory results. Nevertheless the simulation with MLSW model is computationally expensive.
We interpret this MLSW system as a numerical scheme to solve the RNSP equations (2.4)-(2.6), the layers being the numerical discretization along the vertical direction (see [22] for a similar point of view in elastic tubes). In our approach, at the RNSP level, the two superposed layers have different physical properties (viscous/ ideal fluid), and analyzed through asymptotic rescalings of the equations. At the integrated level, the closest formulation is the apparent topography, we have again two superposed layers of different physical nature, and their relative thickness,δδ 1 /h, is not fixed but evolves in time, in contrast with the multi-layer model.

Numerical illustrations
The aim of this last part is to give a few illustrations of the behaviour of the ESW system (3.10). We are aware that more accurate analysis is mandatory, in both the numerical approach (in particular numerical boundary conditions) and the qualitative and quantitative behaviour of the model. Fisrt we give a convergence study based on a steady-state solution. A second step is devoted to a comparison between the ESW solutions and those of the classical viscous shallow water system. Finally we compute the solutions over a small bump, in order to evidence the above mentioned phase-lag of the friction term, and the behaviour of the model with respect to various parameters.
In all test cases, we considered a very thin viscous layer by settingδ = 10 −3 . The Froude number F r is set to unity meaning that the longitudinal velocity was scaled by the reference celerity √ gh 0 . The computation domain was [0, L]. Initial conditions was h(0, x) = h 0 , u e (0, x) = 1, δ 1 (0, x) = 0 ∀x ∈ [0, L]. We set h 0 = 2 for subcritical test cases, so F r 0 0.7, while we used h 0 = 0.5, i.e. F r 0 1.4, for supercritical cases. On the boundary conditions, we considered at x = L a free outflow boundary. At x = 0 for both sub-and supercritical flows we imposed a constant velocity u e = 1 with flat profile, so that δ 1 = 0. For supercritical flows, a constant water depth h = h 0 was imposed as well. In the case of subcritial flows, h was computed using the Riemann invariant of classical shallow water model. This approach is only an approximation for ESW, more in-depth study on numerical boundary condition is of course needed.

Convergence study on a Blasius-like steady solution
In this section we investigate an approximate stationary solution of the ESW model which is an order one perturbation of a basic stationary solution to the ESW system, namely (h 0 , u 0 e , δ 0 1 ), where h 0 and u 0 e are constant solutions of the frictionless shallow water system, and δ 0 1 is the classical Blasius profile for the von Kármán equation on a flat plate. We look for a solution to the ESW system at first order inδ h = h 0 +δh 1 , u e = u 0 e +δu 1 e , δ 1 = δ 0 1 +δδ 1 1 .
A very interesting feature of this solution is that the order one terms are not necessarily stationary. This gives an explicit illustration of the actual interaction between the viscous layer and the perfect fluid. Plugging the expansion (5.1) in (3.10), we recover through the first two equations a standard inviscid shallow water model, for which a basic stationary solution consists in constant h 0 and u 0 e . Now we turn to the stationary von Kármán equation In this equation, H and f 2 depend on Λ 1 = δ 2 1 ∂ x u e =δδ 0 1 ∂ x u 1 e + O(δ 2 ). Therefore at zeroth order H and f 2 are constant, so that we indeed recover the classical Blasius solution We use the solution (h 0 , u 0 e , δ 0 1 ) as the basic solution in the expansion (5.1), and turn now to order one terms. Once again, straightforward computations lead to uncoupling the first two equations, yealding the following linearized shallow water model: Notice that this system has a stationary solution, given by  Figure 8-right we display the results of a mesh convergence study on the gap between δ 1 and its zeroth order approximation δ 0 1 . The convergence study was performed in both sub-and supercritical regimes. First we remark that when the mesh size ∆x is small enough, namely ∆x ≤ 10 −4 , numerical results reach the model error. Indeed, we obtained at this spatial resolution that Next, we notice that the supercritical case converges faster than the subcritical one. This could be explained by the fact that, on the one hand, numerical treatment of the left boundary condition is more accurate in the supercritical case as we have noticed before; on the other hand, it is well known that the HLL-type numerical flux (4.8) is also more accurate in that case.

Impulsively started flow over a flat bed
We turn now to a configuration introduced by Stewartson [51, 50, Sec. 3] as a simple test-case to study unsteady boundary layer solutions. It consists in a semi-infinite flow impulsively started from rest at t = 0 with constant velocity u e , see figure 9. The fluid is injected continuously at x = 0 with a constant velocity, the flow must satisfy the no slip boundary condition for x > 0. The solution exhibits different behaviours depending on two asymptotic regimes: for large t (or small x) we recover the Blasius solution (5.2); conversely, for large x (or small t), convective terms in the momentum equation (2.4) are negligible so the Prandtl system reduces to Stokes' first problem (also called Rayleigh problem by Stewartson). The solution for the velocity profile can be expressed using the erf error function: Stewartson noticed that an approximate integral form can be used to solve this problem. Assuming a fixed profile, e.g. Blasius (constant) value of H and f 2 , von Kármán equation (2.16) can be rewritten in the form and it has to be complemented with the following initial and boundary value conditions: The solution is readily obtained by the method of characteristics: 2f 2 Ht otherwise, Transition zone is found at the characteristic line x = uet H , i.e. viscous layer informations propagate at velocity u e /H as we have seen in equation (4.4). It is worth noticing that this solution with constant u e creates an unbounded viscous layer, without regards of the limitation of the water depth. Hence it is clearly physically not valid for large x or t.
On figure 10-left, we plot on the left side the evolution of the displacement δ 1 along x at different times. The two-regimes behaviour of the unsteady solution (5.5) is qualitatively well recovered: δ 1 is constant in x far from the entrance and increases in time; it reaches the Blasius solution which is steady. But, quantitatively, we observe a gap between the solution of ESW and the Stokes one, which increases in time. This is due to the choice of a fixed Blasius profile in the viscous layer. Regarding equations (5.2) and (5.4), the Blasius and Stokes profile have not exactly the same shape and this leads to different values of H and f 2 . More precisely, in the Stokes region the model predicts δ 1 = 1.067 √ t while the exact solution is δ 1 = 1.128 √ t. To illustrate the transition in profiles from Stokes to Blasius, we propose here to compare the numerical results given by the ESW model with those of the multi-layer scheme described in Section 4.3 (1.7)-(1.11). On figure 10-right, we plot the difference between the numerical solutions and the two asymptotic solutions as function of t/x. We plot on this figureτ b √ x which is x/πt for small t/x (Stokes solution) and becomes 0.322 (the Blasius value) for large t/x. On the solution of ESW, the unsteady part of ESW solution presents a small difference as well. This is simply due to the fact that the Blasius profile and the Stokes one have not exactly the shame shape as we have observed for δ 1 . The value f 2 H given by the Stokes solution is greater than that of the Blasius solution. Good agreement was found on the solution of MLSW: the solution tends from Stokes to Blasius behaviours; the transition between these two regimes is clearly captured.

Flows over a small bump
This final section considers cases with a flat bottom with a bump. The lenght of the bump is such that those compatible with classical shallow water model. Hence, we consider now a domain x ∈ [0, 2] with a small bump of Gaussian form at the center: with α and σ given. Falkner-Skan closure, which includes as well the Blasius one, is used to compute the shape and friction factors. As defined by (2.18), we have to compute the pressure gradient parameter Λ 1 = δ 2 1 ∂ x u e for which accurate approximation for the partial derivative ∂ x u e is required in order to provide a faithful representation of the shape of thin viscous layer. Therefore, we test second order and fourth-order finite difference derivative: Influence of the flow regime. We use first α = 0.01 and σ = 0.1. As in the test case of an impulsively started flow, the displacement thickness develops from zero, at t = 0. It reaches a steady value when the time is large enough, typically when t > xH/u e 6 from characteristic solution (5.5) applied for this case. However, this steady solution is no longer Blasius but is slightly perturbed due to the presence of the bump; this perturbation depends furthermore on the flow regime, as one can observe on figure 11.
For sub-critical case, the flow is accelerated on the upstream side of the bump while it is decelerated on the downstream side. As a consequence, the displacement thickness decreases before the crest and increases after. This behaviour is also well reported by classical shallow water model, even with linearized solution (1.24). However, ESW provides an asymmetric friction due to inertia effect of the fluid, compared with (1.24). More precisely, the friction reaches its maximum before the crest for sub-critical flow while for super-critical case, it becomes maximum after the crest, see figure 11 (right). This is in fact the important phage-lag behaviours that we have found to cope with ESW. Phase-lag reduced for shorter bump. Here we investigate the influence of the length of the bump on the friction computed with ESW model. Returning to sub-critical case, we perform now the same test case but with a shorter bump, by imposing σ = 0.05. Regarding the result given with Falkner-Skan closure, i.e. when the shape factor depends on velocity gradient ∂ x u e , figure 12 (left) shows that the variation of friction is more important than the case with σ = 0.1. The phage-lag is always observed. To highlight the role of the closure imposed in viscous layer, we do again this test but with a constant value of shape factor (e.g. by using Blasius closure). It can be observed on figure 12 (right) that both amplitude and phase-lag of friction are significantly reduced compared to the case with (variable) Falkner-Skan closure. Moreover, constant shape factor does not allow the friction to decrease enough in decelerated region localized at downstream side of the bump. As a consequence, this kind of closure is unable to recover reverse flow whatever the bump shape.
Larger friction for higher bump, reverse flow observation. This last test case is devoted to highlight the possibility of ESW to capture indeed reverse flow if the bump is high enough. We fix the bump length to be σ = 0.1 and change α. Falkner-Skan closure (2.18) is of course imposed. Figure 13 shows the displacement thickness and the friction for increasing α up to negative value of friction (which defines steady boundary layer separation). Associated to boundary layer separation is an increased displacement thickness.
Note that accurate estimation of velocity gradient ∂ x u e is very important to capture this particular phenomena. That is why we have used the fourth-order formula (5.6). Using a second-order approximation for ∂ x u e results just in the incipient separation. Finally, performing this test case with MLSW model allows to find as well the phase-lag, but the amplitude of friction is smaller as observed on figure 14. It seems that numerical solver of MLSW model together with the first-order approximation (4.20) onτ b , and even with 100 layers, is not enough to accurately capture the friction.
Nevertheless, plotting the friction factor f 2 as function of the shape factor H and comparing it with Falkner-Skan closure (2.18) shows that indeed, the reduced shear and shape factors computed from the profilers of MLSW are close to Falkner Skan curve. The smaller the angle, the closer the curves. The agreement is better during the accelerated phase of the flow.

Conclusion
In order to improve "shallow water" models, this paper proposed a novel description of parietal friction for free surface shallow flows (long wave approximation) in large Reynolds number limit. The proposed model relies on a perfect fluid -viscous layer decomposition. It consists in a system of four equations: conservation of mass (1.16), momentum (3.5), an ideal fluid equation (1.13), and the von Kármán equation (3.4). The first two are similar to classical shallow water system with a slight correction on momentum flux, and with a specific friction term.
Two equivalent versions of the model can be obtained by interpreting them under two points of view. In the first one, the viscous layer acts as a new topography for the model, see (3.7). In the second one, this is an example of "Interactive Boundary Layer" or "Viscous Inviscid Interaction", see (3.10), and it seems to be the most convenient choice for numerical purposes. In these models the friction term is no longer an empirical combination of velocity and depth (as in usual laws such as Darcy or Manning) but the result of a viscous layer like approach. A crucial point at this stage is the choice of an apropriate closure for the shape factor and the friction factor in the viscous layer, in order to obtain a closed system after integration.
As it is, the friction term actually depends on the topography, as evidenced by the examples provided in the last part of this work. In particular its maximum is reached before or after the summit of a bump, depending on the criticity of the flow. Also, possible boundary layer separation with recirculation downstream (in subcritical case) of the bump can be observed as in similar cases in litterature ( [37,38,39] in case of flows in pipes, and as in preliminary comparisons with multilayer shallow water [2]). Those two behaviors are impossible to observe in the classical Shallow Water model.
Our proposed approach is restricted to laminar flows, but the ideas developed here can be extended with little modifications to mean turbulent profiles.
The main drawback of this model is the viscous layer/ideal fluid decomposition, which forbids the viscous layer to fill all the water depth far downstream of the bump. Extra modelling in this direction, as well as a more careful study of the closure laws in the viscous layer have now to be worked out and tested. Also, the system is conditionnally hyperbolic, and the numerical study has to be improved, in particular boundary conditions. In spite of these limitations, the model is a good compromise between the classical shallow water system and RNSP equations (discretized with the Multilayer Saint-Venant scheme [2]), much more costly in computation time.