DERIVATION VIA HAMILTON’S PRINCIPLE OF A NEW SHALLOW-WATER MODEL USING A COLOR FUNCTION FOR THE MACROSCOPIC DESCRIPTION OF PARTIAL WETTING PHENOMENA

. This paper presents a new shallow-water type model suitable for the simulation of partially wetting liquid films without the need for very fine resolution of the contact line phenomena, which is particularly suitable for industrial applications. This model is based on the introduction of a color function, propagated at the averaged velocity of the bulk flow, and equal to one where there is a liquid film and zero in the dry zone, which implies a non zero gradient only at the interface. This approach has the advantage of easily locating the interface, allowing to model macroscopically the forces acting at the contact line, which is essential for the simulation of partial wetting phenomena. The formal derivation of this model is based on the principle of least action known as Hamilton’s principle. Here this principle is applied in a full Eulerian form to derive the complete system of equations with the color function. This method proves to be particularly suitable for this type of development and as an illustration it is also applied to recover another model proposed by Lallement et al. [39,73]. Finally, both models are compared from a theoretical point of view and the advantages of the new color function based model are discussed.


Introduction
Situations where one or more fluids meet and interact with a solid surface are very common in our natural or industrial environment [4,12]. Whether it's a drop on the leaves in the morning dew, water flowing from the wings of a duck coming out of the water [52], during a wave soldering operation [4] or in the field of inkjet printing [14,69]. All these situations belong to the physical phenomenon called wetting / dewetting. Despite the multiplication of experimental, theoretical and numerical studies since the middle of the 20 th century, wetting and even more dewetting is still poorly understood.
From a macroscopic point of view, the wetting is described with the interfacial surface tensions. Thus, in the case of a problem with three phases, three surface tensions are defined: γ lg , γ sg and γ sl which correspond to the surface energy per unit area of the interface between the liquid film and the air, the solid wall and the air, the liquid film and the solid wall, respectively. The relation between these three quantities is fundamental to characterize the different types of wetting. For that it is necessary to introduce the spreading parameter: If S > 0 the wetting will be total i.e. the liquid will completely cover the solid. Whereas in the limit S ≤ 0 the wetting will be partial i.e. the liquid will form a contact line with the solid surface and meet this surface with a finite static contact angle θ s . The law characterizing this equilibrium in partial wetting is the well-known Young-Dupré's law illustrated in Figure 1: The surface free energy density of this macroscopic film is then: 1 ONERA/DMPE, Universite de Toulouse F-31055 -France Figure 1. Balance of capillary forces close to the contact line: Young-Dupré law.
By minimizing this energy with a constant volume constraint, the liquid thickness at the center of the liquid film is equal to [13]: If static equilibrium is well understood, the dynamic of wetting phenomena is a more complex subject. This is largely due to the fact that classical hydrodynamic theory based on a continuous approach breaks down when describing a moving contact line. Indeed, this approach leads to the appearance of the well known non-integrable singularity in the wall shear stress and pressure fields due to the application of the classical non-slip boundary condition at the wall [11,31,47]. This singularity is actually the consequence of the strong assumption of continuity of the hydrodynamic theory and highlights the molecular nature of wetting phenomena close to the contact line. Indeed, wetting is closely linked to the surface tension property of a fluid, which intrinsically comes from molecular interactions at the multiphase interfaces. However, using molecular approach, although more physical, is extremely limited especially in the perspective of making large scale simulations. This is why many modeling artefacts have been imagined to circumvent the difficulty of the singularity and to continue using a continuous approach. One of the easiest approach is to relax the non-slip constraint by introducing a slip length at the wall [19,30,48]. Another method consists in relying on the existence of a microscopic precursor film ahead of the macroscopic bulk flow involving the disappearance of a real contact line and preventing the film from advancing on a dry surface. However, this precursor film has a physical validity only in the case of total wetting [3,21,55,60]. Finally, another widely used approach is to relax the sharp interface constraint and consider physical quantities, such as the density with a smooth variation in an artificially thickened interface called diffuse interface. With this approach, it has been shown that there is no need for slip or precursor film to get rid of the singularity [32,65]. In the end, the idea remains the same to introduce a critical thickness (slip length, thickness of the precursor film or of the diffuse interface) below which we know that the hydrodynamic theory collapses. However, it is precisely at these cut scales that the physical phenomena, which explain the partial wetting properties and dewetting behaviour of a fluid, come into play. These scales not being solved by classical hydrodynamic models, different strategies had to be devised to capture the dynamics and the equilibrium state of a partial wetting film. In this study we will focus exclusively on the models in the lubrication approximation [2,51]. Indeed, this approximation is perfectly justified in view of the targeted applications where the film thickness is negligible compared to the characteristic spreading dimension. Among the different available models, in the literature, the majority are not adapted to industrial applications as we will see below.
The main objective of the present work is to develop a new robust and more versatile shallow-water type model for the simulation of partial wetting while being suitable for large-scale simulations such as the formation of rivulets on aircraft wings [76].
The present article is divided into five parts. First, we review existing models for partial wetting applications in the lubrication approximation and highlight our motivation for introducing a color function in a shallow-water type model. Secondly, we present a methodology based on Hamilton's principle in Eulerian form to recover a shallowwater type model that was previously used and verified in previous works [38,73]. Then, with this methodology we establish a new shallow-water type model using a color function approach to distinguish between dry and wet areas. Finally, we prove the consistency of this new model with respect to the description of wetting phenomena at the macroscopic scale and we end with a comparison with another existing physical model.

