G L Richard

M Gisclon

C Ruyer-Quil

J P Vila

Optimization of consistent two-equation models for thin film flows

Keywords: free-surface flows, thin films, capillary flows, solitary waves, interfacial instability

A general study of consistent two-equation models for thin film flows is presented. In all models derived by the energy integral method or by an equivalent method, the energy of the system, apart from the kinetic energy of the mean flow, depends on the mean velocity. We show that in this case the model does not satisfy the principle of Galilean invariance. All consistent models derived by the momentum integral method are Galilean invariant but they admit an energy equation and a capillary energy only if the Galilean-invariant part of the first-order momentum flux does not depend on the mean velocity. We show that, both for theoretical and numerical reasons, two-equations models should be derived by a momentum integral method admitting an energy equation leading to the structure of the equations of fluids endowed with internal capillarity. Among all models fulfilling these conditions, those having the best properties are selected. The nonlinear properties are tested from the speed of solitary waves at the high Reynolds number limit while the linear properties are studied from the neutral stability curves and from the celerity of the kinematic waves along these curves. The latter criterion gives the best consistent way to write the second-order diffusive terms of the model. Optimized consistent two-equation models are then proposed and numerical results are compared to numerical and experimental results of the literature.

Introduction

The practical importance of thin film flows of Newtonian liquids led to many attempts to derive a suitable model for the evolution of the film thickness h and often also for the evolution of the mean velocity or of the discharge, both accurate and easy to solve numerically. A direct numerical simulation of the Navier-Stokes equations for these flows in real situations having a huge numerical cost, models of reduced dimensionality must be used. Many models are derived in the long-wave assumption with an asymptotic method. Consistency in the implementation of this method is required for the model to predict in particular the correct instability threshold of the equilibrium flow with a uniform depth and a parabolic velocity profile.

The definition of the consistency of a model is relative to a set of hypotheses concerning the order of magnitude of the dimensionless numbers of the problem which are, in the case of thin film flows, the Reynolds number, the Froude number, the Weber number and the slope of the bottom wall. In the present case, since the film is thin, it is possible to define a small parameter ε which is the ratio of the typical fluid depth h N to the wavelength L:

ε = h N L 1 (1)
As an asymptotic method is meaningless if there is more than one small parameter, the order of magnitude of all dimensionless numbers must be defined in relation to ε. Then all flow fields, i.e. the pressure and the components of the velocity and of the viscous stress tensor, are expanded as formal power series of the small parameter ε. These fields are supposed to be the exact solution of the Navier-Stokes equations. There is only one asymptotic expansion of each of these quantities depending only on the fluid depth h and on its derivatives. These expansions are found by using the Navier-Stokes equations, the boundary conditions and the constitutive relation of Newtonian fluids (see Richard et al. [START_REF] Richard | A three-equation model for thin films down an inclined plane[END_REF] for more details).

These expansions are then inserted into the equations of the model. If the leading terms in the model are of O [START_REF] Richard | A three-equation model for thin films down an inclined plane[END_REF] and if all terms vanish except for a remainder of O(ε n+1 ) then the model is said to be consistent and n is the order of consistency of the model. Moreover a quantity is said consistent to another one if the relative difference between them is of O(ε) when using these asymptotic expansions.

The continuity equation averaged over the depth can be used to obtain models with only one equation for h. The first one was the Benney equation which is a consistent one-equation model (Benney [START_REF] Benney | Long waves on liquid films[END_REF]) but its solutions can blow up if the equilibrium flow becomes unstable (Pumir et al. [START_REF] Pumir | On solitary waves running down an inclined plane[END_REF]). Attempts to solve this problem led to inaccurate models (Ooshida [4]). It is therefore widely acknowledged that at least two variables are needed to describe accurately the dynamics of waves that appear in films on an inclined plane.

A model with four equations, fully consistent to the second order and giving very accurate results, was proposed by Ruyer-Quil and Manneville [START_REF] Ruyer-Quil | Improved modeling of flows down inclined planes[END_REF] but its complex mathematical structure makes its numerical resolution difficult. For a greater resolution easiness, accurate models with two equations are actively sought for. The two variables are usually h and the average velocity U (or the flow discharge q = hU ). Shkadov [START_REF] Shkadov | Wave flow regimes of a thin layer of viscous fluid subject to gravity[END_REF] proposed the first twoequation model with the assumption that the velocity profile remains parabolic even out of equilibrium. Although it can give accurate results in some situations, the Shkadov model is not consistent because the deviation from the parabolic profile cannot be neglected. These deviations were measured by experiments (Alekseenko et al. [START_REF] Alekseenko | Wave formation on a vertical falling liquid film[END_REF]) and calculated by numerical simulations (Bach and Villadsen [START_REF] Bach | Simulation of the vertical flow of a thin, wavy film using a finite-element method[END_REF], Malamataris et al. [START_REF] Malamataris | Solitary waves on inclined films: Flow structure and binary interactions[END_REF], Malamataris and Balakotaiah [START_REF] Malamataris | Flow structure underneath the large amplitude waves of a vertically falling film[END_REF]). Consistent two-equation models were then derived by Ruyer-Quil and Manneville [START_REF] Ruyer-Quil | Improved modeling of flows down inclined planes[END_REF], Usha and Uma [START_REF] Usha | Modeling of stationary waves on a thin viscous film down an inclined plane at high Reynolds numbers and moderate Weber numbers using energy integral method[END_REF], Mudunuri and Balakotaiah [START_REF] Mudunuri | Solitary waves on thin falling films in the very low forcing frequency limit[END_REF], Abderrahmane and Vatistas [START_REF] Abderrahmane | Improved two-equation model for thin layer fluid flowing down an inclined plane problem[END_REF], Luchini and Charru [START_REF] Luchini | Consistent section-averaged equations of quasi-one-dimensional laminar flow[END_REF]. All these models come down to the same family based on the depth-averaged mass conservation equation and work-energy theorem (or theorem of the kinetic energy) by the energy integral method (EIM). They differ from each other in the second or thirdorder terms but are identical at the first order of approximation. Another family of models is based on the depthaveraged mass conservation equation and momentum balance equation by the momentum integral method (MIM). Because it is harder to derive a consistent model by this method than by the EIM (see Luchini and Charru [START_REF] Luchini | Consistent section-averaged equations of quasi-one-dimensional laminar flow[END_REF]), there are fewer models in this family. The models derived in Noble and Vila [START_REF] Noble | Thin power-law film flow down an inclined plane: consistent shallow-water models and stability under large-scale perturbations[END_REF] and in Lavalle et al. [START_REF] Lavalle | A numerical reduced model for thin liquid films sheared by a gas flow[END_REF] are obtained by the MIM. They are consistent but differ in first-order terms.

Recently, a three-equation model was derived by Richard et al. [START_REF] Richard | A three-equation model for thin films down an inclined plane[END_REF] from the mass, momentum and kinetic energy equations. This approach combines the EIM and the MIM with the help of a third variable, called enstrophy, equal to the variance of the fluid velocity divided by h 2 . This raises the question of the compatibility between the EIM and the MIM for two-equation models. One goal of this paper is to determine whether it is possible to find consistent two-equation models that satisfy both the averaged momentum and kinetic energy equations, thus reconciling both the EIM and the MIM, or whether there is an in-herent theoretical flaw in the two-equation approach that misses some physical features.

A second goal of this paper is to find what consistent two-equation model gives the best results. All consistent models give similar results in the validity domain of the asymptotic expansion. However, in practice the models are commonly used beyond this domain and this means that the models performance must be tested in conditions where the derivation assumptions are not so well ensured. In particular, the linear properties of the model can be tested on the calculation of the neutral stability curves and the nonlinear properties can be evaluated on the speed of the solitary waves at the limit of the high Reynolds numbers. Note that the performance of two-equation models can be acceptable to predict the speed of solitary waves especially at high Reynolds numbers but is significantly poorer for the prediction of the amplitude (Chakraborty et al. [START_REF] Chakraborty | Extreme solitary waves on falling liquid films[END_REF]). The mathematical structure of the equations is an important aspect to consider in order to solve numerically the equations efficiently. The dispersive capillary terms are by far the most complex terms to handle. If some conditions are fulfilled they can be treated with the entropy stable scheme of Noble and Vila [START_REF] Noble | Stability Theory for Difference Approximations of Euler-Korteweg Equations and Application to Thin Film Flows[END_REF], later extended by Bresch et al. [START_REF] Bresch | A generalization of the quantum Bohm identity: Hyperbolic CFL condition for Euler-Korteweg equations[END_REF]. This scheme was successfully used by Lavalle et al. [START_REF] Lavalle | A numerical reduced model for thin liquid films sheared by a gas flow[END_REF] and by Richard et al. [START_REF] Richard | A three-equation model for thin films down an inclined plane[END_REF]. It requires the existence of a capillary energy and the mathematical structure of the equations of compressible fluids with internal capillarity. This leads to an extended model which adds an extra equation relative to the gradient of the fluid depth.

The physical meaning of the momentum and energy equations is discussed in §2 as well as the general form of MIM and EIM models. The problem of the Galilean invariance of the equations is discussed in §3 and the compatibility between the momentum and energy equations is studied in §4. The existence of a capillary energy and the possibility to derive an extended model is discussed in §5. The nonlinear properties are tested in §6 on the speed of the solitary waves. This allows us to select the optimal model as regards its hyperbolic part. The linear properties are tested in §7. The kinematic waves velocity along the neutral stability curve is used to find an optimized way to model the diffusive terms. The effect on the neutral stability curves is presented. Finally the equations are numerically solved in §8 in the conditions of the experiments of Liu and Gollub [START_REF] Liu | Solitary wave dynamics of film flows[END_REF] and of Dietze et al. [START_REF] Dietze | Experimental study of flow separation in laminar falling liquid films[END_REF] and the results are presented and compared to the experimental results and to direct numerical simulations.

Energy and momentum integral methods

Physical meaning of the momentum and energy equations

We study the two-dimensional gravity-driven flow of a Newtonian liquid thin film on a plane inclined at an angle θ with respect to the horizontal. In the Ox-direction parallel to the plane, the component of the fluid velocity is u. In the Oz-direction normal to the plane and directed upward, the velocity component is denoted by w. The fluid depth is denoted by h. The Figure 1 presents the main notations in the case of a typical wavy film with a main hump and several smaller capillary ripples.

The assumption of a thin film allows to define a small parameter ε as the ratio between the typical fluid depth and the wavelength L (or equivalently the typical distance of variations of the variables in the Ox-direction). This long-wave assumption is classically used to derive several shallow-water models, such as the Saint-Venant equations, by a depth-averaging procedure. A justification of the applicability of these models to describe film flows with large depth gradients in the region of wave fronts was recently proposed in Ostapenko [START_REF] Ostapenko | Justification of shallow-water theory[END_REF] with the concept of hydrostatic approximation which generalizes the long-wave approximation.

Denoting by ρ, µ, ν and p, the fluid density, dynamic viscosity, kinematic viscosity and pressure respectively, the constitutive relation of Newtonian fluids is σ = -pI + τ with τ = 2µD where σ, τ and D are respectively the Cauchy stress tensor, the viscous stress tensor and the strain rate tensor and where I is the identity tensor.

The physical governing equations are the mass conservation equation (equation of continuity) and the Navier-Stokes momentum balance equations. The work-energy theorem does not give an independent equation but can be deduced from the Navier-Stokes equations. The fluid is supposed to be incompressible. Denoting by V the velocity field and V = V , the continuity equation is

div V = 0 (2)
The momentum balance equation is

∂ρV ∂t + div (ρV ⊗ V ) = div σ + ρg (3)
and the work-energy theorem is ∂ ∂t

1 2 ρV 2 + div 1 2 ρV 2 V = div (σ • V ) -σ : D + ρg • V (4)
The symbol ⊗ represents the tensor product and the colon denotes the double contraction. The right-hand side of (3) represents the external forces acting on the fluid particle while the right-hand side of (4) represents the power of both the external and internal forces. The latter is expressed by the term -σ : D.

To build a model of reduced dimensionality, these equations can be averaged over the depth taking into account the boundary conditions which are the no-penetration and no-slip conditions on the bottom

w (z = 0) = 0 , u (z = 0) = 0 ( 5
)
the kinematic boundary conditions at the free surface

∂h ∂t + u (z = h) ∂h ∂x = w (z = h) (6)
and the dynamic boundary condition at the free surface

(σ • n) (z = h) = γ R n ( 7
)
where γ is the surface tension, R the radius of curvature of the free surface and n the unit normal vector at the free surface pointing upward. The atmospheric pressure is supposed to be constant and can be taken equal to zero. For any variable A, its averaged value over the depth is defined as

A = 1 h h 0 A dz (8)
The averaged value of u is denoted by U . Averaging the continuity equation [START_REF] Benney | Long waves on liquid films[END_REF] gives the averaged mass equation