Existing models for partial wetting applications in the lubrication approximation
In this section we present existing models in the literature to simulate phenomena of partial wetting. We have classified these different models into two categories. A first one corresponding to the models based on the film thickness lubrication equation incorporating a disjoining pressure term to simulate the partial wetting. The second category corresponds to models based on the shallow-water equations and incorporating a macroscopic contact line force. The purpose of this section is to lay the foundation for the construction of our new model with color function that will be introduced at the end of this section.

Role of molecular interactions in dewetting phenomena
Before presenting the different existing models, it is important to briefly clarify the physical basis of the dewetting phenomenon and the link with molecular interactions. This will allow to understand the construction of the existing models and to justify the physical and theoretical foundations of our new model with color function, which will be detailed later in Section 6.2.
The dewetting can be the consequence of many phenomena such as evaporation [71], surface tension gradient due to spatial changes of temperature [13] or simply the consequence of microscopic molecular forces that break bonds [36,46,63]. The common point between all these phenomena is their molecular origin. Indeed, molecular microscopic interactions generate the destabilization of microscopic liquid films which leads to partial wetting configuration or spontaneous dewetting [7,56]. Thus, the macroscopic energetic formulation Eq.3 does not allow to describe dewetting phenomena because it does not take into account molecular interactions. It is then necessary to introduce a surface interaction potential energy e d which takes into account these effects. In the limit h → 0 this potential energy is equal to the spreading parameter S (Eq.1) to obtain a film energy equal to the surface energy of a dry solid γ sg . And, in the limit h → ∞, this energy is negligible and equal to 0. The introduction of this interaction energy allows to define a free surface energy valid for all film thicknesses: The form of the potential energy e d is fundamental to understand theoretically the phenomena of partial wetting and dewetting. Indeed, the sign of this interaction energy will characterize the attractive or repulsive character of molecular interactions in the liquid film. This will play a fundamental role in the stability of thin films. Thus, characterizing the wettability of a liquid film with only the sign of the spreading parameter S (Eq.1) may appear as a strong approximation. In reality, for microscopic films it also depends on the sign of the interaction potential energy. Four configurations are effectively possible with respect to the sign of e d and S. These different configurations have been widely studied in [8,12,34]. However, if we are only interested in films whose mean thickness is large compared to the inter-molecular force length scale (typically if h > 1µm) only the parameter S plays a role on their wettability. It is therefore clear that, depending on the scale of the phenomena to be treated, it will be possible or not to simplify the expression of the intermolecular forces (see Section 2.3). We represent in Figure 2 the evolution of this free energy for a partially wetting film of various thicknesses. For large thicknesses the gravity energy largely dominates. Its quadratic form leads to a total free energy with a convex form which is characteristic of a stable film. This film will therefore have a unique equilibrium state corresponding to a thickness h c depending on the spreading parameter S (Eq.4) and therefore will not tend to split. For very thin films (h < h i ≈ µm [56]) the energy is dominated by the interaction potential energy which has no choice but to be concave in the case of partial wetting (e d tends to S < 0 within the limit of a zero thickness and to zero for large thicknesses). This concave energy is characteristic of an unstable film. Thus, the film can change its equilibrium state naturally. A characteristic phenomenon of these very thin films is the spinodal dewetting where the film splits into several drops arranged orderly [57,58,75]. Finally, films of intermediate thickness are said to be in a meta-stable state. They will not tend to split in a natural way but with the help of an external force it is possible to force dewetting. This situation is illustrated by the arrows in Figure 2 where we can see that the total free energy of a meta-stable film can be lowered to the dotted line. This dotted line is actually the common tangent to the free energy curve for two different thicknesses: a thickness equal to zero and a thickness equal to h c . This common tangent shows the possibility of a stable configuration where a film of thickness h c coexists with a dry surface [13]. Thus, an initially meta-stable film can become more stable by splitting. This kind of dewetting is called dewetting by nucleation [56,67,74].
This figure illustrates that the molecular interactions play an essential role in the physics of dewetting. In order to model such phenomena, it is thus necessary to use a force that takes into account these molecular interactions in Figure 2. Representation of the different stability configurations of a liquid film in the case of partial wetting. h c represents the characteristic thickness of a macroscopic film at equilibrium and h i represents the characteristic thickness at which molecular effects (modeled by e disj ) predominate. The tangent to the film free energy curve common in two different points, represented in dotted line, illustrates the possible coexistence between two films of different thicknesses [13]. the models. A natural idea is to use a force derived from the theoretical disjoining pressure. Indeed, this pressure introduced by Frumkin and Derjaguin in the middle of the last century [15,16,22], is a pressure that takes into account the interaction between two surfaces (eg solid/liquid and liquid/gas) close enough to interact with molecular long range forces. Its typical expression is of the following form [9,17]: The right hand terms correspond to the pressures associated to the Van der Waals, electrostatic and structural forces, respectively. With A the Hamaker constant, λ s the characteristic length of structural force application, ϵ 0 the dielectric permittivity of the liquid, Φ 1 and Φ 2 the electrostatic potentials associated to each surface, K a factor that reflects the repulsive (K > 0) or attractive (K < 0) property of the structural force. Moreover, in the particular case of a microscopic liquid film of constant non-zero thickness the interaction potential energy e d is directly related to the disjoining pressure by [8,70]: The disjoining pressure has the advantage of adding non-negligible physical ingredients in the modeling of wetting phenomena. However, in its theoretical form (Eq.6) this pressure diverges in the limit where the film thickness becomes zero. This is particularly problematic for dewetting situations where the film, when removed, leaves dry areas. As we will see below, its use for the modeling of dewetting or partially wetting films therefore requires some precautions.
In summary, we note that molecular interactions, whose associated energy is a strictly concave function of the film thickness, play an important role in the process of partial wetting and dewetting. Thus, it is necessary to include terms in the models, such as disjoining pressure, to capture the physics related to molecular interactions to simulate partial wetting and dewetting phenomena.

Presentation of mesoscopic existing models with disjoining pressure
In the literature, the vast majority of models with disjoining pressure are derived in the lubrication approximation [2,51], and are of the form: with ρg z h the hydrostatic pressure, P g the ambient gas pressure, µ the dynamic viscosity of the liquid and F s the external force source term including the effect of the tangential gravity acceleration g t = (g x g y ) T projected in the (x-y) flow plane, and the viscous shear force of the ambient air. The form of the disjoining pressure can be completely different from one model to another [33]. For example, a very interesting approach was developed by Pismen and Pomeau [54,72] who derived a disjoining potential energy by coupling a hydrodynamic approach with a diffuse interface model to study dewetting phenomena. The expression of the disjoining pressure derived from this potential energy presents two very satisfactory advantages. First, relaxing the classical approach of a sharp liquid-gas interface allows to define a disjoining pressure which is no longer divergent when the thickness tends towards zeros. In the same time, as said in the introduction, thanks to this diffuse interface the singularity of the viscous stress is also avoided due to the fact that the location of the contact line becomes indefinite. In many other cases, the main ingredient for avoiding singularities is to consider a precursor film. To introduce this precursor film in the model, it is common to use a form of disjoining pressure which takes into account both the attractive and the repulsive forces, similarly to Eq.6. This microscopic film covers the entire surface of the solid around the macroscopic bulk flow and its thickness is given by the balance between attractive and repulsive components of the disjoining pressure [24]. The motivation of such a model is perfectly understandable but can be discussed in so far as the precursor film is only valid in the case of total wetting. However, this approach has been proven to work for modelling dewetting phenomena as shown by different studies [1,35,61,62].
These models based on a physical representation of molecular interactions, with the disjoining pressure, allow to physically represent configurations of partial wetting and dewetting but remains however not suitable for large scale simulations using industrial codes based on the finite volume method. This is particularly due to the fact that the lubrication equation (Eq.8a) is a fourth order equation. Moreover, using a very thin precursor film or using a diffuse interface requires very fine resolutions to properly resolve the fluid front, which can be very expensive and thus once again not adapted for the simulation of large scale wetting phenomena [37,43,68].