∂h ∂t + ∂hU ∂x = 0 (9)
The momentum and kinetic equation equations are first put into dimensionless form and then averaged over the depth. The characteristic depth and velocity are the depth h N and average velocity u N of the equilibrium flow, also called the Nusselt flow, which has the parabolic velocity profile

u N (z) = g sin θ ν z h N - z 2 (10)
The expression of its average value over the depth is

u N = gh 2 N sin θ 3ν (11)
The small parameter ε is now defined as

ε = h N L (12)
The Reynolds number Re, Froude number Fr and Weber number W e are defined as Re = h N u N /ν, Fr = u N / √ gh N and W e = ρh N u N 2 /γ. We define also the dimensionless number

λ = Re sin θ Fr 2 (13)
It follows from the previous relations that λ = 3. It is convenient to use the following number (see Richard et al. [START_REF] Richard | A three-equation model for thin films down an inclined plane[END_REF])

κ = ε 2 Fr 2 W e (14)
which is supposed to be of O(

ũ = u u N w = w ε u N x = x L z = z h N t = t L u N h = h h N Ũ = U u N p = p ρgh N τxx = Lτ xx u N µ τzz = h N τ zz ε u N µ τxz = h N τ xz u N µ (15)
The equations and the boundary conditions in dimensionless form are given in Appendix A.

In dimensionless form, the averaged mass equation has the same form as [START_REF] Malamataris | Solitary waves on inclined films: Flow structure and binary interactions[END_REF]. Averaging the momentum balance equation gives (see Richard et al. [START_REF] Richard | A three-equation model for thin films down an inclined plane[END_REF])

∂ h Ũ ∂ t + ∂ ∂ x h ũ2 + h2 cos θ 2Fr 2 = λ h -τxz (0) εRe + κ Fr 2 h ∂ 3h ∂ x3 + O ε Re (16)
The third derivative is the capillary term and the terms of O(ε/Re) are diffusive terms. All these terms could be written in conservative form. Apart from these terms, the right-hand side of equation ( 16) represents the external forces acting on the averaged system i.e. the Oxcomponent of the weight and the friction on the bottom. Shearing forces internal to the flow disappeared by the averaging operation. This was of course expected since internal forces cancel in accordance with Newton's Third Law.

Averaging the work-energy theorem over the depth gives (see Richard et al. [START_REF] Richard | A three-equation model for thin films down an inclined plane[END_REF])

∂ ∂ t h ũ2 2 + h2 cos θ 2Fr 2 + ∂ ∂ x h ũ3 2 + h2 Ũ cos θ Fr 2 = 1 εRe λ h - 3 Ũ h Ũ + κ Fr 2 h Ũ ∂ 3h ∂ x3 + O (ε) (17)
The right-hand side of this equation, apart from the capillary term, corresponds to the power of both the external and internal forces acting on the system, in accordance with the work-energy theorem. More precisely, the two terms are the power of the Ox-component of the weight and the dissipative power of the viscous forces internal to the flow. Note that the Oz-component of the weight is included in the left-hand side of this equation as an average potential energy. Because of the no-slip boundary condition, the friction force on the bottom does not do work. Consequently, the physical meaning of the averaged momentum equation and of the averaged kinetic energy equation are not equivalent. In the former appears the friction on the bottom but not the internal viscous forces while it is the reverse in the latter. It follows that, even if the workenergy theorem is not an independent equation with respect to the Navier-Stokes equations, after averaging over the depth, the averaged kinetic equation is an equation independent of the averaged momentum equation. The averaged kinetic energy equation (later called simply "energy equation") is redundant with the averaged mass and momentum equations only for the equilibrium (Nusselt) flow since in this case the bottom friction is equal to the Ox-component of the weight and the viscous dissipation is equal to the power of this component. Consequently, at the equilibrium, τxz (0) Ũ = 3 Ũ 2 / h. However, this relation is not satisfied out of equilibrium. For models of reduced dimensionality derived with an averaging procedure over the fluid depth (or any equivalent method) there are three independent physical equations for mass, momentum and energy. Two-equation models are in reality two-variable models. For such a model to be theoretically coherent, one of the three physical equations has to be redundant with the two others. If not, an incoherence arises from the fact that there is three independent physical equations for only two independent variables.

In hydraulics of open-channel flows, the Saint-Venant equations with friction and a flat velocity profile are theoretically coherent since the energy equation is redundant with the mass and momentum equation. In other terms the averaged kinetic energy equation can be derived from the averaged mass and momentum equations. This property is true as long as the deviations from a flat velocity profile are negligible (Richard et al. [START_REF] Richard | Consistent equations for open-channel flows in the smooth turbulent regime with shearing effects[END_REF]). In the case of viscous thin films down an inclined plane, assuming no deviation from a parabolic velocity profile implies that τxz (0) Ũ = 3 Ũ 2 / h but this leads to models (like Shkadov's model) that are not consistent.

General form of two-equation models

Existing two-equation models are derived either from the averaging over the fluid depth of the mass and kinetic energy equations (energy integral method or EIM) either from the averaging of the mass and momentum equations (momentum integral method or MIM). Some models were not originally presented as derived by a standard depthaveraging method but turn out to be equivalent. For example the two-equation model of Ruyer-Quil and Manneville [START_REF] Ruyer-Quil | Improved modeling of flows down inclined planes[END_REF] is equivalent to the EIM and gives, at the first order, the same model as the models of Usha and Uma [START_REF] Usha | Modeling of stationary waves on a thin viscous film down an inclined plane at high Reynolds numbers and moderate Weber numbers using energy integral method[END_REF], Mudunuri and Balakotaiah [START_REF] Mudunuri | Solitary waves on thin falling films in the very low forcing frequency limit[END_REF], Abderrahmane and Vatistas [START_REF] Abderrahmane | Improved two-equation model for thin layer fluid flowing down an inclined plane problem[END_REF] and Luchini and Charru [START_REF] Luchini | Consistent section-averaged equations of quasi-one-dimensional laminar flow[END_REF].

Taking κ = 0, most consistent two-equation models write at the first order with the mass equation ( 9) and one of the two following equations, the momentum equation for MIM models, written under the general form

∂ h Ũ ∂ t + ∂ ∂ x h Ũ 2 + F 1 h, Ũ + α 1 h2 cos θ 2Fr 2 = C 1 εRe λ h - 3 Ũ h ( 18
)
or the energy equation for EIM models that can be written

∂ ∂ t h Ũ 2 2 + E 2 h, Ũ + α 3 h2 cos θ 2Fr 2 + ∂ ∂ x h Ũ 3 2 + F 2 h, Ũ + α 2 h2 Ũ cos θ Fr 2 = C 2 Ũ εRe λ h - 3 Ũ h (19)
In these expressions, F 1 , E 2 and F 2 are functions of the variables h and Ũ but do not depend on θ whereas C 1 , C 2 , α 1 , α 2 and α 3 are dimensionless constants. The consistency of equations ( 18) and [START_REF] Bresch | A generalization of the quantum Bohm identity: Hyperbolic CFL condition for Euler-Korteweg equations[END_REF] with equations ( 16) and [START_REF] Chakraborty | Extreme solitary waves on falling liquid films[END_REF] imposes some relations between these constants and some conditions on these functions which are examined below.

The MIM model of Lavalle et al. [START_REF] Lavalle | A numerical reduced model for thin liquid films sheared by a gas flow[END_REF] has

F 1 = 2λ 2h5 225 ; C 1 = α 1 = 1 ( 20
)
whereas the MIM model of Noble and Vila [START_REF] Noble | Thin power-law film flow down an inclined plane: consistent shallow-water models and stability under large-scale perturbations[END_REF] has

F 1 = λ 2h5 45 ; C 1 = α 1 = 5 6 (21)
Note that the Shkadov model has F 1 = h Ũ 2 /5, C 1 = α 1 = 1 but is not consistent (see for example Richard et al. [START_REF] Richard | A three-equation model for thin films down an inclined plane[END_REF]).

The EIM models (Ruyer-Quil and Manneville [START_REF] Ruyer-Quil | Improved modeling of flows down inclined planes[END_REF], Usha and Uma [START_REF] Usha | Modeling of stationary waves on a thin viscous film down an inclined plane at high Reynolds numbers and moderate Weber numbers using energy integral method[END_REF], Mudunuri and Balakotaiah [START_REF] Mudunuri | Solitary waves on thin falling films in the very low forcing frequency limit[END_REF], Abderrahmane and Vatistas [START_REF] Abderrahmane | Improved two-equation model for thin layer fluid flowing down an inclined plane problem[END_REF] and Luchini and Charru [START_REF] Luchini | Consistent section-averaged equations of quasi-one-dimensional laminar flow[END_REF]) were originally written under different forms but it can be shown that all can be written at the first order of approximation in the same conservative form with

E 2 = h Ũ 2 10 ; F 2 = 19 70 h Ũ 3 ; α 3 = α 2 = C 2 = 1 (22)
The various EIM models proposed in the literature differs from each other only in the second-order or even thirdorder terms.

To study the consistency of the models the flow variables are expanded as formal power series of the small parameter ε:

ũ = ũ0 + εũ 1 + O ε 2 (23)
w = w0 + ε w1 + O ε 2 (24)
p = p0 + εp 1 + O ε 2 (25) τxz = τ0 + ετ 1 + O ε 2 (26)
The consistent expressions of τxz on the bottom (z = 0) and of the average velocity U = u can be found by inserting these expansions into the dimensionless equations (A.1)-(A.8). In (A.2), the leading order gives:

∂ τ0 ∂ z = -λ (27)
Using the dynamic boundary condition (A.7) τxz ( h) = O(ε) which implies that τ0 ( h) = 0, the integration of this relation gives the expression of τxz (0) at order 0:

τ0 (0) = λ h ( 28
)
Then the constitutive relation of Newtonian fluids (A.4) becomes at order zero

∂ ũ0 ∂ z = τ0 (29)
The integration of this relation gives ũ0 = λz( h -z/2) which, by taking the average of this expression, leads to the expression of the average velocity at order 0

3 Ũ0 h = λ h ( 30
)
The expressions at order 1 are similarly found although the calculations are much more complicated. Note that the derivative ∂ h/∂ t is expressed as a function of h and ∂ h/∂ x by using the depth-averaged mass equation [START_REF] Malamataris | Solitary waves on inclined films: Flow structure and binary interactions[END_REF] which is an exact equation. This yields

∂ h ∂ t = -λ h2 ∂ h ∂ x + O (ε) (31)
The expressions at order 1 are

τ1 (0) = Re h 1 3 λ 2 h3 - cos θ Fr 2 ∂ h ∂ x ( 32
) 3 Ũ1 h = Re h 2 5 λ 2 h3 - cos θ Fr 2 ∂ h ∂ x (33)
These expressions show that, while at order zero τ0 (0) = 3 Ũ0 /h, at order 1 τ1 (0) = 3 Ũ1 /h. Since the term in τxz (0) in equation ( 16) has a factor 1/ε, τ xz (0) must be expanded to the order 1 in [START_REF] Lavalle | A numerical reduced model for thin liquid films sheared by a gas flow[END_REF] to get an equation consistent at the first order. The dissipative term in the energy equation [START_REF] Chakraborty | Extreme solitary waves on falling liquid films[END_REF] writes consistently -3 Ũ 2 / h without any further correction but the friction term -τ xz (0) in the momentum equation ( 16) needs a correction with respect to the corresponding zero-order term -3 Ũ / h. For this reason, as it was pointed by Luchini and Charru [START_REF] Luchini | Consistent section-averaged equations of quasi-one-dimensional laminar flow[END_REF], the EIM can provide consistent models of order one more easily that the MIM. Since Ũ = Ũ0 + ε Ũ1 + O(ε 2 ), the previous expressions of Ũ0 and Ũ1 imply that the right-hand side of ( 18) is consistent to -3C 1 U 1 /(Re h) which can be written

∂ ∂ x - 2 25 C 1 λ 2 h5 + C 1 h2 cos θ 2Fr 2 + O (ε) (34)
We can write that u 2 = U 2 + (u -U ) 2 and the calculation gives

(ũ 0 -Ũ0 ) 2 = λ 2h4 45 ( 35
)
Because the relaxation term in the right hand side of ( 16) is equal to -τ 1 (0)/Re + O(ε), the general form ( 18) is consistent to [START_REF] Lavalle | A numerical reduced model for thin liquid films sheared by a gas flow[END_REF] provided that

∂ ∂ x F 1 + C 1 2 25 λ 2 h5 + (α 1 -C 1 ) h2 cos θ 2Fr 2 = ∂ ∂ x 4λ 2h5 45 + O(ε) (36)
Given that F 1 does not depend on θ, the consistency of equation [START_REF] Noble | Stability Theory for Difference Approximations of Euler-Korteweg Equations and Application to Thin Film Flows[END_REF] implies that

α 1 = C 1 (37)
The same method is applied to the energy equation. The right-hand side of ( 19) is consistent to

∂ ∂ x - 2 105 C 2 λ 3 h7 + 1 12 C 2 λ h4 cos θ Fr 2 (38)
Given that E 2 and F 2 do not depend on θ, the consistency of equation [START_REF] Bresch | A generalization of the quantum Bohm identity: Hyperbolic CFL condition for Euler-Korteweg equations[END_REF] implies that

α 2 = C 2 + 3α 3 4 (39)
We want to determine whether two-equation models can satisfy both the averaged momentum equation and the averaged kinetic energy equation under the general forms [START_REF] Noble | Stability Theory for Difference Approximations of Euler-Korteweg Equations and Application to Thin Film Flows[END_REF] and [START_REF] Bresch | A generalization of the quantum Bohm identity: Hyperbolic CFL condition for Euler-Korteweg equations[END_REF]. The additional constraint that the energy of the system is consistent to the physical energy

e = ũ2 2 + h cos θ 2Fr 2 (40)
could also be imposed.

The discussion in §2.1 shows that there is no physical reason to have a redundancy between the averaged momentum and energy equations. This redundancy is needed only to get a theoretically coherent two-equation model.

Galilean invariance

The principle of Galilean invariance states that the laws of mechanics are the same in all Galilean reference frames. We suppose here that the reference frame R 0 of the bottom wall is Galilean. Let's consider another Galilean reference frame R which is in a constant rectilinear translation in the Ox direction at a velocity v e with respect to R 0 . In R the bottom wall has a velocity v b = -v e and any velocity v R is related to the corresponding velocity in R 0 , denoted by v R0 , by

v R = v R0 + v b . In R, the no-slip boundary condition at the bottom becomes u R (z = 0) = v b . The average velocity of the Nusselt flow in R is u R N = gh 2 sin θ 3ν + v b (41)
Defining the Reynolds number by

Re = h N u R0 N /ν = h N ( u R N -v b )/ν and the Froude number Fr = u R0 N / √ gh n = ( u R N -v b )/
√ gh n , we get λ = 3. The average velocity u R0 N is used to define the dimensionless quantities. In these conditions, the dimensionless average velocity at order 0 in R is ũR

0 = λ h2 3 + ṽb (42)
The material derivative Df /Dt of any quantity f is is defined as ∂f /∂t + U ∂f /∂x. It is Galilean invariant i.e. (Df /Dt) R = (Df /Dt) R0 (this property is not true for the partial derivative with respect to time). The mass conservation equation ( 9) obviously satisfies unconditionally the Galilean invariance. Since the mass equation implies that

∂ h Ũ ∂ t + ∂ h Ũ 2 ∂ x = h D Ũ D t ( 43
)
the momentum equation ( 18) is Galilean invariant if F 1 is Galilean invariant (provided that the mass conservation equation is satisfied). In R, equation ( 18) is written (superscripts R are omitted) with the Galilean-invariant form

∂ h Ũ ∂ t + ∂ ∂ x h Ũ 2 + F 1 + α 1 h2 cos θ 2Fr 2 = C 1 εRe λ h - 3 h Ũ -ṽb (44)
The function F 1 is Galilean invariant if it is a function of h only but also if F 1 depends on Ũ since [START_REF] Kliakhandler | Viscous beads on vertical fibre[END_REF] implies that in this case F 1 depends in fact on Ũ -ṽb . For example, [START_REF] Liu | Solitary wave dynamics of film flows[END_REF] and ( 21) are consistent Galilean invariant MIM models but

F 1 = h 5 Ũ -ṽb 2 , C 1 = α 1 = 5 6 (45)
or

F 1 = 2 h 25 Ũ -ṽb 2 , C 1 = α 1 = 1 (46)
are also consistent Galilean invariant MIM models. The dependance on Ũ -ṽb in the models ( 45) and ( 46) means physically that the Galilean-invariant part of the momentum flux depends on the shearing since Ũ = ṽb implies a shearing effect. Note that the flux in the averaged momentum balance equation ( 16) is the sum of a Galilean-invariant term and a term which is not Galilean invariant. Writing

ũ (x, z, t) = Ũ (x, t) + ũ (x, z, t) ( 47
)
where ũ is the deviation of the velocity to its average value, the momentum flux in equation ( 16) becomes

h Ũ 2 + h ũ 2 + h2 cos θ 2Fr 2 (48)
In this expression, h Ũ 2 is a non-Galilean invariant term and h ũ 2 + h2 cos θ/(2Fr 2 ) is Galilean invariant. The general form [START_REF] Noble | Stability Theory for Difference Approximations of Euler-Korteweg Equations and Application to Thin Film Flows[END_REF] respects the physical distinction between these two parts. The necessity to find a consistent way to write the bottom friction -τ xz (0), which is Galilean invariant, explains that F 1 is not necessarily consistent to h ũ 2 and that α 1 , and thus C 1 (see equation ( 37)), are not necessarily equal to 1. We can conclude that all consistent MIM models are Galilean invariant provided that the mass conservation ( 9) is satisfied which is always the case.

In the same way, the energy ( 40) is the sum of Ũ 2 /2 which is not Galilean invariant and which assumes the role of the kinetic energy of the model and of ( ũ 2 + h2 Fr -2 cos θ)/2 which is Galilean invariant and which can therefore be treated as the potential energy (see Gavrilyuk and Perepechko [START_REF] Gavrilyuk | Variational approach to constructing hyperbolic models of two-velocity media[END_REF] for a discussion on what are the kinetic energy and the potential energy of a system in the context of a variational method). It follows that, in the energy of the model Ũ 2 /2 + E 2 + α 3 h2 Fr -2 cos θ/2, E 2 is Galilean invariant. The same discussion as for F 1 can be followed for E 2 with the conclusion that E 2 should be a function of h and ( Ũ -ṽb ) 2 . This means physically that the potential energy could depend on shearing i.e. on the relative velocity of the liquid with respect to the bottom. A dependance of the potential energy with respect to a relative velocity is already known for some systems, in particular in the case of a two-liquid medium (Gavrilyuk and Perepechko [START_REF] Gavrilyuk | Variational approach to constructing hyperbolic models of two-velocity media[END_REF]).

However the discussion on the Galilean invariance of the averaged work-energy theorem is more complex than for the averaged momentum balance equation. Starting from the general form [START_REF] Bresch | A generalization of the quantum Bohm identity: Hyperbolic CFL condition for Euler-Korteweg equations[END_REF], we write

F 2 = Ũ F 3 + F 4 ( 49
)
where F 3 and F 4 are Galilean invariant. Let's assume that equation ( 19) is written in R 0 . The Galilean velocityaddition formula is Ũ R0 = Ũ R -ṽb . We have

F R0 2 = Ũ R0 F 3 + F 4 and F R 2 = ( Ũ R -v b )F 3 + F 4 .
The relations between the partial derivatives of any quantity G in R 0 and R are

∂G ∂x R0 = ∂G ∂x R = ∂G ∂x (50)
and

∂G ∂t R0 = ∂G ∂t R + v b ∂G ∂x (51)
The averaged work-energy theorem in R can be written

∂ ∂ t h( Ũ R ) 2 2 + E 2 + α 3 h2 cos θ 2Fr 2 R + ∂ ∂ x h( Ũ R ) 3 2 + F R 2 + α 2 h2 Ũ R cos θ Fr 2 = C 2 Ũ R εRe λ h - 3 h Ũ R -ṽb + ṽb ∂ h Ũ R ∂ t R + ∂ ∂ x h( Ũ R ) 2 + F 3 -E 2 + (2α 2 -α 3 ) h2 cos θ 2Fr 2 - C 2 εRe λ h - 3 h Ũ R -ṽb (52)
Consequently, the averaged work-energy theorem is Galilean invariant if the expression between brace brackets in (52) is equal to zero i.e. if the averaged momentum balance equation ( 44) is satisfied with the following conditions:

   F 1 = F 3 -E 2 α 1 = 2α 2 -α 3 C 1 = C 2 (53)
Taking into account the consistency conditions ( 37) and [START_REF] Benjamin | Wave formation in laminar flow down an inclined plane[END_REF] implies that all dimensionless coefficients must be equal and will be thereafter denoted by C 1

C 1 = C 2 = α 1 = α 2 = α 3 (54)
This is a general result: the work-energy theorem is Galilean invariant provided that the momentum balance equation is satisfied, the mass conservation equation being a prerequisite for the Galilean invariance of both the momentum and kinetic energy equations. This result was proved in the case of the Saint-Venant equations and the one-dimensional Teshukov equations (Teshukov [START_REF] Teshukov | Gas-dynamic analogy for vortex free-boundary flows[END_REF]) by Ostapenko [START_REF] Ostapenko | Conservation laws of shallow water theory and the Galilean relativity principle[END_REF] who established a hierarchy of the conservation laws for the Saint-Venant equations: the mass conservation is unconditionally Galilean invariant; the momentum balance equation is Galilean invariant provided that the mass conservation equation is satisfied; the energy balance equation is Galilean invariant if the mass and momentum equations are satisfied. A similar result is obtained.

Compatibility between the energy and momentum equations

The conditions of compatibility between the energy and the momentum equations for two-equation models are derived in Appendix B. Adding the Galilean-invariance conditions (54) to the compatibility condition (B.8) leads to

∂E 2 ∂ Ũ = 0 (55)
Consequently, E 2 depends only on h. In all models where E 2 depends also on U , even under the form U -v b , there is a contradiction. For example, EIM models [START_REF] Ostapenko | Justification of shallow-water theory[END_REF] are not Galilean invariant since they are not compatible with the momentum equation.

The second compatibility condition (B.14) between the energy and momentum equations added to the Galileaninvariance conditions (54) writes

∂F 1 ∂ h = h ∂ 2 E 2 ∂ h2 (56)
Since E 2 depends only on h, the integration of this relation yields

F 1 = hE 2 ( h) -E 2 ( h) + K( Ũ ) (57)
where dE 2 /d h is denoted by E 2 and where K is a function of Ũ that does not depend on h. A study of the consistency shows that F 1 is consistent to α 4 λ h5 where α 4 is a dimensionless constant. At order 0, Ũ is consistent to λ h2 /3. To be consistent, K should be on the form α 5 Ũ 5/2 / √ λ where α 5 is a dimensionless constant. No existing model features such a term. Accordingly models where the dependency of F 1 with h and Ũ is written as α 6 h Ũ 2 , where α 6 is a dimensionless constant, are explicitly excluded.

Further a term in Ũ 5/2 behaves badly if U becomes negative. All these models are derived assuming that the Nusselt velocity profile is slightly perturbed. In these conditions, the average velocity U has the sign of λ which is the sign of the slope. Thus a term with the form U/λ remains positive if the perturbations to the Nusselt profile are small. In practice however the average velocity can become negative even with a positive slope. This is the case with the phenomenon of flow reversal which can be observed in the first depth minimum of the capillary waves (Tihon et al. [START_REF] Tihon | Solitary waves on inclined films: their characteristics and the effects on wall shear stress[END_REF], Malamataris et al. [START_REF] Malamataris | Solitary waves on inclined films: Flow structure and binary interactions[END_REF]). An artificial solution would be to write |U | but otherwise, these terms must be equal to zero for the model to be acceptable.

In practice, we can assume that K( Ũ ) = 0 and thus

∂F 1 ∂ Ũ = 0 (58)
All models where F 1 depends on Ũ are incompatible with the energy equation. For example, the MIM models (45) and ( 46) are Galilean invariant but they don't admit an energy equation.

The consistency of the right-hand side of momentum equation [START_REF] Seppecher | Etude des conditions aux limites en théorie du second gradient: cas de la capillarité[END_REF] and of the other terms of the momentum equation ( 18) combined with the condition (58) leads to

F 1 = 1 9 - C 1 10 4 5 λ 2 h5 + O(ε) (59)
The positivity of F 1 implies that

C 1 < 10 9 ( 60
)
The integration of equation ( 56) then gives

E 2 = 1 9 - C 1 10
λ 2h5 5 + O(ε) (61)
with the same positivity condition as (60). The two compatibility conditions (B.9) and (B.10) reduce to

∂F 2 ∂ h = Ũ ∂E 2 ∂ h + ∂F 1 ∂ h (62)
and

∂F 2 ∂ Ũ = h ∂E 2 ∂ h (63)
This allows to find the expression of F 2

F 2 = 1 9 - C 1 10 λ 2 h5 Ũ + O(ε) (64)
The conclusion of this study is that consistent Galileaninvariant two-equation models which are compatible both with the momentum equation and the energy equation are MIM models whose E 2 , F 1 and F 2 flux terms are of the form (61), (59) and (64) respectively. This family of MIM models will be called thereafter momentum integral with mathematical energy (MIME) models. The general conservative form of these models can be obtained by writing

E 2 = λ 2 h 5 4K ; F 1 = λ 2 h 5 K ; F 2 = 5λ 2 h 5 U 4K ( 65
)
where the dimensionless number K is

K = 225 20 -18C 1 (66)
The models of Lavalle et al. [START_REF] Lavalle | A numerical reduced model for thin liquid films sheared by a gas flow[END_REF] and of Noble and Vila [START_REF] Noble | Thin power-law film flow down an inclined plane: consistent shallow-water models and stability under large-scale perturbations[END_REF] are MIME models with respectively

C 1 = 1 ; E 2 = λ 2h5 450 ; F 1 = 2λ 2 h 5 225 ; F 2 = λ 2h5 Ũ 90 (67)
and

C 1 = 5 6 ; E 2 = λ 2h5 180 ; F 1 = λ 2 h 5 45 ; F 2 = λ 2h5 Ũ 36 (68)
If, in addition, we impose the constraint that the energy of the model is consistent to the physical energy, there are more conditions to satisfy. The calculation of the consistency of the energy e (40) at order 0 shows that e is consistent to

2λ 2h4 30 + h cos θ 2Fr 2 + O (ε) (69)
This implies that E 2 is consistent to λ 2h5 /90. From the expression (61), we get C 1 = 5/9. On the other hand, the consistency of the term in cos θ implies that C 1 = 1. These two results are of course contradictory and no consistent Galilean-invariant two-equation model has an energy consistent to the physical energy while being compatible with the momentum and energy equations except for one model which has inconveniently terms in Ũ 5/2 and Ũ 7/2 in the fluxes (see Appendix C).

Capillary energy

In this section, all tilde symbols are dropped to lighten the notations. The conservative hyperbolic part of the model admits a mathematical entropy if there is a convex function η of the conservative variables h and q = hU that satisfies the additional conservation law

∂η ∂t + ∂Φ ∂x = 0 ( 70
)
where Φ is a function of h and q. For Galilean-invariant MIM models which are compatible with the energy equation, thus which satisfy the condition ∂F 1 /∂U = 0, the mathematical entropy is the energy of the system. Other MIM models admit a mathematical entropy but this entropy is not homogeneous to an energy. For example, a MIM model with C 1 = 5/6 and F 1 = hU 2 /5 admits the mathematical entropy η = q 2 h 4/5 + 125 198

h 11/5 cos θ Fr 2 (71)
which is not an energy (see Appendix D). In the classical two-equation models of isentropic gas dynamics, the basic conservation laws are the mass and momentum conservation laws and the energy conservation law is a convex extension of the basic laws which implies that the energy is the mathematical entropy of the system. In the case of the two-equation models of thin film flows, this property is only obtained with the hyperbolic part of MIME models. From this point of view these MIME models are more methodologically grounded than the EIM models or other MIM models. The capillary terms in the averaged momentum equation [START_REF] Lavalle | A numerical reduced model for thin liquid films sheared by a gas flow[END_REF] and in the averaged work-energy theorem [START_REF] Chakraborty | Extreme solitary waves on falling liquid films[END_REF] exhibit the structure of an internal capillarity (Casal [START_REF]La capillarité interne[END_REF], Casal and Gouin [START_REF] Casal | Relation entre l'équation de l'énergie et l'équation du mouvement en théorie de Korteweg de la capillarité[END_REF]). The capillary term in ( 16) can be written in the conservative form

κ Fr 2 h ∂ 3 h ∂x 3 = κ Fr 2 ∂ ∂x h ∂ 2 h ∂x 2 - 1 2 ∂h ∂x 2 (72)
and the capillary term in [START_REF] Chakraborty | Extreme solitary waves on falling liquid films[END_REF] writes

κ Fr 2 hU ∂ 3 h ∂x 3 = - ∂ ∂t 1 2 κ Fr 2 ∂h ∂x 2 - ∂ ∂x U 2 κ Fr 2 ∂h ∂x 2 + κ Fr 2 ∂ ∂x hU ∂ 2 h ∂x 2 - U 2 ∂h ∂x 2 -h ∂h ∂x ∂U ∂x (73)
The fluid depth of the depth-averaged models is analogous to the density in the equations of compressible fluid mechanics. It is then possible to define a capillary specific energy

e c = κ Fr 2 1 2h ∂h ∂x 2 (74)
This capillary energy depends on the gradient of h (here in the one-dimensional case) in exactly the same manner as the capillary energy of the fluids endowed with internal capillarity depends on the gradient of the density (the description of capillarity with such an energy is due to Korteweg [START_REF] Korteweg | Sur la forme que prennent les équations du mouvements des fluides si l'on tient compte des forces capillaires causées par des variations de densité considérables mais connues et sur la théorie de la capillarité dans l'hypothèse d'une variation continue de la densité[END_REF]). It is then possible to defined a generalized stress tensor (that reduces to a scalar in this onedimensional case) for the internal capillary effects that can be written

σ = κ Fr 2 h ∂ 2 h ∂x 2 - 1 2 ∂h ∂x 2 (75)
There is an extra term h(∂h/∂x)(∂U/∂x) in the flux of the energy equation that corresponds to the interstitial working of Dunn and Serrin [START_REF] Dunn | On the thermodynamics of interstitial working[END_REF] although Casal and Gouin [START_REF] Casal | Relation entre l'équation de l'énergie et l'équation du mouvement en théorie de Korteweg de la capillarité[END_REF] showed that this term appears naturally in the framework of the theory of the second gradient developed by Germain [START_REF] Germain | La méthode des puissances virtuelles en mécanique des milieux continus[END_REF] (see also Gatignol and Seppecher [START_REF] Gatignol | Modelisation of fluid-fluid interfaces with material properties[END_REF], Seppecher [START_REF] Seppecher | Etude des conditions aux limites en théorie du second gradient: cas de la capillarité[END_REF]).

The capillary energy can be written W 2 /2 by defining the variable

W = κ Fr 2 h ∂h ∂x (76)
This variable can be treated as an independant variable if we add an initial condition such that (76) is satisfied at t = 0 (Gavrilyuk and Gouin [START_REF] Gavrilyuk | Trends in applications of mathematics to mechanics[END_REF]). The variable W satisfies an equation derived from the mass conservation. This is at the basis of an extended model that can be solved by a numerical scheme for which the nonlinear stability can be proved in the case of a Rusanov Riemann solver (Noble and Vila [START_REF] Noble | Stability Theory for Difference Approximations of Euler-Korteweg Equations and Application to Thin Film Flows[END_REF]). Without this treatment, nothing can be proved for the nonlinear stability and upwind schemes are always linearly unstable whereas centred schemes are linearly stable if they are implicit and linearly unstable if they are explicit (Noble and Vila [START_REF] Noble | Stability Theory for Difference Approximations of Euler-Korteweg Equations and Application to Thin Film Flows[END_REF]).

The possibility to derive an extended model for a given two-equation model and the proof of the nonlinear stability of the numerical scheme require the existence of a capillary energy. A MIM model admits a capillary energy if and only if ∂F 1 /∂U = 0 (see Appendix D). This is the same condition as the compatibility condition between the momentum and energy equations (the dubious models with terms in U 5/2 being excluded). In practice, we can assert that a MIM model admits a capillary energy if and only if it satisfies the energy equation. The specific capillary energy is (Appendix D)

e c = C 1 κ Fr 2 1 2h ∂h ∂x 2 (77)
This capillary energy is equal to the physical capillary energy (74) if C 1 = 1, which corresponds to the model of Lavalle et al. [START_REF] Lavalle | A numerical reduced model for thin liquid films sheared by a gas flow[END_REF]. However, even if C 1 = 1, an extended model can be derived with the extra variable

W = C 1 κ hFr 2 ∂h ∂x (78)
The equation of W can be derived from the mass equation.

It can be written (Noble and Vila [START_REF] Noble | Stability Theory for Difference Approximations of Euler-Korteweg Equations and Application to Thin Film Flows[END_REF])

∂hW ∂t + ∂hU W ∂x = - ∂ ∂x C 1 κ Fr 2 h 3/2 ∂U ∂x (79)
If ∂F 1 /∂U = 0, the quantity W can be defined as in (78) and it satisfies the same equation (79) but no capillary energy can be found. More generally, no function of h and ∂h/∂x corresponding to capillarity can be added to the mathematical entropy for a MIM model with ∂F 1 /∂U = 0. The proof is given in Appendix D. This means that not only the structure of an internal capillarity is lost in the model but no generalization to a purely "mathematical" capillary energy can be found. The extended model of Noble and Vila [START_REF] Noble | Stability Theory for Difference Approximations of Euler-Korteweg Equations and Application to Thin Film Flows[END_REF] is not applicable and the theorems of nonlinear stability cannot be proved. Fulfilling the condition (58) is thus of major importance for a two-equation model for a reliable treatment of capillarity.

Speed of solitary waves for the hyperbolic system

The nonlinear properties of models can be tested on the calculation of the speed of solitary waves down an inclined plane (Chakraborty et al. [START_REF] Chakraborty | Extreme solitary waves on falling liquid films[END_REF]). All consistent models give approximately the same results when the Reynolds number is small but they differ widely one from each other if the Reynolds number becomes high (here "high" can be only of the order of 30). A direct numerical simulation (DNS) calculation made by Chakraborty et al. [START_REF] Chakraborty | Extreme solitary waves on falling liquid films[END_REF] allows to compare the accuracy of each model in the high Reynolds number limit. The four-equation model of Ruyer-Quil and Manneville [START_REF] Ruyer-Quil | Improved modeling of flows down inclined planes[END_REF] gives very accurate values whereas some models do not give any value at all because they do not allow solitary waves solutions beyond a certain limit value of the Reynolds number (Demekhin et al. [START_REF] Demekhin | Mathematical models of the theory of viscous liquid films[END_REF], Ruyer-Quil and Manneville [START_REF] Ruyer-Quil | Improved modeling of flows down inclined planes[END_REF]).

At a Reynolds number value high enough, the speed of the system depends only on the hyperbolic part of the equations of the model. The viscous diffusive terms and the capillary dispersive terms become negligible if Re is high. To calculate the speed of solitary waves for high Reynolds numbers, it is enough to take the first-order model with κ = 0.

We consider first MIM models. The model is rewritten in dimensional form. To lighten the notation, we will write simply F 1 in the dimensional equations. Note that, if F1 is the dimensionless form (used in the previous sections without tilde), then

F 1 = h N u N 2 F1
. The mass and momentum equations are

∂h ∂t + ∂hU ∂x = 0 (80) ∂hU ∂t + ∂ ∂x hU 2 + F 1 + C 1 2 gh 2 cos θ = C 1 gh sin θ -3ν U h (81)
The study of the solitary waves is done in the spirit of Dressler's theory of roll waves (Dressler [START_REF] Dressler | Mathematical solution of the problem of rollwaves in inclined opel channels[END_REF]). We look for periodic discontinuous solutions which are stationary in the Galilean reference frame propagating at the speed c of the waves. A solitary wave is then found as the limit of these periodic waves when the period goes to infinity.

Thus the solutions depend only on the variable ξ = x -ct.

The mass equation (80) leads to the conservation of the relative flow discharge m

h (U -c) = m = constant < 0 ( 82
)
Then the momentum equation yields

dh dξ = C 1 gh sin θ -3νU/h ∂F 1 ∂h + C 1 gh cos θ - m h 2 ∂F 1 ∂U - m 2 h 2 (83)
Following closely the discussion in Dressler [START_REF] Dressler | Mathematical solution of the problem of rollwaves in inclined opel channels[END_REF], there is a critical point where the denominator vanishes. The characteristics of the systems are given in Appendix E. The expression (E.3) shows that if

∂F 1 /∂U = 0 (condition ( 58
∂F 1 ∂U (85) a = 1 4h 2 ∂F 1 ∂U 2 + ∂F 1 ∂h + C 1 gh cos θ (86)
In this case, a is the velocity of the surface waves and U + a (not U ) is an effective fluid velocity that depends on shearing (see §3). A generalized relative Froude number (in the reference frame of the wave) can be defined as

Fr * = (c -U -a )/a instead of the usual Fr = (c -U )/a.
This generalized relative Froude number is equal to 1 at the critical point.

As for Dressler's roll waves, the numerator must vanish at the critical point for the slope to remain finite. This gives the condition, where the subscripts C stand for the quantities at the critical point

U C = gh 2 C sin θ 3ν (87)
The equation ( 83) is written in dimensionless form with the following scaling:

H = h h C X = ξ h C U = U U C C = c U C M = m h C U C F 1 = F 1 h C U 2 C (88)

Moreover, we define a critical Reynolds number Re

C = h C U C /ν, a critical Froude number Fr C = U C / √ gh C and the dimensionless number λ C = Re C sin θ/Fr 2
C . The condition (87) implies that λ C = 3. The constant M can be calculated at the critical point:

M = - 1 2 ∂F 1 ∂U C - 1 4
∂F 1 ∂U 2 C + ∂F 1 ∂H C + C 1 cos θ F 2 C ( 89
)
The dimensionless equation can be written

dH dX = C 1 Re C λ C H -3 U/H ∂F 1 ∂H + C 1 H cos θ Fr 2 C - M H 2 ∂F 1 ∂U - M 2 H 2 (90)
For a solitary wave, when X → -∞, the solution goes to the fixed point

(H ∞ , U ∞ ) such that λ C H 2 ∞ = 3U ∞ . This reduces to U ∞ = H 2
∞ . From (82), we write U = C + M/H and C = 1 -M since M is a constant. This yields the equation H 3

∞ -(1-M )H ∞ -M = 0 which has the evident solution H ∞ = 1 which should be discarded. We thus obtain the equation

H 2 ∞ + H ∞ + M = 0. The positive solution is H ∞ = 1 2 -1 + √ 1 -4M (91)
The velocity of the solitary wave c normalized by the speed of the kinematic waves 3U

∞ is equal to C/(3U ∞ ) or C/(3H 2 ∞ ). Its expression is c 3U ∞ = 4 3 1 -M -1 + √ 1 -4M 2 (92)
The Froude number defined by Chakraborty et al. [START_REF] Chakraborty | Extreme solitary waves on falling liquid films[END_REF], that we will denote by F r , is F r = 3Fr ∞ / √ cos θ where

Fr ∞ = U ∞ / √ gh ∞ . It follows that F r -2 = 1 9 cos θ Fr 2 C 1 H 3 ∞ (93)
A vertical slope corresponds to F r -2 = 0 and the instability threshold to F r -2 = 0.4. The velocity of the solitary wave c/(3U ∞ ) is calculated for values of F r -2 equal to 0 ; 0.1 ; 0.2 ; 0.3 and 0.4 (in the latter case, it should be 1). The results are presented in Table 1. The value of C 1 and the expression of F 1 are given in each case. The results of the DNS of Chakraborty et al. [START_REF] Chakraborty | Extreme solitary waves on falling liquid films[END_REF] and those of the four-equation model of Ruyer-Quil and Manneville [START_REF] Ruyer-Quil | Improved modeling of flows down inclined planes[END_REF] are given for comparison.

We consider now the EIM model [START_REF] Ostapenko | Justification of shallow-water theory[END_REF]. The same approach is followed. Writing E 2 = hE 3 and

F 2 = hU E 3 + F 3 , we get dh dξ = gh sin θ -3νU/h Ξ ( 94
)
where slopes with an excellent precision, taking the DNS values as reference. For a vertical wall, the velocity of the solitary wave calculated with this model is (Ruyer-Quil and Manneville [START_REF] Ruyer-Quil | On the speed of solitary waves running down a vertical wall[END_REF])

Ξ = gh cos θ + m U ∂E 3 ∂h + 1 U ∂F 3 ∂h - m 2 h 2 U ∂E 3 ∂U - m h 2 U ∂F 3 ∂U - m 2 h 2 (95) with E 3 = U 2 /
c 3U ∞ = 1 + 1 √ 6 + 1 2 + 2 3 2.556 (97)
For this model, the shape factor defined by

α S = u 2 U 2 (98)
is constant and equal to 6/5 which is the value of a parabolic velocity profile. This model is consistent while being reminiscent of the Shkadov model. It will be called thereafter the corrected Shkadov model. Its major drawback is that, condition (58) being unfulfilled, it does not satisfy the energy equation and does not admit a capillary energy. Consequently, it does not admit an extended model in the manner of Noble and Vila [START_REF] Noble | Stability Theory for Difference Approximations of Euler-Korteweg Equations and Application to Thin Film Flows[END_REF] and nothing can be proved about the nonlinear stability of the numerical scheme.

A MIME model gives exactly the same value (97) of c/(3U ∞ ) for a vertical wall if

C 1 = 31 30 - √ 6 45 0.9789 (99)
The value of K is then K = 315 -90 √ 6 94.55. This model will be called thereafter the equivalent MIME model. For non-vertical walls, the velocities calculated with the equivalent MIME model are slightly less precise because the dependance of the velocity on cos θ in the expression (89) is the same in all MIM models except for the value of C 1 . The difference between the equivalent MIME model and the corrected Shkadov model with regard to the non-vertical walls lies only in the difference between the values of C 1 . The relative deviation on the velocity reaches 7.2% for F r -2 = 0.3. On the other hand, the equivalent MIME model satisfies the energy equation and admits a capillary energy and an extended model to treat capillarity. It has thus a better mathematical structure and its numerical resolution is more reliable. For this model, the shape factor is not constant.

Other variants give values which deviate more or less from the DNS results. The Galilean-invariant compatible model (see Appendix C) predicts velocities which are far too low. The EIM model gives values slightly too high with a relative deviation compared to the DNS of 7% for a vertical slope, decreasing to 2.5% for F r -2 = 0.3. The MIME model with C 1 = 1 (Lavalle et al. [START_REF] Lavalle | A numerical reduced model for thin liquid films sheared by a gas flow[END_REF]) has a relative deviation of 10% on a vertical wall, increasing to 15% for F r -2 = 0.3 due to the high value of C 1 in the term in cos θ. In practice, most applications of thin film flows pertaining to vertical walls, the accuracy on the value of the velocity on a vertical wall is the most important. Therefore, the corrected Shkadov model and the equivalent MIME model, which give exactly the same value for a vertical wall, will be favoured in the following.

It can be noted that the values of the velocities are very sensitive to small variations of C 1 . These velocities are calculated outside the validity domain of the asymptotic methods used to derive the models. As a consequence, these values depend on the extrapolation of the models outside their validity domain, which explains the huge differences between a model and another.

Diffusive terms and linear properties

The viscous terms, particulary the diffusive terms, are important to consider for laminar flows of viscous thin films. These diffusive terms can be written in different ways. In the averaged momentum balance equation ( 16), the linear terms of O(ε/Re) are consistent to (see Richard et al. [START_REF] Richard | A three-equation model for thin films down an inclined plane[END_REF])

ε Re ∂ ∂x 15 4 h ∂U ∂x (100)
In this section, all tilde symbols are dropped to lighten the notations. In the averaged energy equation ( 17), the linear terms of O(ε/Re) are consistent to (Richard et al. [START_REF] Richard | A three-equation model for thin films down an inclined plane[END_REF])

ε Re ∂ ∂x 9 2 hU ∂U ∂x (101)
In the EIM model of Ruyer-Quil and Manneville [START_REF] Ruyer-Quil | Improved modeling of flows down inclined planes[END_REF] and of Usha and Uma [START_REF] Usha | Modeling of stationary waves on a thin viscous film down an inclined plane at high Reynolds numbers and moderate Weber numbers using energy integral method[END_REF], the linear diffusive terms can be written (neglecting all nonlinear terms)

Note that the last two terms of this expression write the same way if the equations are linearized. For the linear properties, the linear diffusive terms can be written either with one term having a second derivative of U or with two terms, one having a second derivative of U and another one having a second derivative of h. In a general way, both for MIM or for EIM models, we can replace consistently

∂ ∂x h ∂U ∂x (104) by α ∂ ∂x h ∂U ∂x + 2 (1 -α) ∂ ∂x U ∂h ∂x ( 105
)
where α is a number. In the second term of the right-hand side of this expression, λh 2 /3 (or h 2 by taking λ = 3) can replace U without any effect on the linear properties. In the energy equation, these expressions are multiplied by U . We consider first MIME models. The general form of the momentum equation of these models can be written

∂hU ∂t + ∂ ∂x hU 2 + λ 2 h 5 K + C 1 h 2 cos θ 2Fr 2 = C 1 εRe λh - 3U h + C 1 κ Fr 2 h ∂ 3 h ∂x 3 + ε Re 9 2 C 1 α ∂ ∂x h ∂U ∂x + 2 (1 -α) ∂ ∂x U ∂h ∂x + ε Re C 1 3U h ∂h ∂x 2 (106)
The equations are linearized around the Nusselt flow solution considering small perturbations h and U around the equilibrium values:

h U = 1 1 + h U (107)
The perturbations are taken of the form

h U = B 1 B 2 e ik(x-ct) (108)
where B 1 , B 2 are amplitudes, k is the wavenumber and c the phase velocity. The linearized equations (taking ε = 1) can be written

h (1 -c) + U = 0 (109) h ik 45 K + C 1 cos θ Fr 2 - 6C 1 Re + 9C 1 Re (1 -α) k 2 + C 1 κ Fr 2 ik 3 + U ik (1 -c) + 3C 1 Re + 9 2 C 1 Re αk 2 = 0 (110)
The determinant of the system is equal to zero in order to get a non-trivial solution. This leads to the dispersion equation

3C 1 Re 3 -c -3k 2 + 3 2 k 2 α (3 -c) + i k 1 - 45 K -2c + c 2 - 3C 1 Re k cot θ + κk 3 sin θ = 0 (111)
The solution can be written

c = 3 + ik 6 5 Re -cot θ + O k 2 (112) 13
This solution implies that, at the long-wave limit (k → 0), the celerity of the kinematic waves is c = 3. The Nusselt flow is stable if Im(c) < 0. For k → 0, the stability occurs if Re < Re c where Re c is the critical Reynolds number. Equation (112) gives the same critical value as that calculated from the Orr-Sommerfeld equation by Benjamin [START_REF] Benjamin | Wave formation in laminar flow down an inclined plane[END_REF] and Yih [START_REF] Yih | Stability of liquid flow down an inclined plane[END_REF], which is

Re c = 5 6 cot θ (113)
For Re > Re c , the neutral (or marginal) stability condition is Im(c) = 0. With this condition, the curve k as a function of Re is the neutral stability curve. On this curve, c is real and the dispersion relation yields two real equations corresponding to the real and imaginary parts. The real part gives the expression of the velocity of the kinematic waves along the neutral stability curve

c = 3 - 6k 2 2 + 3αk 2 (114)
For k = 0, we get c = 3 as expected. If k → ∞, then c → c lim where

c lim = 3 - 2 α (115)
This limit value depends only on α i.e. the way the linear diffusive terms are written. We calculated this limit value for the Orr-Sommerfeld equations with the AUTO software (Doedel [START_REF] Doedel | AUTO07p: continuation and bifurcation software for ordinary differential equations[END_REF]) for three fluids: water, alcohol and castor oil. The Kapitza number is defined by

Ka = γ ρ gν 4 -1/3 (116)
The three fluids were chosen because they have widely different Kapitza numbers, respectively 3 500, 529 and 0.45, and thus very different capillary properties. The calculations were done for several slopes ranging from θ = 6.4 • to θ = 90 • . Some results are given in Figure 2. The graphs of the kinematic velocity as a function of the wavenumber along the neutral stability curve form a batch of curves very close one from each other whatever the slope or the Kapitza number (Figure 2a). The common limit value is

c lim = 3/2. Figure 2(b) giving ln(c -3/2) as a function of ln(k) shows that, for k 1, c -3/2 = O(k -2
). The only solution giving the same limit value c lim = 3/2 as the Orr-Sommerfeld equations when k → ∞ is

α = 4 3 (117)
However, along the neutral stability curve, for a two-equation model, the wavenumber k does not tend towards infinity when Re → ∞ but to a finite value k ∞ and the velocity of the kinematic waves tends toward c ∞ . These values can be found from the dispersion relation (111). We find c ∞ = 1 + 45/K, or The positivity condition (60) applies to this expression. This limit does not depend on α i.e. not on how the diffusive terms are written. The limit value of the wavenumber is

c ∞ = 1 + 2 1 - 9C 1 10 (118) # Model α c lim c ∞ k ∞ 1 CS/EqMIME
k 2 ∞ = 6C 1 10 -9αC 1 + 10 1 -9C 1 /10 (119)
A similar study was conducted for the corrected Shkadov model (see §6) and for the EIM models of Ruyer-Quil and Manneville [START_REF] Ruyer-Quil | Improved modeling of flows down inclined planes[END_REF] (α = 6/5) and Mudunuri and Balakotaiah [START_REF] Mudunuri | Solitary waves on thin falling films in the very low forcing frequency limit[END_REF] (α = 449/420). As in §6, the corrected Shkadov model gives exactly the same results as the equivalent MIME model. For these models, the cases α = 1, α = 6/5 and α = 4/3 were considered. The value of c lim depends only on α and thus it is the same for all MIM and EIM models that share the same value of α. The values for the different studied models are presented in Table 2. The curves of the velocity of the kinematic waves as a function of the wavenumber along the neutral stability curve for each studied model are presented in Figure 3. These curves are limited to the limit value k ∞ which is different in each case. The curve calculated from the Orr-Sommerfeld equation has no limit value of k.

We can wonder whether it is possible to find a model giving the same values of c lim , c ∞ and k ∞ as the Orr-Sommerfeld equation i.e.

c lim = 3/2, c ∞ = 3/2 and k → ∞ if Re → ∞.
The first condition implies that α = 4/3 whereas the two last conditions are satisfied if and only if C 1 = 25/24. For this model, the curve of the kinematic wave velocity on Figure 3 would not be interrupted prematurely and would follow closely the Orr-Sommerfeld curves for high values of k. Unfortunately, the nonlinear properties of this model are disappointing. The velocity c/(3U ∞ ) calculated for this model as in §6 is 2+ √ 3 3.73 on a vertical wall and 2.09 for F r -2 = 0.3 which are both too high than the DNS values with a relative deviation of about 45%. Further, values of k on the neutral stability curve higher than the value k ∞ = 1.857 of the equivalent MIME model with α = 4/3 can be reached only for very high Reynolds numbers that are usually outside the domain of application of these models.

The examination of Figure 3 shows that the value α = 4/3 gives the best result followed closely by α = 6/5. Note that the curves calculated from the Orr-Sommerfeld equation depends slightly on the Kapitza number whereas the curves calculated from the models do not depend on surface tension at all. They do not depend on the slope either whereas the Orr-Sommerfeld curve depends also slightly on the slope. The neutral stability curves are presented in Figure 4 along with the Orr-Sommerfeld curve. The effect of the value of α is studied in Figure 4(a) for the equivalent MIME model in the case of water on a vertical wall and in Figure 4(b) in the case of a glycerin-water solution for θ = 6.4 • (as in Liu and Gollub [START_REF] Liu | Solitary wave dynamics of film flows[END_REF]). Three values of α are studied: 1, 6/5 and 4/3. The value α = 4/3 gives the best results in all cases. The effect of the value of C 1 for α = 4/3 is then shown in Figure 4(c) in the case of water on a vertical wall with the values C 1 = 1, C 1 = 25/24 and the value (99) of the equivalent MIME model. The value C 1 = 1 together with α = 4/3 gives a very good agreement with the Orr-Sommerfeld curve while the curve is too high with C 1 = 25/24 and slightly too low for the equivalent MIME model. This trend is also observed with other fluids and slopes (not presented here). The corrected Shkadov model and the equivalent MIME model both for α = 4/3 and the EIM models of Ruyer-Quil and Manneville [START_REF] Ruyer-Quil | Improved modeling of flows down inclined planes[END_REF] and Mudunuri and Balakotaiah [START_REF] Mudunuri | Solitary waves on thin falling films in the very low forcing frequency limit[END_REF] are compared in Figure 4(d). The equivalent MIME model (α = 4/3) is slightly better than the corrected Shkadov model with the same value of α which, in turn, is better than the EIM models. The fluid parameters are given in Table 3.

The models of Usha and Uma [START_REF] Usha | Modeling of stationary waves on a thin viscous film down an inclined plane at high Reynolds numbers and moderate Weber numbers using energy integral method[END_REF] and of Abderrahmane and Vatistas [START_REF] Abderrahmane | Improved two-equation model for thin layer fluid flowing down an inclined plane problem[END_REF] are not studied here because they add only third-order terms to the model of Ruyer-Quil and Manneville [START_REF] Ruyer-Quil | Improved modeling of flows down inclined planes[END_REF] that change very little the neutral stability curve (see Richard et al. [START_REF] Richard | A three-equation model for thin films down an inclined plane[END_REF]).

All models being consistent, the instability threshold [START_REF] Liu | Solitary wave dynamics of film flows[END_REF] and the values for castor oil are taken from Kliakhandler et al. [START_REF] Kliakhandler | Viscous beads on vertical fibre[END_REF].

ρ ν γ Ka kg•m -3 mm 2 •s -1 mN•m -1 - Glycerin-
at k = 0 is equal to the value calculated from the Orr-Sommerfeld equation. The curves are nearly identical to the Orr-Sommerfeld equation when k is small i.e. in the validity domain of the asymptotic method used to derive them. When k (which is of the order of ε) and Re become bigger, differences arise between models because these values are beyond the validity range of the asymptotic method. Note that the curves presented in Figure 4 extend to values of the Reynolds number which are unusually high for this kind of models (until Re = 800 for water on a vertical wall and until Re = 400 for the glycerinwater solution at θ = 6.4 • ). This allows to show clearly the differences between the models but for more usual values of the Reynolds number, the deviations are much more moderate.

Writing the diffusive terms as in (105) with α = 4/3 leads to a clear improvement of all linear properties of the models and this will be kept thereafter. The best neutral stability curves are then obtained with a MIME model and C 1 = 1 but the nonlinear properties, estimated from the velocity of a solitary wave (see §6), are better with the equivalent MIME model and even better with the corrected Shkadov model. On the other hand, as pointed out in §6, the mathematical structure and thus the reliability of the numerical resolution are better with a MIME model than with the corrected Shkadov model. It follows that, while none of these three models is clearly superior to the others on all points, the equivalent MIME model with α = 4/3 seems to be the best compromise: it gives the right solitary wave velocity on a vertical wall and acceptable linear properties and it has a good mathematical structure which allows to derive an extended model for the numerical resolution.

Numerical results

Since the MIME models admit a capillary energy, an extended MIME model including the variable W (78) can be written and the results of Noble and Vila [START_REF] Noble | Stability Theory for Difference Approximations of Euler-Korteweg Equations and Application to Thin Film Flows[END_REF] and Bresch et al. [START_REF] Bresch | A generalization of the quantum Bohm identity: Hyperbolic CFL condition for Euler-Korteweg equations[END_REF] apply. It can be proved that this numerical scheme is nonlinearly stable in the case of a Rusanov Riemann solver. An explicit scheme is nonlinearly stable under a Courant-Friedrichs-Levy (CFL) condition (Noble and Vila [START_REF] Noble | Stability Theory for Difference Approximations of Euler-Korteweg Equations and Application to Thin Film Flows[END_REF]) and a semi-implicit scheme is unconditionally nonlinearly stable (Bresch et al. [START_REF] Bresch | A generalization of the quantum Bohm identity: Hyperbolic CFL condition for Euler-Korteweg equations[END_REF]). This consideration and the previous sections show that the twoequation models which are the most suited to a numerical resolution both for the expected accuracy of the results and for the reliability of the numerical scheme are the MIME model with C 1 = 1 and the equivalent MIME model. On the other hand, while it is expected that the corrected Shkadov model gives accurate results, no extended model can be written for lack of a capillary energy and consequently no proof of the nonlinear stability of the numerical scheme can be found.

The equations of the extended MIME model are written in dimensional form in order to make comparisons with experimental results. The viscous terms are optimized by taking α = 4/3 (see §7). These equations can be written

∂h ∂t + ∂hU ∂x = 0 (120) ∂hU ∂t + ∂ ∂x hU 2 + g 2 h 5 sin 2 θ Kν 2 + C 1 2 gh 2 cos θ = C 1 gh sin θ -3ν U h + 6C 1 ν ∂ ∂x h ∂U ∂x -3C 1 ν ∂ ∂x U ∂h ∂x + 3C 1 ν U h ∂h ∂x 2 + γC 1 ρ ∂ ∂x h 3/2 ∂W ∂x (121) ∂hW ∂t + ∂hU W ∂x = - γC 1 ρ ∂ ∂x h 3/2 ∂U ∂x (122)
where the quantity W writes in dimensional form

W = γC 1 ρh ∂h ∂x (123)
From the previous sections, the parameters giving the best results are

C 1 = 1 ; K = 225 2 (124)
or

C 1 = 31 30 - √ 6 45 ; K = 315 -90 √ 6 ( 125
)
More details on the numerical scheme are given in Appendix F.

The numerical resolution of these MIME models was used to simulate the experiments of Liu and Gollub [START_REF] Liu | Solitary wave dynamics of film flows[END_REF]. The fluid used is a glycerin-water solutions whose physical properties are ν = 6.28 × 10 -6 m 2 • s -1 , ρ = 1070 kg • m -3 and γ = 6.7 × 10 -2 N • m -1 . The inclination angle is θ = 6.4 • and the Reynolds number is Re = 19.3. A small perturbation was applied at the channel's inlet with a forcing frequency of 1.5 Hz. The Nusselt flow being unstable, waves appear and grow. After a transient evolution, a periodic train of waves is reached. The final periodic solution depends only on the frequency of the inflow perturbation and not on its amplitude albeit the length of the transient evolution does depend on its amplitude.

The numerical results for these two MIME models are presented in Figure 5(a) with the experimental measures of Liu and Gollub [START_REF] Liu | Solitary wave dynamics of film flows[END_REF]. Although there is no similar numerical scheme, the corrected Shkadov model was solved in the particular case of a periodic regime with TRIFLOW, a Python solver for partial differential equations (Cellier [43]). This solver uses a second-order finite difference method, an implicit centred scheme and the method of lines. The EIM model was solved in exactly the same way for comparison purpose. The result obtained for the corrected Shkadov model is compared to the equivalent MIME model solved with the numerical scheme of Appendix F in Figure 5(b). The equivalent MIME model, the EIM model of Ruyer-Quil and Manneville [START_REF] Ruyer-Quil | Improved modeling of flows down inclined planes[END_REF], the direct numerical simulation (DNS) of Malamataris et al. [START_REF] Malamataris | Solitary waves on inclined films: Flow structure and binary interactions[END_REF], the version of the three-equation model of Richard et al. [START_REF] Richard | A three-equation model for thin films down an inclined plane[END_REF] and the experimental measures of Liu and Gollub [START_REF] Liu | Solitary wave dynamics of film flows[END_REF] are presented in Figure 5(c).

A difficulty arises because there is a marked discrepancy between the experimental measures and the DNS. The wavelength and thus the wave velocity predicted by the DNS is smaller than the experimental values and the wave amplitude calculated by the DNS is significantly greater. The amplitudes calculated by all models are smaller than the DNS value. The wave velocity of the corrected Shkadov model is the same as that of the DNS while the equivalent MIME model predicts a velocity slightly higher although the difference seems negligible. This was expected since it was shown in §6 that the corrected Shkadov model gives the same wave velocities as the DNS and that the equivalent MIME model gives the same velocity on a vertical wall but slightly greater velocities for non-vertical slopes. The even higher velocity of the EIM model and of the MIME model with C 1 = 1 agrees also with the results of §6. Compared to the experimental values, the equivalent MIME model gives a smaller wave velocity and a slightly higher wave amplitude. The version of the three-equation model proposed in Richard et al. [START_REF] Richard | A three-equation model for thin films down an inclined plane[END_REF] gives the closest agreement with the experimental measures and the highest difference with the DNS. The capillary ripples are correctly predicted for all models.

A second simulation was done in the conditions of the experiments of Dietze et al. [START_REF] Dietze | Experimental study of flow separation in laminar falling liquid films[END_REF]. A liquid film of dimethylsulphoxide (DMSO) (kinematic viscosity : ν = 2.85 × 10 -6 m 2 • s -1 , density ρ = 1098.3 kg • m -3 and surface tension γ = 48.4 mN • m -1 ) is falling on a vertical wall (θ = 90 • ). The Reynolds number is Re = 15.0 and the frequency of the inflow perturbation is 16 Hz. In this experiment, the air was subjected to an aerostatic pressure drop so that the effect of the gas on the film was negligible.

The results of the equivalent MIME model and the MIME model with C 1 = 1 are compared with the experimental measures of Dietze et al. [START_REF] Dietze | Experimental study of flow separation in laminar falling liquid films[END_REF] in Figure 6. Both models predict wave amplitudes significantly smaller than the experimental values. The overall curves are slightly shifted towards the smaller depths. The result of the DNS of Dietze et al. [START_REF] Dietze | Experimental study of flow separation in laminar falling liquid films[END_REF] (not reproduced in this paper) shows a good agreement with the experimental measures. The wave amplitude predicted by the DNS is slightly higher than the measures but the difference seems negligible.

▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ ▫ 0.

Conclusion

Many consistent two-equation models lack theoretical coherence. The EIM models do not satisfy the averaged momentum equation and are not Galilean invariant. The MIM models are Galilean-invariant but only the MIME models satisfy the averaged work-energy theorem and include a capillary energy, which is needed to write an extended model in the manner of Noble and Vila [START_REF] Noble | Stability Theory for Difference Approximations of Euler-Korteweg Equations and Application to Thin Film Flows[END_REF] and to prove the nonlinear stability of the numerical scheme. However the energy of the MIME model is not consistent to the physical energy and is thus called a mathematical energy.

The linear and nonlinear properties of consistent twoequation models were studied in order to find an optimized model. It was found that no two-equation model is better than all others in every respect. For all models, the viscous diffusive terms can be written in an optimized form by the choice (117). This value improves the neutral stability curves and the variations of the kinematic waves velocity along these curves. The three following models are clear improvements of two-equation models:

-The corrected Shkadov model is a consistent version of the well-known Shkadov model. It gives wave velocities in very good agreement with the DNS for all slopes but it lacks a capillary energy and does not satisfy the energy equation.

-The equivalent MIME model gives the same wave velocity as the corrected Shkadov model on a vertical wall but slightly faster ones on non-vertical walls. On the other hand, it satisfies the energy equation and admits a capillary energy. It is well suited to numerical resolutions since an extended model can be derived.

-The MIME model with C 1 = 1 gives neutral stability curves in excellent agreement with the Orr-Sommerfeld curve. The predicted wave velocities are a little bit too fast but, as the equivalent MIME model, it is well suited to numerical resolutions. This is the model of Lavalle et al. [START_REF] Lavalle | A numerical reduced model for thin liquid films sheared by a gas flow[END_REF] improved by the choice (117) for the viscous terms.

All these models are Galilean invariant. It is assumed that the viscous terms are written in the optimal way as described in §7. A common drawback of these models is that the predicted wave amplitudes are not always accurate. This problem cannot be solved within the scope of two-equation models. This justifies the interest of more complicated models that could improve the accuracy on the wave amplitudes. On the other hand, the simplicity of the MIME models and the existence of an extended model with the variable W (78) allow an easy and reliable numerical resolution.

Appendix A. Equations and boundary conditions in dimensionless form

Continuity equation:

∂ ũ ∂ x + ∂ w ∂ z = 0 (A.1)
Navier-Stokes momentum balance equation in the Ox-direction:

∂ ũ ∂ t + ∂ ũ2 ∂ x + ∂ ũ w ∂ z = λ εRe - 1 Fr 2 ∂ p ∂ x + ε Re ∂ τxx ∂ x + 1 εRe ∂ τxz ∂ z (A.2)
Navier-Stokes momentum balance equation in the Oz-direction:

ε 2 Fr 2 ∂ w ∂ t + ∂ ũ w ∂ x + ∂ w2 ∂ z = -cos θ - ∂ p ∂ z + εFr 2 Re ∂ τzz ∂ z + εFr 2 Re ∂ τxz ∂ x (A.3)
Constitutive relation of Newtonian fluids:

τxz = ∂ ũ ∂ z + ε 2 ∂ w ∂ x ; τxx = -τ zz = 2 ∂ ũ ∂ x (A.4)
No-penetration and no-slip conditions at the bottom:

w (z = 0) = 0 ; ũ (z = 0) = 0 (A.5)

Kinematic boundary condition at the free surface:

∂ h ∂ t + ũ h ∂ h ∂ x = w h (A.6)
Dynamic boundary conditions:

τxz h + εRe Fr 2 p h -ε 2 τxx h ∂ h ∂ x + κ εRe Fr 2 ∂ h ∂ x ∂ 2h ∂ x2   1 + ε 2 ∂ h ∂ x 2   -3 2 = 0 (A.7) p h + εFr 2 Re ∂ h ∂ x τxz h - εFr 2 Re τzz h + κ ∂ 2h ∂ x2   1 + ε 2 ∂ h ∂ x 2   -3 2 = 0 (A.8)
Work-energy theorem :

∂ ∂ t ũ2 2 + ε 2 w2 2 + ∂ ∂ x ũ ũ2 2 + ε 2 w2 2 - x sin θ εFr 2 + z cos θ Fr 2 + pũ Fr 2 - ε Re (τ xx ũ + τxz w) + ∂ ∂ z w ũ2 2 + ε 2 w2 2 - x sin θ εFr 2 + z cos θ Fr 2 + p w Fr 2 - 1 εRe τxz ũ - ε Re τzz w = - 4ε Re ∂ ũ ∂ x 2 - 1 εRe ∂ ũ ∂ z 2 - 2ε Re ∂ ũ ∂ z ∂ w ∂ x - ε 3 Re ∂ w ∂ x 2 (A.9)

Appendix B. Conditions of compatibility between the momentum and energy equations for two-equation models

The energy equation ( 19) can be written

U ∂ h Ũ ∂ t + ∂ ∂ x h Ũ 2 + α 3 h2 cos θ 2Fr 2 + ∂E 2 ∂ t + ∂Fr 2 ∂ x + (α 2 -α 3 ) ∂ ∂ x h2 U cos θ Fr 2 = C 2 Ũ εRe λ h - 3 Ũ h (B.1)
On the other hand, the momentum equation ( 18) is written

∂ h Ũ ∂ t + ∂ ∂ x h Ũ 2 + α 3 h2 cos θ 2Fr 2 = C 1 εRe λ h - 3 Ũ h - ∂F 1 ∂ x + (α 3 -α 1 ) ∂ ∂ x h2 cos θ 2Fr 2 (B.2)
A compatibility condition between the momentum and energy equation is thus

∂E 2 ∂ t + ∂F 2 ∂ x - Ũ ∂F 1 ∂ x + C 2 -α 3 4 h2 cos θ Fr 2 ∂ Ũ ∂ x + C 2 + α 3 -2C 1 2 h Ũ cos θ Fr 2 ∂ h ∂ x = C 2 -C 1 εRe Ũ λ h - 3 Ũ h (B.3)
Given that E 2 , F 2 and F 1 are functions of the two variables h and Ũ , this condition writes

∂E 2 ∂ h ∂ h ∂ t + ∂E 2 ∂ Ũ ∂ Ũ ∂ t + ∂F 2 ∂ h ∂ h ∂ x + ∂F 2 ∂ Ũ ∂ Ũ ∂ x - Ũ ∂F 1 ∂ h ∂ h ∂ x - Ũ ∂F 1 ∂ Ũ ∂ Ũ ∂ x + C 2 -α 3 4 h2 cos θ Fr 2 ∂ Ũ ∂ x + C 2 + α 3 -2C 1 2 h Ũ cos θ Fr 2 ∂ h ∂x 20 = C 2 -C 1 εRe Ũ λ h - 3 Ũ h (B.4)
The derivative of the depth with respect to the time is given by the averaged mass equation ( 9)

∂ h ∂ t = - h ∂ Ũ ∂ h - Ũ ∂ h ∂ x (B.5)
The derivative of the average velocity with respect to the time is given by the averaged momentum equation ( 18)

∂ Ũ ∂ t = - Ũ ∂ Ũ ∂ x - 1 h ∂F 1 ∂ h ∂ h ∂ x - 1 h ∂F 1 ∂ Ũ ∂ Ũ ∂ x -C 1 cos θ Fr 2 ∂ h ∂ x + 1 εRe C 1 h λ h - 3 Ũ h (B.6)
The compatibility condition can then be written

∂ h ∂ x - Ũ ∂E 2 ∂ h - 1 h ∂E 2 ∂ Ũ ∂F 1 ∂ h -C 1 cos θ Fr 2 ∂E 2 ∂ Ũ + ∂F 2 ∂ h - Ũ ∂F 1 ∂ h + C 2 + α 3 -2C 1 2 h Ũ cos θ Fr 2 + ∂ Ũ ∂ x - h ∂E 2 ∂ h - Ũ ∂E 2 ∂ Ũ - 1 h ∂F 1 ∂ Ũ ∂E 2 ∂ Ũ + ∂F 2 ∂ Ũ - Ũ ∂F 1 ∂ Ũ + C 2 -α 3 4 h2 cos θ Fr 2 + ∂E 2 ∂ Ũ C 1 εRe h λ h - 3 Ũ h = C 2 -C 1 εRe Ũ λ h - 3 Ũ h (B.7)
This implies firstly that

∂E 2 ∂ Ũ = C 2 -C 1 C 1 h Ũ (B.8)
and then the two following relations

∂F 2 ∂ h = Ũ ∂E 2 ∂ h + Ũ C 2 C 1 ∂F 1 ∂ h + C 2 -α 3 2 h Ũ cos θ Fr 2 (B.9) ∂F 2 ∂ Ũ = h ∂E 2 ∂ h + C 2 -C 1 C 1 h Ũ 2 + Ũ C 2 C 1 ∂F 1 ∂ Ũ - C 2 -α 3 4 h2 cos θ Fr 2 (B.10) Writing ∂ ∂ Ũ ∂F 2 ∂ h = ∂ ∂ h ∂F 2 ∂ Ũ (B.11) leads to - h ∂ 2 E 2 ∂ h2 + C 2 C 1 ∂F 1 ∂ h = (α 3 -C 2 ) h cos θ Fr 2 (B.12)
Since E 2 and F 1 do not depend on θ, this relation implies that

C 2 = α 3 . (B.13)
Then equation [START_REF] Benjamin | Wave formation in laminar flow down an inclined plane[END_REF] gives α 2 = C 2 . Consequently we obtain the condition

h ∂ 2 E 2 ∂ h2 = C 2 C 1 ∂F 1 ∂ h (B.14)

Appendix C. Consistency to the physical energy

We want to derive all consistent models compatible with the momentum and energy equations and with an energy consistent to the physical energy. The energy e (40) is consistent to

2λ 2h4 30 + h cos θ 2Fr 2 + O (ε) (C.1)
The consistency of the term in cos θ implies that α 3 = 1. The compatibility condition (B.13) gives then C 2 = 1 and the consistency condition (39) yields α 2 = 1. The integration of the relation (B.8) gives

E 2 = 1 C 1 -1 h Ũ 2 2 + K 1 h (C.2)
where K 1 ( h) is a function of the one variable h. The consistency of the physical average energy (C.1) implies that E 2 is consistent to λ 2h5 /90. We can thus calculate the expression to whom K 1 is consistent. Since K 1 depends only on h, there is only one solution, which is

K 1 h = 6C 1 -5 90C 1 λ 2 h5 (C.3)
From equation (B.14), we find

∂F 1 ∂ h = 2 9 (6C 1 -5) λ 2 h4 (C.4)
and then, by integration

F 1 h, Ũ = 2 45 (6C 1 -5) λ 2 h5 + K 2 Ũ (C.5) where K 2 ( Ũ ) is a function of the one variable Ũ . Calcu- lating that F 1 is consistent to 2 25 λ 2 h5 10 9 -C 1 (C.6)
and since K 2 depends only on Ũ , there is only one solution

K 2 Ũ = 2 5 3 λ Ũ 5/2 7 - 39 5 C 1 (C.7)
Finally, F 2 can be found with the relations (B.9) and (B.10):

F 2 h, Ũ = 2 C 1 3 λ 1 - 39 35 C 1 Ũ 7/2 + 1 2 h Ũ 3 1 C 1 -1 + λ 2h5 Ũ 3 1 - 5 6C 1 (C.8)
The terms in Ũ

Appendix D. Mathematical entropy and capillary energy

The conservative hyperbolic part of the model admits a mathematical entropy if there is a convex function η of the conservative variables h and q that satisfies the additional conservation law

∂η ∂t + ∂Φ ∂x = 0 (D.1)
where Φ is a function of h and q. This conservation law can be written The characteristic velocities of this system are λ = AU ± (A 2 -A)U 2 + a 2 with a 2 = dP/dh. The system is unconditionally hyperbolic if A 1. This condition is thus supposed to be satisfied in the following. Equation The function η is a convex function of h and q if and only if -1 m 0. If A 1, the root m + with the plus sign in the above expression always satisfies this condition whereas the root m -with the minus sign never satisfies it. The model thus admit a mathematical entropy with the exponent m + . For a model with A = 6/5 and a 2 = 5h cos θ/(6Fr 2 ), we get m + = -4/5 and the expression of the mathematical entropy is If additionally C 1 = 1 the capillary terms (73) of the averaged work-energy theorem [START_REF] Chakraborty | Extreme solitary waves on falling liquid films[END_REF] are recovered. This is the case of the MIME model with C 1 = 1 and F 1 = 2λ 2 h 5 /225 (which implies that A = 1 and m = -1) used by Lavalle et al. [START_REF] Lavalle | A numerical reduced model for thin liquid films sheared by a gas flow[END_REF].

η = q 2 h 4/

Appendix E. Hyperbolicity

Taking κ = 0, at the first order of approximation there is no dispersive and no diffusive terms. The model should then be unconditionnally hyperbolic. The mass [START_REF] Malamataris | Solitary waves on inclined films: Flow structure and binary interactions[END_REF] and momentum [START_REF] Noble | Stability Theory for Difference Approximations of Euler-Korteweg Equations and Application to Thin Film Flows[END_REF] equations of MIM models can be written in dimensional form

∂U ∂t + A 1 ∂U ∂x = S 1 (E.1)
where U = [h, U ] T , S 1 = [0, C 1 (g sin θ -3νU/h 2 )] T and

A 1 = U h 1 h ∂F 1 ∂h + C 1 g cos θ U + 1 h ∂F 1 ∂U (E.2)
Here F 1 is written in dimensional form. The characteristics of system (E.1) are

U + 1 2h ∂F 1 ∂U ± 1 4h 2 ∂F 1 ∂U 2 + ∂F 1 ∂h + C 1 gh cos θ (E.
3)

The system (E.1) is unconditionally hyperbolic if

∂F 1 ∂h 0 (E.4)
Both the MIM models ( 20) and ( 21) fulfil this condition.

Appendix F. Numerical scheme

An explicit numerical scheme is used. It is nonlinearly stable under CFL condition. The hyperbolic part of the equations is solved by a classical finite volume method. If X is any one of the variables h, U , q = hU , W , hW , qW or if it is the flux

F = hU 2 + g 2 h 5 sin 2 θ Kν 2 + C 1 2 gh 2 cos θ (F.1)
then the value of X at cell i and at iteration n is denoted by X n i . The intercell fluxes are denoted by X n i+1/2 and X n i-1/2 . They are calculated by the Rusanov Riemann solver. Note that for the hyperbolic system where c is the shock velocity, m is h(U -c) which is constant through the shock and where [[X]] denotes the jump through the shock of the quantity X i.e. the difference X 2 -X 1 between the values of X on both sides of the shock. Since the third equation is a transport equation for W which is totally separated from the two other equations, the shock wave amplitude has no theoretical upper bound.

Of course because of the diffusive terms, there will be no shock in the whole system. The viscous terms are handled by a finite difference method. The notations X n i+1/2 and X n i-1/2 refer to the average values X n i+1 + X n i /2 and X n i + X n i-1 /2 respectively. The time step is denoted by ∆t and the cell width by ∆x. The numerical scheme can be written

h n+1 i = h n i + ∆t ∆x q n i-1/2 -q n i+1/2
(F.6)

q n+1 i = q n i + ∆t ∆x F n i-1/2 -F n i+1/2
+ ∆t C 1 gh n i sin θ -3ν

U n i h n i +6C 1 ν ∆t (∆x) 2 h n i+1/2 U n i+1 -U n i -h n i-1/2 U n i -U n i-1 -3C 1 ν ∆t (∆x) 2 U n i+1/2 h n i+1 h n i -U n i-1/2 h n i -h n i-1 + 3C

Figure 1 :

 1 Figure 1: Notations used in the text.

  )), then the characteristic velocities are of the form U ± a. The critical point corresponds to the point where m = -ha or c = U + a where a = ∂F 1 ∂h + C 1 gh cos θ (84) is the velocity of the surface waves. The critical point splits the wave into a region where the flow is subcritical (c < U + a) and a region where it is supercritical (c > U + a). If ∂F 1 /∂U = 0, the characteristic velocities (E.3) are of the form U + a ± a and the critical point corresponds to the point where m = -h(a + a ) (or c = U + a + a ) with a = 1 2h

ln c 3 2 Figure 2 :Figure 3 :

 223 Figure 2: (a) Velocity of the kinematic waves as a function of the wavenumber along the neutral stability curve for the Orr-Sommerfeld equations. Thick solid curve: water θ = 6.4 • ; thin dashed curve: alcohol θ = 6.4 • ; thick dashed curve: alcohol θ = 90 • ; dotted curve : castor oil θ = 90 • ; thin solid curve: castor oil θ = 6.4 • . (b) Logarithmic plot of c -3/2 as a function of the wavenumber along the neutral stability curve. Thick solid curve: castor oil θ = 90 • : thin solid curve; castor oil θ 6.4 • ; The dashed line and the dotted line have a slope of -2.

Figure 4 :

 4 Figure 4: Neutral stability curves. Dashed curve: Orr-Sommerfeld. (a), (c) and (d) water for θ = 90 • ; (b) glycerin-water solution for θ = 6.4 • ; (a) and (b) equivalent MIME model with α = 1 (thin solid curve), α = 6/5 (dotted curve), α = 4/3 (thick solid curve). (c) MIME model with α = 4/3 and C 1 = 1 (thin solid curve), C 1 = 25/24 (dotted curve) and for the equivalent MIME model (thick solid curve). (d) equivalent MIME model with α = 4/3 (thick solid curve); corrected Shkadov model with α = 4/3 (thin solid curve); EIM model of Ruyer-Quil and Manneville [5]: dotted curve; EIM model of Mudunuri and Balakotaiah [12]: dot-dashed curve.

5 :

 5 Depth profiles in the conditions of the experiments of Liu and Gollub [20]. The abscissa origins are arbitrary. (a) Thick solid curve: equivalent MIME model; thin solid curve: MIME model with C 1 = 1; dashed curve: experimental measures of Liu and Gollub [20]. (b) Thick curve: equivalent MIME model; thin curve: corrected Shkadov model. (c) Thick solid curve: equivalent MIME model; thin curve: DNS; thin dotted curve: EIM model of Ruyer-Quil and Manneville [5]; thick dashed curve: experimental measures of Liu and Gollub [20]; thin dashed curve: three-equation model of Richard et al. [1].

Figure 6 :

 6 Figure 6: Depth profile in the conditions of the experiments of Dietze et al. [21]. Thick curve: equivalent MIME model; thin curve: MIME model with C 1 = 1; squares: experimental measures of Dietze et al. [21].

Fr 2 h m q ∂h ∂x 2 + 2 ∂h ∂x 2 (

 222 can be fulfilled only if n = 2. Then r(h) = C 1 κh m+1 /Fr 2 and the second condition yieldsΨ = -(m + 1)C 1 κ Υ(h x , h) (D.22)where Υ is an unknown function of h x and h. be satisfied only if m = -1 i.e. if A = 1.There is no solution in the case A = 6/5. If A = 1, we get finally r(h) = C 1 κ/Fr 2 ,η c = C 1 κ Fr

2 +

 2 g 2 h 5 sin 2 θ Kν 2 + C1 2 gh 2 cos θ hU W   = 0 (F.2) the characteristics are U ; U ± 5g 2 h 4 sin 2 θ Kν 2 + C 1 gh cos θ (F.3) and the shock relations would be [[h (U -c)]] = 0 ; [[W ]

  1) (if second-order terms are included in the model then κ = O(ε) to avoid unnecessary complicated capillary terms). The Reynolds number, the Froude number and the angle θ are assumed to be of O(1). The components of the viscous stress tensor are defined as τ = τ xx e

x ⊗ e x + τ xz e x ⊗ e z + τ xz e z ⊗ e x + τ zz e z ⊗ e z , e x and e z being the unit vectors in the Ox and Oz directions. The dimensionless quantities are denoted with a tilde and defined as follows:

Table 1 .

 1 10 and F 3 = 6hU 3 /35. The denominator vanishes at the critical point for m = -h(a + a ) where: The MIM model with C 1 = 5/6 and F 1 = hU 2 /5 gives the right values of the velocity of the solitary wave for all

	a =	3 14	U ; a =	37 196	U 2 +	5 6	gh cos θ	(96)

The same procedure as for the MIM models is then followed to calculate c/(3U ∞ ). The results are presented in

Table 1 :

 1

Values of c/ (3U∞) for five values of F r -2 ranging from 0 (vertical wall) to 0.4 (instability threshold) for various types of models [Momentum integral method (MIM), momentum integral with mathematical energy (MIME), momentum-energy compatible models (Compatible) (see Appendix C), energy integral method (EIM), four-equation model of Ruyer-Quil and Manneville

[START_REF] Ruyer-Quil | Improved modeling of flows down inclined planes[END_REF]

(RQM4)] and for a DNS. The RQM4 and DNS values are taken from Chakraborty et al.

[START_REF] Chakraborty | Extreme solitary waves on falling liquid films[END_REF]

. For the MIM, MIME and Compatible models, the values of C 1 and F 1 are given.

Table 2 :

 2 Limit values of the kinematic waves velocity when k → ∞, c lim , and of this velocity and of the wavenumber when Re → ∞, c∞ and k∞ respectively (see text). Corrected Shkadov model or equivalent MIME model: CS/EqMIME

Table 3 :

 3 Fluid properties (density, kinematic viscosity, surface tension and Kapitza number) used in the present study. The values for the glycerin-water solution are taken from Liu and Gollub

	water	1134	6.28	67.0	238
	Water (20 • C)	998	1.007	72.8	3375
	Castor oil	961	440	31.0	0.45

  5/2 and Ũ 7/2 can be cancelled by taking C 1 = 35/39 but this model does not satisfy the Galilean invariance since C 1 = C 2 . The only one Galilean-invariant model has C 1 = 1 but retains the undesirable terms in Ũ 5/2 and Ũ 7/2 .

  It can be proved that η depends on q only asη = h m q 2 + η p (h) (D.8)where the function η p depends only on h and m is any exponent. With this expression, the compatibility condition givesd 2 η p dh 2 = 2a 2 h m

											The compatibility condition yields
											a 2 -A	q 2 h 2	∂ 2 η ∂q 2 =	∂ 2 η ∂h 2 + 2A	q h	∂ 2 η ∂h∂q	(D.7)
											(D.9)
											and
											m =	1 2	1 -4A ± (1 -4A) 2 -8A	(D.10)
											(D.2)
	becomes
	∂h ∂x		A	q 2 h 2 -a 2 ∂η ∂q	+	∂Φ ∂h
				+	∂q ∂x		-	∂η ∂h	-2A	q h	∂η ∂q	+	∂Φ ∂q	= 0 (D.5)
	This implies that
	        	∂Φ ∂h ∂Φ ∂q	q 2 h 2 h + 2A = a 2 -A = ∂η ∂h q	∂q ∂η	∂η ∂q	(D.6)

  The system admits a capillary energy if there is a function η c depending on h and ∂h/∂x such that, by including capillarity, the additional conservation law becomes is a function of h, q, ∂h/∂x, ∂q/∂x and ∂ 2 h/∂x 2 . The second equation of the model is now , q x and h xx denote ∂h/∂x, ∂q/∂x and ∂ 2 h/∂x 2 respectively. This implies that , h, q) (D.[START_REF] Bresch | A generalization of the quantum Bohm identity: Hyperbolic CFL condition for Euler-Korteweg equations[END_REF] where Ψ is another unknown function depending on h x , h and q only. The condition (D.16) becomes Since Ψ does not depend on h x nor on h xx , we get the conditions

	where h x ∂Φ c ∂h xx = -2C 1	κ Fr 2 h m+1 q				(D.17)
	and
	∂Φ c ∂q x	= nr(h)	∂h ∂x	n-1				(D.18)
	Consequently the expression of Φ c can be written
	Φ c = -2C 1	κ Fr 2 h m+1 q	∂ 2 h ∂x 2
	+ nr(h) + Ψ(h x ∂q ∂h ∂x n-1 ∂q ∂x ∂x (n -1) dr dh ∂h ∂x n ∂x 2 + n(n -1)r(h) ∂ 2 h	∂x ∂h	n-2
	-2C 1	κ Fr 2 h m+1 ∂ 2 h ∂x 2 +	∂Ψ ∂q	= -	∂Ψ ∂h x	∂ 2 h ∂x 2
				+ 2C 1	κ Fr 2 (m + 1)h m q	∂ 2 h ∂x 2	∂h ∂x	-	∂Ψ ∂h	∂h ∂x	(D.20)	5 +	125 198	h 11/5 cos θ Fr 2	(D.11)
	                  	n(n -1)r(h) ∂Ψ ∂q + (n -1) ∂Ψ ∂h x ∂ 2 h ∂x 2 + ∂Ψ dr ∂h ∂x dh ∂h ∂h ∂x	n-2 ∂h ∂x = 2C 1 = 2C 1 n = 0 κ Fr 2 (m + 1)h m q κ Fr 2 h m+1	∂ 2 h ∂x 2	∂ ∂t where Φ c ∂q (η + η c ) + ∂t + ∂ ∂x A q 2 ∂x ∂ h + a 2 (h) = C 1 (Φ + Φ c ) = 0 Fr 2 h κ	∂ 3 h ∂h 3	(D.12) (D.13)
														It follows that η c and Φ c satisfy
														∂η c ∂t	+	∂Φ c ∂x	+ C 1	κ Fr 2 h	∂ 3 h ∂h 3	∂η ∂q	= 0	(D.14)
														Taking
														η c = r(h)	∂h ∂x	n	(D.15)
														where n is an exponent and r a function of h only, we find
														the following condition
														dr dh	∂q ∂x	∂h ∂x	n + nr(h)	∂ 2 q ∂x 2	∂h ∂x	n-1 -2C 1	κ Fr 2 h m+1 q	∂ 3 h ∂x 3
														=	∂Φ c ∂h x	∂ 2 h ∂x 2 +	∂Φ c ∂q x	∂ 2 q ∂x 2 +	∂Φ c ∂h	∂h ∂x	+	∂Φ c ∂q	∂q ∂x	+	∂Φ c ∂h xx	∂ 3 h ∂x 3
														(D.16)

Acknowledgments

The work was supported by the Service Hydrographique et Océanographique de la Marine (SHOM) and by the Université Savoie Mont Blanc (project ECOPP).