Presentation of macroscopic existing models with contact line force
For macroscopic films, using the disjoining pressure, which describes microscopic molecular interactions, is questionable. Indeed, different studies have shown that it is not necessary to model all the physics at the microscopic scale to simulate macroscopic partially wetting or dewetting films. The purpose of this section is therefore to present existing shallow-water type models with macroscopic contact line force that are particularly adapted for large scale simulations.
We consider an incompressible liquid film of thickness h(x, t) in the z-direction that spreads in the (x-y) perpendicular plane. The characteristic thickness of the film h 0 is considered negligible against the characteristic spreading dimension L 0 : h 0 /L 0 ≪ 1. Following this long wave assumption, we suppose that the dynamics of the liquid film can be described with shallow-water equations. These equations are obtained by integration of the classical Navier-Stokes equations in the z-direction and by assuming the mechanical equilibrium of the film in this direction. Contrary to the case of a single lubrication equation such as Eq.8 it implies to use an additional momentum equation. This change has the great advantage, for spatial discretization purpose, to lower the capillary force term order from four to three. However, shallow-water equations in their classical form do not allow the simulation of partially wetting films. For this purpose, one idea is to to add a contact line force F cl , concentrated at the contact line, in the momentum equation. Thus, the shallow-water equations for partial wetting applications takes the following form: with γ lg κ(x, y) the Laplace pressure related to the curvature (Eq.8c) of the interface, u the averaged velocity of the bulk flow, u ⊗ u the tensor product of velocity, averaged over the thickness of the film, τ a the air friction, τ w the wall friction, n cl is the unit vector perpendicular to the contact line oriented outside the liquid and δ cl is the Dirac measure concentrated at the contact line position. For the closure of the model, it is necessary to provide an expression for the averaged tensor product and the wall shear stress. For this the following simplification is made (which is not limiting in the lubrication regime as inertial forces play a negligible role as justified in [73]): As far as the computation of the wall shear stress is concerned, the following 3D Poiseuille velocity profile is assumed, which can be theoretically justified in the lubrication regime and assuming a no-slip boundary condition at the wall.
Using Eq.11a and Eq.11b the expression for the wall shear stress is given by: It is now necessary to choose a formulation for the contact line force F cl . Meredith et al. [44,45] used this approach in a finite volume framework. The reinterpretation of their method at continuous level leads to the following expression for the tangential force: Using the Young Dupré's law to introduce the static contact angle θ s , a physical parameter easier to estimate experimentally than the surface tensions, the contact line force takes the following form: At discrete level, in the mesh cell containing the contact line, Meredith et al. used the following approximation : with ∆x the mesh size. Hence, the transformation from a lineic force to a surface force is made by simply approximating the Dirac distribution by the function f (x) = 1/∆x in the cell containing the contact line and f (x) = 0 elsewhere. Moreover, the application of the force only close to the contact line is ensured by the use of a critical thickness criterion which allows to distinguish the dry and wet areas. As shown by the authors, this approach allows to correctly simulate some partial wetting situations but suffers from two major drawbacks. First it leads to a singular force term in the right-hand-side of the momentum equation, prone to numerical instabilities. Second it requires a criterion on the film thickness to locate the contact line whose value can be difficult to calibrate and may be dependent on the configuration. To avoid the first drawback, an other approach has been developed in our group a few years ago. It consists of replacing the singular contact line force F cl δ cl n cl by a pressure term h∇Π cl that vanishes except in the vicinity of the contact line. The contact line pressure Π cl takes the following form [38,39,73]: with h * a cutoff length parameter used to locate the contact line. This expression is derived from a potential energy V cl detailed later in Section 3. This pressure term can be recasted in a conservative form in Eq.9b following the transformation [73]: We notice that this expression is very similar to Eq.14, especially with the presence of the spreading parameter in factor. The exponential term allows to locate the force at the interface. Indeed for h * going to 0, of the wall and 0 in the wetted area (h ≫ h * ). So formally, for h * small enough, we have: It shows that model Eq.17 can be seen as a regularization of model Eq.14. It is also strongly related to disjoining pressure models except that the potential energy expression V cl on which it relies is not derived from physical aspects related to molecular interactions. It is only asymptotically valid for small (h ≪ h * ) and large film thickness (h ≫ h * ). But again, like for the model proposed by Meredith et al., this model involves an adjustable parameter, h * , which is difficult to calibrate and may be sensitive to the mesh size and the numerical scheme [38].
To be completely exhaustive, a recent study, made by Härdi et al. [28], used the same approach as Meredith but in a Smoothed-Particle Hydrodynamic (SPH) framework. It seems to be a very promising mesh-free SPH method which presents some advantages. Firstly, the surface distribution of the lineic contact line force F cl is intrinsically linked to the smoothing processes of the SPH discretization. Another advantage is the facility to locate the force at the contact line. Indeed, it is not necessary to introduce a thickness criterion but simply locate the particles at the interface and apply the contact line force only to them. However, this last point highlights the main drawback of this method which is its computational cost. SPH discretization could be very expensive, especially to locate the particles that define the contact line, despite the optimized algorithm developed specifically for this operation [18,27].
Finally, it would be very interesting to develop a new model with contact line force without any thickness criterion and which could be used with more classical and versatile approaches such as the finite volume method. This explains our motivation to introduce a color function.

Presentation of a new macroscopic model with color function
To overcome the various limitations of the existing models we propose a new model based on a color function formulation. The idea is to use an interface indicator function illustrated in Figure 3 and defined as: This color function is propagated at the averaged velocity of the bulk flow with a simple transport equation: The main advantage of this function is that its gradient is zero everywhere except at the interface (see Figure   Figure 3. Illustration of the color function principle. 3) which allows to easily locate the contact line force. Following this consideration, we can replace the Dirac distribution in Eq.9b by: Moreover, if alpha is approximated by a smooth function (which is always the case at numerical level due to artificial diffusion), the transition from singular contact line force (lineic force) to smooth contact line force (surface force) is naturally realized with the gradient of the color function, following the same idea as Brackbill's work [5]. Finally, we propose the following new formulation for the contact line force: The difficulty is now to couple the transport equation Eq.20 of this color function with the System.9. In the following sections a methodology, based on the application of Hamilton's principle, is presented.

Hamilton's principle applied to macroscopic contact line force models
A key element of our macroscopic modeling is that the dewetting or partial wetting of a liquid film is simply the consequence of a more favorable energy state. The underlying assumption that we actually make is that it is possible to recover this state without solving the molecular scales but by applying only a macroscopic force. It is therefore interesting to use a method based on a variational principle involving the film energy to derive macroscopic models with contact line force. This is the purpose of this section where we present an approach based on Hamilton's principle to re-demonstrate the shallow-water system with contact line force introduced by Lallement et al. [38,73]. Before that we recall issues related to the Hamiltonian fluid dynamics.

Hamiltonian fluid dynamics: Hamilton's principle
Hamilton's principle is a fundamental principle of Hamiltonian mechanics published in 1834. Also known as the principle of least action, this principle states that the dynamics of any physical non-dissipative system can be determined by a single function called action [26,40]. It is therefore a very effective principle for the formulation of the equations that govern many conservative systems. This variational principle can take many forms depending on the physical domain studied. Originally this principle was based on a discrete reformulation of a physical problem using Lagrangian coordinates. It is this approach called particle-mechanics that was first used in fluid mechanics to re-demonstrate the equations of motion for a perfect fluid [20,53,66]. The strategy was to label each fluid particles, continuously distributed in space, with their initial positions. The position of each particle is known over time thanks to a function depending on the initial position which implies a complete mapping of the domain. This approach presents interesting aspects such as the fact that the conservation of mass is inherent to the labelling process of the particles [23,59]. However, as mentioned by Bretherton [6], this full Lagrangian version of Hamilton's principle for fluid mechanics has not really been used due to the fact that this formulation remains complex and is not necessarily the most appropriate for the description of a fluid. Herivel and Lin [29,41] have therefore proposed an Eulerian formulation of this variational principle adapted to fluid mechanics. The system is then no longer described by the position of each fluid particle but by the velocity at a point in the domain regardless of the particle present at that moment. However, this simpler approach requires some precautions. Indeed, the conservation of mass is no longer ensured and must therefore be constrained. Moreover, full Eulerian Hamilton's principle in its classic form does not allow the consideration of rotational flow. One of the great advances is due to Lin who proposed new constraints in the definition of the problem allowing to circumvent this restriction [41,42,50]. To summarize, the main difficulty of the application of Hamilton's principle lies in the formulation of the action associated to the physical system. In other words, to formulate a consistent expression of the energy.
The objective in this section is therefore to show that applying Hamilton's principle in an Eulerian representation allows to re-demonstrate the shallow-water equations with contact line force previously obtained by a simple analogy with other film models based on a disjoining pressure term by Lallement et al. (described in Section.2.3) [38,73]. This is also the method that we will apply later to derive the system involving the color function to model the capillary forces at the contact line.

Energy associated to a liquid film
As mentioned before, the first step in the application of Hamilton's principle is to give an expression of the energy associated to the liquid film to deduce its Lagrangian. To this end, we begin by considering an infinitesimal element of fluid of width dx, depth dy, thickness h(x) and average velocityū as shown in Figure 4. For shallow-water equations we start by defining the kinetic and potential gravity energy density per unit area, considering a constant density ρ for the clarity of the demonstration: In the rest of the paper the average velocity u will be simply noted u.
To avoid having an infinite energy in the case of a partially wetted infinite solid substrate, we prefer to apply a change of reference by setting the surface energy of the dry substrate to zero which implies to substract γ sg from the surface energy between the solid and the liquid. We thus obtain a slightly modified expression, compared to Eq.3, for the potential energy V c related to the capillary forces: Thus, we define the capillary surface energy density by: with V 0 c the 0 order term dominating at large scale, i.e. when the curvature of the film is negligible and V 1 c the first order term taking into account the curvature. This expression does not take into account the fact that the film does not cover the whole solid. In order to adress the case of a partially wetting film, Lallement et al. proposed to add to V c an additional term, noted V cl , which is activated only when the film thickness is less than a critical thickness noted h * . Its expression is the following: This potential energy is directly related to the contact line pressure (Eq.16) by: By construction, we see that when h and ∇h tend to 0 simultaneously, V c + V cl tend to 0 as well, which is the desired result in the absence of liquid on the wall.
Finally, the total surface energy associated to a liquid film is:

Application of Eulerian Hamilton's principle to redemonstrate an existing model with contact line force
We now apply the different steps of Hamilton's principle to re-demonstrate the momentum equation with contact line force of Lallement et al. [73]. Taking into account the energy expressions derived above, the Lagrangian associated to a partially wetting liquid film can be written: In addition, the film mass (or volume) conservation leads to the following equation : The constrained action of the film can be deduce as the spatial and temporal integration of this constrained Lagrangian between two fixed points in time and space (x i , t i ) and (x f , t f ): (31) where Ψ is the Lagrange multiplier associated to the constrains (Eq.30). After some development and integration by parts described in the appendix (Eq.111-115) and assuming that variations are null at the borders of the integration domain, the following expression for the action variation is obtained: To make the equations more compact and facilitate the calculations, it is convenient to introduce the following notations: and for a vector b: We apply the least action principle i.e. δJ f ilm = 0 and the fundamental lemma of calculus of variations to Eq.32 to obtain the following system of equations: Note: By developing Eq.35b we notice that the velocity field is necessarily irrotational (u = ∇Ψ). This nonphysical result is however perfectly consistent with the Eulerian formulation of the action that we have used (Eq.31). Indeed, as said above in Section 3.1, to avoid this limitation to irrotational velocity fields it would have been necessary to take into account the Lin constraints [41]. However, this does not prevent us from continuing our reasoning and to finally obtain the correct system of equations at the end as we will see.
Before continuing the calculation, using Eq.34a and Eq.34b, note that (see Appendix Eq.116 and Eq.117 for details): Then, we take the time derivative of Eq.35b: Using Eq.35a and Eq.35b: Moreover, note that: Thus, for the terms related to the kinetic energy in Eq.38: Finally, the momentum equation in the non-conservative form is given by: Or in detailed form: To switch to the conservative form we need to make some additional computations. First, note that V g and V cl depend only on the variable h, therefore ∇ (V g + V cl ) = ∂Vg ∂h + ∂V cl ∂h ∇h. Thus: Finally for the curvature term: A detailed development of the second right hand side term is given in the appendix (Eq.118) to obtain: This results in the following momentum equation written in conservative form: Or by replacing each term with its full expression: Finally, this is exactly the same system that was studied in [38,73]. This shows that this system derives from a principle of least action and that the existence of a conservative form of the momentum equation results from the invariance of the Lagrangian by translation in space.

Derivation of the new shallow-water model based on a color funtion
The objective of this section is to apply Hamilton's principle in Eulerian form, with a new expression of the film energy density using the color function α introduced in Section 2.4. As α may be interpreted as the film presence indicator function, the new total film energy density (per unit wall area) can be redefined by:

Application of Eulerian Hamilton's principle
Similarly we define the new Lagrangian as: In addition to the constraint of conservation of mass we add here a constraint on the dynamics of the color function in order to impose that α must satisfy the motion equation Eq.20, to express the fact that alpha is just transported by the mean film flow. For this we introduce a new Lagrange multiplier Φ related to this constraint. Finally, the action related to this new constrained Lagrangian is defined as.
After some development and integration by parts that we describe in the appendix (Eq.119-125) the following expression for the action variation is obtained: The least action principle i.e. δJ f ilm = 0 coupled to the fundamental lemma of calculus of variations lead to the following system of equations (using the notations Eq.33): Note: Contrary to the previous model, by developing Eq.52b the velocity field is no longer irrotational (u = ∇Ψ − ϕ∇α). The color function evolution is governed by Eq.52d which is exactly the same as a Lin's constraint [6]. Moreover the work of Seliger & Whitman [64] has shown that it is possible to use a variational principle with a single Lin constraint leading to a velocity field exactly of the above form which corresponds to Clebsch coordinate [10]. Thus, the Eulerian variational principle applicated to the color function based model is automatically generalizable to any kind of velocity field without the need for additional constraints.
Using equations Eq.36a and Eq.36b and taking the time derivative of Eq.52b we get: Using Eq.(52a), Eq.52b and Eq.52c, we find: Using the equation Eq.39, the terms related to the kinetic energy in Eq.54 are simplified: Finally the momentum equation in the non-conservative form is given by: Thus: To obtain the conservative form, we note first that V g depends only on the variable h, therefore ∇V g = ∂Vg ∂h ∇h. Thus: Then, if we consider constant surface energies, as in this study, the capillary potential energy V 0 c is constant, thus: Finally, for the curvature term: A detailed development of the second right hand term is given in the appendix (Eq.126) to obtain: Assembling all these terms results in the momentum equation in conservative form: Or by replacing each term with its full expression:

Source terms
As described in Section 2.3 the typical external forces that can be encountered are the wall shear stress τ w = µ ∂u ∂z (z=0) , the air shear stress τ a = µ g ∂ug ∂z (z=h) , the external pressure P g gradient and the tangential component of the gravity force g t if we consider the flow on an inclined plane. We can gather all these forces in a source term τ s added to the momentum equation such as in Eq.9b:

Final system
From Eq.50 and Eq.63, the final model reads as follows:

Model consistency
The aim of this section is to check that the new proposed model is physically consistent. First we will show that it is possible to recover the expected macroscopic behavior of static liquid films near the contact line (Young Dupré condition). Besides, we will prove that, in absence of external forces, the film energy is conserved locally, in the sense that the energy density per unit area satisfies a conservation equation (for smooth solutions).

Static equilibrium
If we consider a two-dimensional static film ( ∂ ∂t = 0, ∂ ∂y = 0 and u = 0 identically). The conservation of momentum equation (Eq.63) is then written : If we integrate this expression between two points x 0 and x 1 located on either side of the contact line which we suppose to be located at x = 0 as illustrated in Figure5. Knowing that by hypothesis h is identically zero for x > 0, Figure 5. Description of the geometry close to the contact line.
we obtain: Suppose that in the neighborhood of the contact line, h(x) tends to 0 and that ∂h ∂x (x) → p is finite when x → 0 − as in Figure 5. By making x 0 tend to 0 − , it gives : Assuming further that α(x 1 ) = 0 (dry area) and α(0 − ) = 1 (wet area), we finally obtain: Let θ be the angle between the contact line and the x axis. By definition of ∂h ∂x , it follows: The equation Eq.69 can therefore be rewritten in the following form: which is none other than the Young-Dupré relation (Eq.2). This proves the consistency of the model with the macroscopic condition of equilibrium of a static film.
If we are only interested in the macroscopic behaviour of a liquid film without necessarily representing very faithfully the shape of the film close to the contact line, the necessity of modeling the capillary forces related to the curvature is questionable. Especially since these terms can present difficulties during the numerical resolution [49]. This is why it is interesting to look at the behaviour of our model without the term related to the curvature of the interface, i.e. the term V 1 c in the expression of the capillary energy (see Eq.25). In this case, we can no longer assume that h is continuous (equal to 0) close to the contact line, otherwise according to Eq.67 it would follow: which is impossible except for S = 0, i.e. a static contact angle of 90 degrees. It is therefore necessary to suppose h(0 − ) > 0 to obtain: which gives, assuming α(0 − ) = 1 and α(x 1 ) = 0: This equation has a real solution only for S ≤ 0 (partially wetting film), which is the physically necessary condition for the existence of a static contact line. In fact this value of h(0 − ) is nothing else than the theoretical thickness of a flat partially wetting film at equilibrium far from the contact line (see Eq.4). Taking into account the curvature terms is therefore essential to find the correct shape of the film at the contact line and in particular the static contact angle. On the other hand, if we are only interested in very large scale phenomena in front of the film thickness, it may be sufficient to take into account only the term related to the contact force at the contact line. This simplifies the model and gets rid of the numerical difficulties associated with the curvature term.

Thermodynamic consistency
To show formally that our model implies the conservation of energy (for smooth solutions and without external forces), we first derive the equation satisfied by p = ∇h. For this, we rewrite the mass conservation equation by using the transport equation of the color function and takes the gradient of this new conservation equation: i.e α ∂h ∂t + α∇ · (hu) + h ∂α ∂t + u · ∇α = 0 i.e ∂h ∂t + ∇ · (hu) = 0 And thus: Then we introduce Legendre's like transformations of potential energies Eq.23 and Eq.25: Moreover, to facilitate the readability of the demonstration we take ρ = 1 without loss of generality. Thus we can rewrite our model under the following form: We introduce the thermodynamic variable φ conjugated to the variable U = (αh, αhu, α, p) T and rewrite the energy (Eq.48) of a liquid film: To deduce the energy equation we compute the scalar product between the System.78 and φ: We begin by developing each of these terms in the following subsections. Mass conservation term.
The details of the calculations are given in the Appendix Eq.127.
Momentum term.
Color function term.
Film thickness gradient term.
From equations Eq. 82, 83, 84 and 85, we now gather the corresponding physical terms together. Convective terms from Eq.82 and Eq.83.

− |u|
The details of the calculations are given in the Appendix Eq.128. Gravity terms from Eq.82, Eq.83 and Eq.84.
Second order internal capillary force terms from Eq.83 and Eq.85.
First order internal capillary force terms from Eq.83, Eq.84 and Eq.85.
To facilitate the calculation we will use Einstein summation convention and set φ p = αϕ p with ϕ p = for the development of the third term: We note the tensor related to the internal capillary forces (linked to the film curvature) is a symmetrical tensor. Its components read: Moreover, V 1 c depending only on the variable p, we write: It follows: By combining with Eq.89 it only remains for the curvature terms: Zero order internal capillary force terms from Eq.83 and Eq.84. Recall that L c,0 = −V 0 c . Then: Gathering the above intermediate results and using that: we get that: Finally we obtain the following conservative energy equation:

Comparison between the color function based model and the disjoining pressure like model
Let us now compare the new color function based model with the contact line force model developed by Lallement et al [38,73]. We will see that despite their great similarities in terms of physical content, the color function model brings great advantages, especially on the hyperbolic character of the backbone system without film curvature effects.

Location of the interface
If we go back to the contact line force model developed in [73] (see Eq.17) and make the parallel with the color function we may write: We note that α cl (h) = 0 for h ≪ h * and α cl = 1 for h ≫ h * . This means that α cl can be considered as a color function which depends on the film thickness h and parameterized by h * , h * being the thickness threshold that distinguishes the dry from the wet area of the wall. In the model based on the color function, it is not necessary to define such a threshold because we use the additional variable α that we transport to locate the position of the contact line at each instant. We will see that the fact of being able to locate the interface without an adjustable non-physical parameter is not only an advantage for numerical purposes but also from the theoretical point of view.

Hyperbolicity analysis
For numerical purposes it is interesting to compare the hyperbolicity of the two models in the limit of negligible curvature effects as in this limit both models can be written as first order systems of conservation laws. In this case, from Section 3, 1D system from Lallement et al. [73] can be written under the following quasi-linear form.
with U = (h, ρhu) T the vector of conservative variables (see Section5.2) and J the Jacobian matrix of the flux φ cl : defined by: The characteristic polynomial of J is given by: We find the eigenvalues of the system by examining its discriminant which can be written: This shows that the hyberbolicity property is directly related to the convexity of the total energy of the system which is a well known result from the Godunov-Mock theorem [25]. The gravity energy is strictly convex, so it makes the system hyperbolic and therefore contributes to the system stability. On the other hand, for the contact line surface energy density we have: Thus, the contact line energy is strictly concave. Due to this property, for small thicknesses such as (h ≈ h * ) the discriminant becomes negative which leads to imaginary eigenvalues and thus to a local loss of hyperbolicity close to the contact line. For thicknesses such as h ≫ h * the contact line energy term is negligible compared to gravity and the discriminant is positive, which leads to a system with two distinct real eigenvalues characteristic of a hyperbolic system: However, it does not imply that the full system with first order capillary terms is ill-defined. Indeed, the curvature third order term (not taken into account in this study) stabilizes the system thanks to the damping of high frequency instabilities by the capillary forces [38]. Looking at the discriminant Eq.104, to have a completely hyperbolic system, whatever the thickness of the film, it would have been necessary to have a convex energy V cl . However, V cl is a function that must tend to S in the limit of zero thickness and to zero for very large thickness (see Section 3.2). Thus, for partial wetting configuration (S < 0) this function has no choice but to be non-convex, contrary to the case of total wetting (S > 0), as illustrated in Figure 6. This is reminiscent of the concave form of the disjoining energy of the microscopic films discussed in introduction (cf Figure 2) which leads to spinodal dewetting. However, the films we aim to study are sufficiently thick to neglect the long ranges interactions, except very close to the contact line, and therefore spinodal dewetting is very far from our scale study. But, it is interesting to understand we show here that the strategy used in the construction of this model to capture the physics close to the contact line is closely related to the physics of spinodal dewetting and, thus, that the macroscopic construction of this contact line force model is based on microscopic physical foundations, to capture macroscopic dewetting and partial wetting phenomena.
If we now consider the model with color function, the problem is a little different. Indeed as seen previously, the addition of the variable α allows to get rid of an additional contact line energy (like V cl depending on the thickness and necessarily concave) to model partial wetting. This completely changes the study of the system hyperbolicity. We write the corresponding 1D system under the following quasi-linear form: with U = (αh, αρhu, α) T : The characteristic polynomial of A is given by: For h > 0, there are 3 real distinct eigenvalues whose expression is: Thus, contrary to the model with contact line force, the one with color function is hyperbolic whatever the thickness of the film which is a great advantage for the development of a robust and accurate numerical solver suitable for this model.

Summary
To summarize, we gather all that has been demonstrated in this section in Table 1. This highlights the main theoretical improvements brought by this new formulation with color function in comparison with contact line force model.
Dewetting energy term  Table 1. Summary of the main differences between the two shallow-water type models for the simulation of partial wetting phenomena.

Conclusion
In this study we have highlighted theoretical difficulties related to the study of partial wetting phenomena which partly explain why it is still an active subject of research. However, their understanding and modeling are essential because situations where a liquid film wets a solid surface are very often encountered in many industrial applications. The physics of wetting having microscopic origins the majority of existing models are based on the fine resolution of what happens close to the contact line. These models are of course not suitable for typical scales encountered in most industrial applications. This is why our interest here was to derive a new macroscopic model valid in the lubrication regime, that allows to simulate the global behaviour of a partially wetting liquid film without the need for very fine scales of resolution. This kind of model already exists and can be placed in the category of contact line force models. Their ability to simulate wetting and dewetting phenomena have been proven on some test cases. However, these models suffer from the need of a thickness criterion to concentrate the effect of the contact line force at the interface. The major advantage of the model we proposed in this paper is that it does not depend on any adjustable parameter to localize the interface. It is based on the use of a color function, transported at the average speed of the film and whose gradient is non-zero only at the contact line. Moreover, we have shown that in the absence of higher order terms related to the curvature of the air-film interface, the new model is hyperbolic.
Beyond this new approach for the modelling of partial wetting phenomena; we have also presented a methodology to formulate a thermodynamic consistent model. This methodology is based on a Hamiltonian approach of fluid mechanics. Hamiltonian mechanics was initially based on a discrete Lagrangian approach of the physical systems but has since been adapted to be used with Eulerian formalism. Thus, in this study we have used Hamilton's principle in a full Eulerian form to recover the disjoining pressure like model already introduced by Lallement et al. [73] and we then used this new methodology to derive our new color function based model. The interest of this approach, apart from its elegance, lies in the fact that it allows to automatically derive a backbone model (i.e. without dissipative and external force terms) which is consistent with the principle of conservation of energy, the only input being a "reasonable" expression for the free energy of the film as a function of the chosen thermodynamic variables.
In a nutshell, the main contribution of this work is to provide a new macroscopic shallow-water type model for partial wetting applications particularly suitable for numerical simulations. The next step of this work is to develop a finite volume numerical scheme adapted to this model and to perform simulations of partial wetting phenomena. This work is currently in progress and will be presented in a future paper.

Acknowledgements
The financial support from the ONERA french aerospace lab and the doctoral MEGEP school is gratefully acknowledged.