Christine Bernardi
email: bernardi@ann.jussieu.fr

Sarra Maarouf
email: maarouf@ljll.math.upmc.fr

Driss Yakoubi
email: yakoubi@giref.ulaval.ca

Finite element discretization of two immiscible fluids with surface tension

We consider a model for the flow of two immiscible fluids where the surface tension between both of them is taken into account. We first write the variational formulation of the problem and investigate its well-posedness. Next, we consider a finite element discretization of it and prove optimal a priori error estimates. Numerical experiments confirm its good properties.

Résumé. Nous considérons un modèle pour l'écoulement de deux fluides immiscibles où la tension de surface entre les deux fluides est prise en compte. Nous écrivons tout d'abord la formulation variationnelle de ce problème et vérifions qu'il est bien posé. Puis nous en proposons une discrétisation par éléments finis et prouvons des estimations d'erreur a priori qui s'avèrent optimales. Quelques expériences numériques confirment ses bonnes propriétés.

Introduction.

Numerical methods for solving free surface problems are of great importance in many engineering applications and still a challenging field. Problems with free surfaces appear in immiscible multi-fluids problems, motion of glaciers, fluid structure interactions, blood in moving arteries and many other domains.

For two immiscible Newtonian fluids, the time-dependent Navier-Stokes equations govern the fluid motion in both fluids. Each fluid has a constant density and viscosity. As we are dealing with flow with more than one fluid, the surface tension at the fluid interface has to be accounted for. The surface tension is modeled as a body force concentrated at the interface by employing the Continuum Surface Force (CSF) model of Brackbill et al. [START_REF] Braekbill | Zemach -A continuum method for modeling surface tension[END_REF]. The CSF model allows us to treat the dynamic boundary condition at the interface implicitly. So, the main difficulty in this problem is that, of course, each fluid flows through a time-dependent domain and that the interface is involved in the Navier-Stokes equations via the surface tension of the fluids. Moving interfaces can be handled with the Level Set method introduced S. Osher and al. [START_REF] Osher | Sethian -Fronts propagating with curvature dependent speed: algorithms based on Hamilton-Jacobi formulations[END_REF], [START_REF] Sussman | Osher -A level set approach for computing solutions to incompressible two-phase flow[END_REF], see also [START_REF] Osher | Fedkiw -Level Set Methods and Dynamic Implicit Surfaces[END_REF] for a more general review of all these results.

It should be noted that the level set method is one of several interface-tracking techniques used routinely for two immiscible incompressible flow simulations. Others are the front-tracking method, see for instance [START_REF] Popinet | Zaleski -A front-tracking algorithm for accurate representation of surface tension[END_REF], [START_REF] Tryggvason | Jan -A front-tracking method for the computations of multiphase flow[END_REF], the boundary integral method [START_REF] Hou | Boundary integral methods for multicomponent fluids and multiphase materials[END_REF], and the volume-of-fluid (VOF) method [START_REF] Braekbill | Zemach -A continuum method for modeling surface tension[END_REF], [START_REF] Lafaurie | Zanetti -Modelling merging and fragmentation in multiphase flows with SURFER[END_REF], [START_REF] Meier | Smith -A novel technique for including surface tension in PLIC-VOF methods[END_REF]. Unfortunately, the main advantages of the level set method compared with other techniques are • the attractive simplicity of its mathematical formulation and computing, and • the ability of the method to simulate complex interfacial flows with strong surface tension effects.

In the level set method, the surface tension force, which is essentially a surface force, is traditionally modeled as a distributed body force, though concentrated in a band around the interface and arranged in such a way that the force has a maximum on the interface and decays rapidly with distance from it, see [START_REF] Sussman | Osher -A level set approach for computing solutions to incompressible two-phase flow[END_REF], [START_REF]Sethian -Level Set Methods and Fast Marching Methods: Evolving Interfaces[END_REF]. Furthermore, in [START_REF] Doyeux | Ismail -Simulation of two-fluid flows using a finite element/level set method. Application to bubbles and vesicle dynamics[END_REF], the size of this band ε can be taken proportional to the mesh size h as ε ≈ 1.5 h. Thus, the variation of the surface tension across the interface may be described in terms of a regularized (smoothed) delta function with compact support. This approach removes the interface singularity from the standard continuum fluid flow equations, and ultimately allows the surface tension to be modeled using standard numerical techniques on Eulerian grids. The jump in phase properties across the interface, such as density and viscosity, is correspondingly modeled using a regularized Heaviside function.

So the full model consists in the time-dependent Navier-Stokes equations and a transport equation for the interface. A detailed description of a similar model can be found in [START_REF] Frey | Ta -A numerical scheme based on level set method and anisotropic mesh adaptation for solving two-fluid flows[END_REF], see also the references therein.

To our knowledge, the literature that studied the subject has been focussed on numerical validation and computational topics, while numerical analysis, for instance using fnite element method does not appear,

In this work, we investigate the main properties of the model and in particular the existence of a solution in the non-realistic case where the velocity satisfies homogeneous boundary conditions. Next, we propose a discretization of it by the characteristics method in time (indeed each equation contains a convection term) and standard conforming finite elements in space. We perform the a priori analysis of the discrete problem and prove nearly optimal estimates of the error. We finally present some numerical experiments that confirm the interest of both the model and its discretization.

The outline of the paper is as follows:

• In Section 2, we explain and write the model.

• In Section 3, we write its variational formulation and investigate its well-posedness.

• Section 4 is devoted to the description of the discretization and the well-posedness of the discrete problem.

• Its a priori analysis is performed in Section 5.

• In Section 6, we present some numerical experiments.

Mathematical description.

Let Ω be a bounded connected open set in R d , d = 2 or 3, with a Lipschitz-continuous boundary ∂Ω. Let also T f be a positive real number. If we assume a domain Ω with two immiscible fluids F i , i = 1, 2 , then the time dependent subdomains Ω i (t), t ∈ (0, T f ) are bounded by an external boundary ∂Ω and by a dynamic interface Γ(t), see Figure 1. We assume that both fluids are homogeneous and therefore the physical properties are constant in each Ω i which are a bounded connected domain with a Lipschitz-continuous and connected boundary ∂Ω i (t):

Ω = Ω 1 (t) ∪ Ω 2 (t), Ω 1 (t) ∩ Ω 2 (t) = ∅. (2.1)
We also define the interface between the two fluids Then, the standard model for isothermal two immiscible flows can be described by the incompressible Navier-Stokes system

Γ(t) = ∂Ω 1 (t) ∩ ∂Ω 2 (t). ( 2
   ρ ∂u ∂t + u • ∇u + ∇p -∇ • (ν∇u) = f + f σ in Ω × (0, T f ), ∇ • (ρu) = 0 in Ω × (0, T f ), (2.3)
which contain an additional force term f σ due to the surface tension σ at the free interface Γ(t). The unknowns are the velocity u (more precisely, each u i = u| Ω i is the velocity of the fluid F i , i = 1, 2) and the pressure p, f denotes any body force such as gravitational acceleration. Note that the density ρ as well as the viscosity ν are variable and discontinuous, that is

ρ =    ρ 1 in Ω 1 × (0, T f ), ρ 2 in Ω 2 × (0, T f ), and ν =    ν 1 in Ω 1 × (0, T f ), ν 2 in Ω 2 × (0, T f ),
where each ρ i and ν i are nonnegative constants. Thus, the second equation in (2.3) can equivalently be written

∇ • u = 0 in Ω × (0, T f ). (2.4)
Surface tension effects are taken into account through the following force balance at the interface Γ(t)

u 1 = u 2 on Γ(t). and (ν 1 ∇u 1 -p 1 I) • n + (ν 2 ∇u 2 -p 2 I) • n = σ κ n on Γ(t).
Here n is the unit normal at the interface pointing for instance into Ω 2 , σ ≥ 0 is the surface tension coefficient and κ = ∇ • n is the curvature of the interface Γ(t) . The first condition implies continuity of the velocity across the interface, whereas the second describes the force balance on Γ(t). To handle the curvature term, it is often to rewrite it as a volume force, that means

f σ = σ κ n δ σ
where δ σ is the Dirac delta function localizing the surface tension forces to the interface Γ(t). In Lafaurie et al. [START_REF] Lafaurie | Zanetti -Modelling merging and fragmentation in multiphase flows with SURFER[END_REF], the authors propose an other shape but still equivalent to this force by introducing the projection operator in the tangent plan: In ⊗ n and then writing

f σ = -∇ • T, such that T = σ(I -n ⊗ n) δ σ .
According to the applied Continuous Surface Force (CSF) approach, see [START_REF] Braekbill | Zemach -A continuum method for modeling surface tension[END_REF], we introduce the following smooth regularization of the Heaviside function for a small ε > 0,

H ε (ψ) =      0 if ψ < -ε, 1 2 1 + ψ ε + 1 π sin( πψ ε ) if -ε ≤ ψ ≤ ε, 1 if ψ > ε.
(2.5)

On the other hand, following [START_REF] Sussman | Osher -A level set approach for computing solutions to incompressible two-phase flow[END_REF], a level set approach is employed to capture the interface Γ(t). So, we consider the oriented distance function ϕ to Γ(t) (also called level set function). The interface between the two fluids is then identified as the zero level of ϕ, more precisely

Γ(t) = x ∈ Ω; ϕ(x, t) = 0 .
Then, the smoothed density ρ and the viscosity ν are now given by

ρ(ϕ) = ρ 1 + (ρ 2 -ρ 1 )H ε (ϕ), (2.6)
ν(ϕ) = ν 1 + (ν 2 -ν 1 )H ε (ϕ). (2.7)
Since ϕ is a distance function (so that ∇ϕ is a normal vector to the interface Γ(t)) and using the well-known expressions, this tensor T can also be written as a function of ϕ, as follows:

T (ϕ) = σ dH ε (ϕ) dϕ (I -∇ϕ ⊗ ∇ϕ), (2.8)
where ∇ϕ ⊗ ∇ϕ stands for the tensor with coefficients

∂ x i ϕ ∂ x j ϕ, 1 ≤ i, j ≤ d.
We recall that the unit normal n to the interface is classically obtained via ϕ : On the curve Γ(t) with equation ϕ = 0, n = ∇ϕ |∇ϕ| .

Finally, the governing equations which describe two immiscible fluids with surface tension can be writing as a coupled time-dependent Navier-Stokes equations with transport equation:

           ρ(ϕ) ∂u ∂t + u • ∇u + ∇p -∇ • (ν(ϕ)∇u)) + ∇ • T (ϕ) = f in Ω × (0, T f ), ∇ • u = 0 in Ω × (0, T f ), ∂ϕ ∂t + u • ∇ϕ = 0 in Ω × (0, T f ).
(2.9)

The last equation means that the interface Γ(t) is convected by the fluid, see [12, eq. ( 7)].

A very similar model is considered in [9, eq. (2.9)], see also [START_REF]Milcent -Une approche eulérienne du couplage fluide-structure, analyse mathématique et applications en biomécanique[END_REF].

To make the problem complete, we enforce suitable boundary and initial conditions. The Dirichlet condition is imposed to the velocity and the level set function

u = u D on ∂Ω and ϕ = ϕ D on Γ u , (2.10)
where Γ u denotes the part of boundary where the fluid goes in

Γ u = {x ∈ ∂Ω; u(x, t) • n(x) < 0},
and the initial conditions obviously read u(x, 0) = u 0 (x) and ϕ(x, 0) = ϕ 0 (x) for a.e. x ∈ Ω.

(2.11)

The existence of Γ u and the boundary condition on ϕ come from the nature of the third line of (2.9) which is a hyperbolic equation. However, in the case of homogeneous Dirichlet condition on ∂Ω for u, no essential boundary conditions are needed on the level set function.

Certain compatibility conditions must be satisfied by the given data, mainly

∇ • u 0 = 0 in Ω and u D (•, 0) • n = u 0 • n on ∂Ω, (2.12)
in order to ensure the smoothness of the solution. From now on, we take ϕ negative on Ω 1 (t) and positive on Ω 2 (t), whence some compatibility conditions on ϕ 0 and ϕ D .

Variational formulation and well-posedness.

Of course, our dream would be that problem (2.9) provided with the initial conditions (2.10) and the boundary conditions (2.11) is well-posed. However this seems impossible, for the two following reasons: 1) very few existence results are known for the transport equation (third line in (2.9)) when the velocity u does not vanish on the boundary and they require too much regularity of the domain; 2) the uniqueness of the solution of the time-dependent Navier-Stokes equation (first and second lines in (2.9)) is actually unknown in dimension d = 3. So, we restrict ourselves to a simpler case.

In what follows, we use the whole scale of Sobolev spaces W m,p (Ω), with m ≥ 0 and 1 ≤ p ≤ +∞, equipped with the norm • W m,p (Ω) and seminorm | • | W m,p (Ω) , with the usual notation H m (Ω) when p = 2. We also need the space H 1 0 (Ω) of functions in H 1 (Ω) vanishing on ∂Ω. For any separable Banach space E equipped with the norm • E , we denote by C 0 (0, T ; E) the space of continuous functions from [0, T ] with values in E. For each integer m ≥ 0, we also introduce the space H m (0, T ; E) as the space of measurable functions on ]0, T [ with values in E such that the mappings: v → ∂ t v E , 0 ≤ ≤ m, are square-integrable on ]0, T [. We refer to [START_REF] Adams | Fournier -Sobolev Spaces[END_REF] for the main properties of all these spaces.

Let us first consider the transport equation. In the case where u D is not zero, the only existence result can be found in [6, Thm VI.1.6] but it requires too much regularity of the domain Ω and of the velocity u. So, even if the case u D • n = 0 was treated in [START_REF] Diperna | Lions -Ordinary differential equations, transport theory and Sobolev spaces[END_REF]Section IV.4], from now on we assume that

u D = 0. (3.1)
With this condition, it is rather easy to prove the existence of a solution via the characteristics method, see [2, Thm II.7.6]. We are thus in a position to state the result concerning this equation.

Proposition 3.1. Assume that the function u is divergence-free and satisfies

u ∈ C 0 (0, T f ; H 1 0 (Ω) d ) ∩ C 0 (0, T f ; W 1,∞ (Ω) d ). (3.2)
Then, for any datum ϕ 0 in L 2 (Ω), the problem

   ∂ϕ ∂t + u • ∇ϕ = 0 in Ω × (0, T f ), ϕ(•, 0) = ϕ 0 in Ω, (3.3
)

admits a unique solution ϕ in C 0 (0, T f ; L 2 (Ω)). Moreover, this solution satisfies sup 0≤t≤T f ϕ(•, t) L 2 (Ω) ≤ ϕ 0 L 2 (Ω) . (3.4)
Proof. We establish successively the existence, the uniqueness and the stability property (3.4) of the solution.

1) Existence: Owing to the property (3.2) of u, applying the Cauchy-Lipschitz theorem [START_REF]Schatzman -Analyse numérique[END_REF]Th. 21.1] yields that there exists for every pair (

x, t) in Ω × [0, T f ], a characteristic function X = X(x, t; •) in C 0 (0, T ; R d ) solution of the ordinary differential equation dX dt = u • X in Ω × (0, T f ), X(x, t; t) = x in Ω.
Thus, it is readily checked that the function ϕ defined by

ϕ(x, t) = ϕ 0 (X(x, 0; t)), (3.5)
is a solution of problem (3.3).

2) Uniqueness: Let ϕ 1 and ϕ 2 be two solutions of problem (3.3). Setting ϕ = ϕ 1ϕ 2 and multiplying the difference of the two equations by ϕ gives 1 2

∂ϕ 2 ∂t + 1 2 u • ∇ϕ 2 = 0.
Since u is divergence-free and vanishes on the boundary, integrating on Ω yields 1 2

d ϕ 2 L 2 (Ω) dt = 0, whence, since ϕ(•, 0) is zero, ϕ(•, t) is zero. Thus, ϕ 1 and ϕ 2 coincide.
3) Stability: Multiplying as previously problem (3.3) by ϕ and integrating on Ω leads to

d ϕ 2 L 2 (Ω) dt = 0,
whence property (3.4).

It follows from formula (3.5) that the regularity of the solution ϕ only depends on that of u and ϕ 0 . So we state the regularity properties that we need to study the full problem.

Corollary 3.2. If the assumptions of Proposition 3.1 hold, for any real number p, 1 ≤ p < ∞, and any datum ϕ 0 in W 1,p (Ω), the solution ϕ of problem (3.3) belongs to C 0 (0, T f ; W 1,p (Ω)).

We now consider the Navier-Stokes equations. We first write their variational formulation:

Find u in L 2 (0, T ; H 1 0 (Ω) d ) ∩ C 0 (0, T ; L 2 (Ω) d ) and p in L 2 (0, T ; L 2 (Ω)) such that u(•, 0) = u 0 , (3.6)
and, for a.e. t in [0, T ],

∀v ∈ H 1 0 (Ω) d , Ω ρ(ϕ)(x) ∂u ∂t + u • ∇u (x)v(x) dx - Ω (∇ • v)(x)p(x) dx + Ω ν(ϕ)(x)∇u(x) : ∇v(x) dx - Ω T (ϕ)(x) : ∇v(x) dx = f , v , ∀q ∈ L 2 (Ω), Ω (∇ • u)(x)q(x) dx = 0, (3.7)
where •, • denotes the duality pairing between H 1 0 (Ω) d and its dual space H -1 (Ω) d . The arguments for proving its equivalence with the first two lines of (2.9) are fully standard, so we skip them.

First, we note from the definitions (2.6) and (2.7) that, for any ϕ,

min{ρ 1 , ρ 2 } ≤ ρ(ϕ) ≤ max{ρ 1 , ρ 2 }, min{ν 1 , ν 2 } ≤ ν(ϕ) ≤ max{ν 1 , ν 2 }. (3.8)
Next, since dH ε dϕ is bounded independently of ϕ, for any ϕ in L 4 (0, T f ; W 1,4 (Ω)), the quantity T (ϕ) belongs to L 2 (0, T f ; L 2 (Ω) d×d ). So all the terms in problem (3.7) are well-defined. Moreover, thanks to all these properties, the existence of a solution for problem (3. Let Ω be a bounded smooth domain in R 3 , and let p be a real number > 3. We consider data (i) ϕ 0 in W 2,p (Ω) such that |∇ϕ 0 | is larger than a positive constant in a neighbourhood of {x ∈ Ω; ϕ 0 (x) = 0};

f ; W 1,4 (Ω)). Then, for any datum f in L 2 (0, T f ; H -1 (Ω) d ) and u 0 in L 2 (Ω) d satisfying ∇ • u 0 = 0 in Ω and u 0 • n = 0 on ∂Ω, (3.9
(ii) f equal to σd dH ε (ϕ) dϕ , (iii) u 0 in W 2,p (Ω) d ∩ W 1,p 0 (Ω) d which is divergence-free in Ω.
Then, there exists a positive number T * only depending on the initial data such that problem

(3.3) -(3.6) -(3.7) has a solution in Ω × [0, T * ]. Moreover, this solution satisfies ϕ ∈ L ∞ (0, T * ; W 2,p (Ω)), u ∈ L ∞ (0, T * ; W 1,p 0 (Ω) d ) ∩ L p (0, T * ; W 2,p (Ω) d ), ∇p ∈ L p (0, T * ; L p (Ω) d ).
Even if neither the assumptions of this proposition nor the results are very realistic, it leaves the hope that the problem we work with is not fully incoherent.

The discrete problem and its well-posedness.

To describe the discrete problem, we first need some notation. Next, we write it and prove its well-posedness.

Some notation.

Since we intend to work with non-uniform time steps, we introduce a partition of the interval [0,

T f ] into subintervals [t n-1 , t n ], 1 ≤ n ≤ N , with 0 = t 0 < t 1 < • • • < t N = T f .
We denote by τ n the time step t nt n-1 , by τ the N -tuple (τ 1 , . . . , τ N ) and by |τ | the maximum of the τ n , 1 ≤ n ≤ N .

We assume that Ω is a polygon (d = 2) or a polyhedron (d = 3). For 0 ≤ n ≤ N , let (T n h ) h be a regular family of triangulations of Ω (by triangles or tetrahedra), in the sense that, for each h:

• Ω is the union of all elements of T n h ; • The intersection of two different elements of T n h , if not empty, is a vertex or a whole edge or a whole face of both of them;

• The ratio of the diameter h K of any element K of T n h to the diameter of its inscribed circle or sphere is smaller than a constant independent of h and n. As usual, h n stands for the maximum of the diameters h K , K ∈ T n h . We also denote by V n h the set of all vertices of elements of T n h . In what follows, c, c , . . . are generic constants which may vary from line to line but are always independent of the τ n and h n .

We have decided to work with the Taylor-Hood finite element introduced in [START_REF] Hood | Taylor -A numerical solution of the Navier-Stokes equations using the finite element technique[END_REF]]. So, denoting by P k (K) the space of restrictions to K of polynomials with d variables and total degree ≤ k, we consider the finite element spaces

Y n h = v h ∈ H 1 (Ω); ∀K ∈ T n h , v h | K ∈ P 2 (K) , X n h = (Y n h ) d ∩ H 1 0 (Ω) d , (4.1)
and

M n h = q h ∈ H 1 (Ω); ∀K ∈ T n h , q h | K ∈ P 1 (K) . (4.2)
We denote by I n h the Lagrange interpolation operator at the nodes of V n h with values in M n h .

Discretization of the transport equation.

For reasons which appear later on, we still work with homogeneous boundary condition on the velocity, i.e. we assume that (3.1) holds. We have decided to use the characteristics method to discretize the transport equation, as introduced and firstly analyzed in [START_REF] Pironneau | On the transport-diffusion algorithm and its applications to the Navier-Stokes equations[END_REF].

Assuming that the datum ϕ 0 is continuous on Ω, we simply define

ϕ 0 h = I 0 h ϕ 0 . (4.3)
Next, at time t n , assuming that u n-1 h and ϕ n-1 h are known, we define for all x in V n h

x n h = x -τ n u n-1 h (x), (4.4)
thus, we define φn h by φn

h (x) = ϕ n-1 h (x n h ), (4.5)
and finally

ϕ n h = I n h φn h . (4.6)
The simplicity of this algorithm is obvious. Note simply that: (i) If x belongs to ∂Ω, the fact that u n-1 h belongs to H 1 0 (Ω) d implies that x n h is equal to x, hence belongs to Ω; (ii) otherwise, τ n must be chosen small enough for x n h to be in Ω. This condition is made precise later on.

Remark 4.1. Since the function ϕ n h is continuous on Ω, we can define the curve

Γ n h = x ∈ Ω; ϕ n h (x) = 0 . (4.7)
We hope that this curve is a good approximation of the interface Γ(t) at time t = t n . This will be investigated later on.

Discretization of the Navier-Stokes equations.

This discretization relies on two arguments: using the characteristics method for handling the convection term and a Galerkin method for the rest of the equation. As previously, assuming that the datum u 0 is continuous on Ω, we simply define

u 0 h = I 0 h u 0 . (4.8)
Next, at time t n , 1 ≤ n ≤ N , assuming that u n-1 h and ϕ n h are known, still relying on equation (4.4), we define ũn-1 h as the interpolate in (M n h ) d of the values u n-1 h (x n h ) at all nodes of V n h which belong to Ω. Then, the discrete problem reads:

Find u n h in X n h and p n h in M n h such that ∀v ∈ X n h , Ω ρ(ϕ n h )(x) u n h -ũn-1 h τ n (x)v(x) dx - Ω (∇ • v)(x)p n h (x) dx + Ω ν(ϕ n h )(x)∇u n h (x) : ∇v(x) dx - Ω T (ϕ n h )(x) : ∇v(x) dx = f (t n ), v , ∀q ∈ M n h , Ω (∇ • u n h )(x)q(x) dx = 0. (4.9)
4.4. Well-posedness of the discrete problem.

To prove the well-posedness of problem (4.4) -(4.5) -(4.6), we must check that the x n h defined in (4.4) belongs to Ω. This is proved in [27, Prop. 1], see also [START_REF] Pironneau | Tabata -Stability and convergence of a Galerkin-characteristics finite element scheme of lumped mass type[END_REF]Rem. 3], when the following condition holds

τ n < 1 u n-1 h W 1,∞ (Ω) d . (4.10)
However, we have no a priori control on the quantity u n-1 h W 1,∞ (Ω) d . So we prefer to use the following modified algorithm: 1) if the quantity defined in (4.4) belongs to Ω, take x n h equal to it; 2) otherwise, take x n h equal to

x n h = x -t u n-1 h (x), (4.11)
for the smallest t in [0, τ n ] such that xt u n-1 h (x) belongs to ∂Ω. Note that the modification introduced above is only used for a small number of nodes x which are very near to ∂Ω. We are thus in a position to state the next proposition. . We also observe that the quantity T (ϕ n h ) belongs to L 2 (Ω) d×d . So the following result is readily checked. Proposition 4.3. For any datum u 0 continuous on Ω, problem (4.9) has a unique solution (u n h , p n h ) in X n h × M n h , up to an additive constant on the pressure p n h . Proof. When ϕ n h (and T (ϕ n h )), ũn-1 h and f are known, this problem results in a square linear system. So it suffices to prove that the only solution of the problem

∀v ∈ X n h , Ω ρ(ϕ n h )(x) u n h (x) τ n v(x) dx - Ω (∇ • v)(x)p n h (x) dx + Ω ν(ϕ n h )(x)∇u n h (x) : ∇v(x) dx = 0, ∀q ∈ M n h , Ω (∇ • u n h )(x)q(x) dx = 0,
is zero, up to an additive constant on the pressure p n h . By taking v equal to u n h , combining the two previous equations and using (3.8), we easily derive that u n h is zero. Thus, the fact that p n h is a constant follows from the inf-sup condition, see for instance [14, Chap. II, Thm 4.2],

∀q ∈ M n h , Ω q(x) dx = 0, sup v∈X n h Ω (∇ • v)(x)q(x) dx v H 1 (Ω) d ≥ β q L 2 (Ω) , (4.12)
where the constant β is positive.

Remark 4.4. The main advantage of the discretization that we propose is that it uncouples the computation of the ϕ n h and of the u n h . However, this can lead to a lack of convergence for too large time steps or in the case of time oscillations of the solution. To remedy this, we can use an iterative algorithm as follows: at each time t n , 1) define ϕ n,0 h as equal to ϕ n-1 h and u n,0 h as equal to

u n-1 h ; 2) setting x n,k h = x -τ n u n,k-1 h (x),
(or the modified version (4.11)), we define ϕ n,k h as the interpolate in M n h of the values ϕ n,k-1 (x n,k h ) at all nodes of V n h and similarly ũn,k-1 h as the interpolate in M n h of the values u n,k-1 (x n,k h ) at all nodes of V n h which belong to Ω, next we solve

Find u n,k h in X n h and p n,k h in M n h such that ∀v ∈ X n h ; Ω ρ(ϕ n,k h )(x) u n,k h -ũn,k-1 h τ n (x)v(x) dx -λ k Ω (∇ • v)(x)p n,k h (x) dx + λ k Ω ν(ϕ n,k h )(x)∇u n,k h (x) : ∇v(x) dx -λ k Ω T (ϕ n,k h )(x) : ∇v(x) dx = λ k f (t n ), v , ∀q ∈ M n h , Ω (∇ • u n,k h )(x)q(x) dx = 0.
The λ k are positive parameters, chosen in order to ensure the consistency of the algorithm.

3) For a fixed integer K, we set: ϕ n = ϕ n,K and u n = u n,K , and go back to step 1.

Of course, the parameter K must be chosen small enough, namely equal to 2 or 3, in order that the cost of this algorithm remains reasonable. Its analysis is nearly the same as for the previous problem (4.4) -(4.5) -(4.6) -(4.9), so that we skip it.

A priori error analysis.

In this section, we still work with u D = 0 and we admit that, for all x in V n h , the x n h defined in (4.4) belong to Ω (indeed, it is not so difficult in practice to satisfy condition (4.10): when u n-1 h W 1,∞ (Ω) d is larger, it is possible to work with a smaller τ n in order to enforce this condition).

We first state a stability property for the discrete transport equation. Lemma 5.1. For any datum ϕ 0 continuous on Ω, the next property holds for

1 ≤ n ≤ N ϕ n h L ∞ (Ω) ≤ ϕ 0 L ∞ (Ω) . (5.1)
Proof. Since each ϕ n h belongs to M n h , hence is piecewise affine, we have

ϕ n h L ∞ (Ω) = max x∈V n h |ϕ n h (x)|. (5.2)
Thus it follows from (4.5) and (4.6) that

ϕ n h L ∞ (Ω) ≤ ϕ n-1 h L ∞ (Ω)
.

By noting that

I 0 h ϕ 0 L ∞ (Ω) ≤ ϕ 0 L ∞ (Ω)
, the desired result is derived by induction on n.

The error estimate is a little more tricky, even in the L ∞ -norm. We refer to [25, Prop. 1] for a similar result in a different case. Proposition 5.2. Assume that the solution ϕ of problem (3.3)

satisfies ϕ ∈ C 0 (0, T f ; W 2,∞ (Ω)) ∩ W 1,∞ (0, T f ; L ∞ (Ω)).
(5.3)

Then the following estimate holds for

1 ≤ n ≤ N ϕ(•, t n ) -ϕ n h L ∞ (Ω) ≤ C h 2 n + τ n (1 + u n-1 h L ∞ (Ω) d ) + ϕ(•, t n-1 ) -ϕ n-1 h L ∞ (Ω) , (5.4)
where the constant C only depends on the regularity of ϕ.

Proof. We start with the triangle inequality

ϕ(•, t n ) -ϕ n h L ∞ (Ω) ≤ ϕ(•, t n ) -I n h ϕ(•, t n ) L ∞ (Ω) + I n h ϕ(•, t n ) -ϕ n h L ∞ (Ω) .
Bounding the first term in the right-hand side relies on standard results [4, Chap. IX, Lemme 1.1]. To bound the second one, we use (5.2) and for each x in V n h , we compute

I n h ϕ(x, t n ) -ϕ n h (x) = ϕ(x, t n ) -ϕ n-1 h (x n h ) = ϕ(x, t n ) -ϕ(x n h , t n-1 ) + ϕ(x n h , t n-1 ) -ϕ n-1 h (x n h ).
From the definition of x n h and the Lipschitz continuity of ϕ, we derive

|I n h ϕ(x, t n ) -ϕ n h (x)| ≤ cτ n (1 + u n-1 h L ∞ (Ω) d ) + ϕ(•, t n-1 ) -ϕ n-1 h L ∞ (Ω)
where c stands for the Lipschitz constant of ϕ on Ω×]0, T f [. Combining all this gives the desired estimate.

In view of the quantity T (ϕ), we also need an estimate for ϕ(•, t n )ϕ n h W 1,4 (Ω) . Proposition 5.3. Assume that, for a real number p ≥ 1, the solution ϕ of problem (3.3) satisfies ϕ ∈ C 0 (0, T f ; W 2,p (Ω)).

(5.5)

Then the following estimate holds for

1 ≤ n ≤ N ϕ(•, t n ) -ϕ n h W 1,p (Ω) ≤ C h n + Ch -1 n,min τ n (1 + u n-1 h L ∞ (Ω) d ) + ch -1 n,min ϕ(•, t n-1 ) -ϕ n-1 h L ∞ (Ω) , (5.6)
where h n,min stands for the minimum of the diameters h K , K ∈ T n h . Proof. We start once more from the triangle inequality

ϕ(•, t n ) -ϕ n h W 1,p (Ω) ≤ ϕ(•, t n ) -I n h ϕ(•, t n ) W 1,p (Ω) + I n h ϕ(•, t n ) -ϕ n h W 1,p (Ω) .
and use [4, Chap. IX, Lemme 1.2] to evaluate the first term. On the other hand, denoting by ψ x the Lagrange function associated with each node in V n h , switching to the reference element to evaluate the norm of this function in W 1,p (Ω) and using the fact, due the regularity of the family of triangulations, the support of each ψ x only intersects the support of a finite number of other ones, we derive

I n h ϕ(•, t n ) -ϕ n h p W 1,p (Ω) ≤ c x∈V n h |I n h ϕ(x, t n ) -ϕ n h (x)| p h d-p K x .
where K x is any element of T n h containing x. The maximmum of the |I n h ϕ(x, t n )ϕ n h (x)| has been evaluated in the previous proof. We derive from the regularity of the family of triangulations that

x∈V n h h d K x ≤ c meas(Ω),
which leads to the desired result.

We are now interested in the evaluation of the error issued from the discrete Navier-Stokes equations. We refer to [START_REF] Boukir | Razafindrakoto -A high-order characteristics/finite element method for the incompressible Navier-Stokes equations[END_REF]Section 4.3] for this evaluation in a simpler case. We begin with the terms involving ϕ n h and, for simplicity, we denote by ε n 1 and ε n 2 the right-hand sides of estimates (5.4) and (5.6), respectively.

≤ n ≤ N ρ(ϕ(., t n )) -ρ(ϕ n h ) L ∞ (Ω) + ν(ϕ(., t n )) -ν(ϕ n h ) L ∞ (Ω) ≤ ε n 1 ε , (5.7)
and

T (ϕ(., t n )) -T (ϕ n h ) L 2 (Ω) d×d ≤ c ε n 1 ε 2 + ε n 2 ε . (5.8)
Proof. The two estimates in (5.7) come from the fact that H ε is Lipschitz-continuous (the norm of its derivative is bounded by 1 ε ). To prove (5.8), we use the triangle inequality

T (ϕ(., t n )) -T (ϕ n h ) L 2 (Ω) d×d ≤ σ dH ε (ϕ(•, t n )) dϕ - dH ε (ϕ n h ) dϕ (I -∇ϕ ⊗ ∇ϕ)(•, t n ) L 2 (Ω) d×d + σ dH ε (ϕ n h ) dϕ (∇(ϕ(•, t n ) -ϕ n h ) ⊗ ∇ϕ(•, t n )) L 2 (Ω) d×d + σ dH ε (ϕ n h ) dϕ (∇ϕ(•, t n ) ⊗ ∇(ϕ(•, t n ) -ϕ n h )) L 2 (Ω) d×d + σ dH ε (ϕ n h ) dϕ (∇(ϕ(•, t n ) -ϕ n h )) ⊗ ∇(ϕ(•, t n ) -ϕ n h )) L 2 (Ω) d×d .
We observe that dH ε dϕ is also Lipschitz-continuous, with bounded Lipschitz constant (by π 2ε 2 ). We use Proposition 5.2 for the first term, Proposition 5.3 (with p = 4 for instance) for the last three terms, and obtain the desired result.

To go further, we subtract the first line of (4.9) from the first line of equation (3.7) at time t n , which gives for any v in

X n h , Ω ρ(ϕ)(x, t n ) ∂u ∂t + u • ∇u (x, t n ) -ρ(ϕ n h )(x) u n h -ũn-1 h τ n (x) v(x) dx - Ω (∇ • v)(x) p(x, t n ) -p n h (x) dx + Ω ν(ϕ)(x, t n )∇u(x, t n ) -ν(ϕ n h )(x)∇u n h (x) : ∇v(x) dx = Ω T (ϕ(x, t n )) -T (ϕ n h )(x) : ∇v(x) dx.
(5.9)

The key idea consists in inserting in this equation the mean value

M n (x) = 1 τ n t n t n-1 ∂u ∂t + u • ∇u (x, t) dt.
Indeed, it is equal to

M n (x) = u(x, t n ) -u(X(x, t n ; t n-1 )) τ n ,
thus easier to compare to the discrete one, and standard arguments yield that, if u belongs to H 2 (t n-1 , t n ; L 2 (Ω) d ) and C n denotes its norm in this space

( ∂u ∂t + u • ∇u (•, t n ) -M n L 2 (Ω) d ≤ C n τ n . (5.10)
From now on and for simplicity, we forget the dependency of the estimates with respect to ε (so that the righthand side of estimate (5.8) reads c ε n 2 ). Proposition 5.5. Assume that the solution ϕ of problem (3.3) satisfies (5.3) and that the solution (u, p) of problem

(3.6) -(3.7) is such that u ∈ W 1,∞ (Ω×]0, T f [) ∩ H 2 (0, T f ; L 2 (Ω) d ) ∩ L ∞ (0, T f ; H 2 (Ω) d ), p ∈ L ∞ (0, T f ; H 2 (Ω)).
(5.11)

Then, the following estimate holds for

1 ≤ n ≤ N τ -1 n u(•, t n ) -u n h L 2 (Ω) d + |u(•, t n ) -u n h | H 1 (Ω) d ≤ C (h n + τ n ) + c ε n 2 , (5.12)
where the constant C only depends on the regularity of ϕ, u and p.

Proof. Since many terms in (5.9) have already been evaluated, we present it in an abridged way. In (5.9), we add and subtract I n h u in the first terms and take v equal to I n h u(•, t n )-u n h . We observe from the second equations of (3.7) and (4.9) that the term involving the pressure can be written as

- Ω ∇ • (I n h u(•, t n ) -u n h ) (x) p(x, t n ) -p n h (x) dx = Ω ∇ • (u -I n h u)(•, t n ) (x) p(x, t n ) -p n h (x) dx + Ω ∇ • u n h (x)(p -I n h p)(x, t n ) dx,
and we integrate by parts the first term. All this leads to a bound of the quantity

τ -1 n I n h u(•, t n ) -u n h 2 L 2 (Ω) d + |I n h u(•, t n ) -u n h | 2 H 1 (Ω) d
by terms either evaluated in Lemma 5.4 or in (5.10) or involving interpolation error. We conclude by using triangle inequalities and interpolation estimates.

Remark 5.6. Owing to the inf-sup condition (4.12), an error estimate for the pressure p(•, t n )p n h can also be derived, we skip it for brevity. Unfortunately, the term ε n 2 contains error terms at time t n-1 , so we need an induction on n to conclude. We denote by h the maximum of the h n , 0 ≤ n ≤ N . Theorem 5.7. Assume that the solution (ϕ, u, p) of problem Then, the following estimates hold for

(3.3) -(3.6) -(3.7) is such that ϕ ∈ C 0 (0, T f ; W 2,∞ (Ω)), u ∈ W 1,∞ (Ω×]0, T f [) ∩ H 2 (0, T f ; L 2 (Ω) d ) ∩ C 0 (0, T f ; H 2 (Ω) d ), p ∈ L ∞ (0, T f ; H 2 (Ω)).
1 ≤ n ≤ N ϕ(•, t n ) -ϕ n h L ∞ (Ω) ≤ C n (h 2 + |τ |), (5.15
)

and u(•, t n ) -u n h H 1 (Ω) d ≤ C n (h + |τ |), (5.16)
where the constant C n depends on the regularity of ϕ, u and p, and also on n.

Proof. We proceed by induction on n. 1) For n = 0, estimates (5.15) and (5.16) are obvious since ϕ 0 h and u 0 h are the interpolates of ϕ 0 and u 0 , respectively. 2) Assume that (5.15) and (5.16) hold for n -1. We use the inverse inequality (with η > 0)

u n-1 h L ∞ (Ω) d ≤ w h L ∞ (Ω) d + h 1-d 2 -η n-1,min ( u(•, t n-1 ) -w h H 1 (Ω) d + u(•, t n-1 ) -u n-1 h H 1 (Ω) d ),
where w h is any approximation of u(•, t n-1 ) in X n-1 h

, and thus we derive from (5.14) and the induction hypothesis that u n-1 h L ∞ (Ω) d is bounded. Then, inserting the induction hypothesis in (5.4) and (5.12) gives (5.15) and (5.16), respectively.

Despite the technicity of the proofs, estimates (5.15) and (5.16) are fully optimal. But condition (5.14) is too restrictive: it involves both a hard Courant-Friedrichs-Lévy condition and the uniform regularity of the families of triangulations. Fortunately the same arguments as in the previous proof lead to the following result. (5.17)

Then, the following estimates hold for

1 ≤ n ≤ N ϕ(•, t n ) -ϕ n h L ∞ (Ω) ≤ C n (h 2 + |τ |), (5.18
)

and u(•, t n ) -u n h H 1 (Ω) d ≤ C n h -1 n,min (h 2 + |τ |), (5.19)
where the constant C n depends on the regularity of ϕ, u and p, and also on n. Assumption (5.17) is much more realistic and the convergence of the discretization holds for standard properties of the family of triangulations.

Remark 5.9. The fact that the constant C n in estimates (5.15), (5.16), (5.18) and (5.19) depends on n can be avoided either by weakly strengthening conditions (5.14) and (5.17) or working with specific family of parameters: for instance, it can be assumed that the quantity h 2 n + τ n decreases with n. Remark 5.10. We refer to [START_REF] Bernardi | Verfürth -A posteriori error analysis of the method of characteristics[END_REF] for a detailed a posteriori analysis of the method of characteristics. Since the same arguments hold for our discretization, we have decided not to perform the a posteriori analysis of problem (4.4) -(4.5) -(4.6) -(4.9).

Numerical results.

In the absence of analytical solutions, which are very hard to come by for this model, its validation is best accomplished with numerical benchmarking. In this section, we present three well-known test problems for the computation of two immiscible fluids with surface tension. All the examples are described in two dimensions and using FreeFem++ software, see [START_REF] Hecht | New development in FreeFem++[END_REF]. The level set method is applied to solve them with the ability to capture the changes of the topology correctly. So, the interface is captured by level set function which is defined as a signed distance function described in Section 2. Note that, all computations are performed using:

• the symmetric operator ∇ • (ν(∇u + ∇u t )) in the Navier-Stokes equation,

• the characteristics method with fixed time step, • the parameter g denotes the gravitational acceleration.

Rayleigh-Taylor instability.

The widely used test problem for numerical simulations for two-fluid flow is the Rayleigh-Taylor instability, see for instance [START_REF] Tryggvason | Numerical simulation of the Rayleigh-Taylor instability[END_REF]. When a layer of heavier fluid is placed on top of another lighter layer in a gravitational field with gravity pointing downward, the initial planar interface is unstable. Any disturbance will grow to produce spikes of heavier fluid moving downwards and bubbles of lighter fluid moving upwards. This is the so-called Rayleigh-Taylor instability. In this test case a heavy fluid is placed on the top of a light fluid and the initial position of the perturbed interface between two fluids is

ϕ(x) = tanh y -2 -0.1 cos(2π x) 0.01 √ 2
The computational domain is the rectangle ]0, a[×]0, 4a[ where a = 1 is the width. The density difference is normally represented by the Atwood number

At = ρ 2 -ρ 1 ρ 2 + ρ 1 ,
where ρ 2 > ρ 1 correspond to the heavier and lighter fluids, respectively. We also introduce the following Reynolds number

Re = ρ 1 a 3 2 g 1 2
ν .

According to [START_REF] Tryggvason | Numerical simulation of the Rayleigh-Taylor instability[END_REF], the governing equations are made dimensionless by using the following references: dx, and the rise velocity u 2 (x c , y c ; t) which is the velocity component in the direction opposite to the gravitational vector. For both test cases, we use the mesh adaptation in order to capture efficiently the interface between flows.

ρ ref = ρ 1 , a ref =
Test case 1 results: • the maximum rise velocity has an approximate magnitude of 0.24 and occurs between times t = 0.9 and t = 0.94 • the evolution of the center of mass can be described as a linear function of time and approaches y c = 1.09 toward the end of both simulations.

Test case 2 results:

As mentioned above, in the test case 2 , we increase the Eötvös number to Eo = 125. Figures 10 and11 describe the evolution in time of the shape deformation obtained in this work and in [START_REF] Hysing | Tobiska -Quantitative benchmark computations of two-dimensional bubble dynamics[END_REF]. So, we can observe that under t = 2.2 (the top panel of Figures 10 and11), the results are similar. Figure 12 shows all meshes during the computation. Nevertheless, the center of mass of both simulations is similar, see Figures 13 and14. Concerning the maximum rise velocity, the results in [START_REF] Hysing | Tobiska -Quantitative benchmark computations of two-dimensional bubble dynamics[END_REF] shows that there exist two velocity maxima, the first occurring at time t = 0.7332 with a magnitude of 0.2524 and the second one at t = 2.0705 with a slightly smaller magnitude of 0.2434, see the right side of Figure 12. Whereas in our work, Figure 13 (right side) shows two velocity maxima. The first is obtained at time ≈ 0.7 with a magnitude ≈ 0.25. But, the maximum of the second rise velocity is close to 0.24 at time t ≈ 1.7. Figure 16 shows the time evolution of the deformation of the bubble with respect to the tension surface coefficient σ. We observe that when σ increases from σ = 0 to σ = 20, the shape deformation of the bubble decreases and we remark that the shapes for the small σ = 0.75 and 1.5 have the same geometries. However, for σ = 5, 10 or 20, all results are different especialy when t = 8.

. 2 )Figure 1 .

 21 Figure 1. Two examples of domain Ω.

  7) is derived by exactly the same arguments as for the standard Navier-Stokes equations. So we refer to [31, Chap. III, Thms 3.1 & 3.2] and [13, Chap. V, Thms 1.4 & 1.5] for the next result. Proposition 3.3. Assume that the function ϕ belongs to L 4 (0, T

  ) problem (3.6) -(3.7) has at least a solution (u, p). Moreover, in dimension d = 2, this solution is unique, up to an additive constant on the pressure.Unfortunately, it seems difficult to prove, specially in dimension d = 3, that this solution satisfies the regularity properties required in Proposition 3.1 (even if they could be weakened). Up to our knowledge, the only existence result for the full problem (3.3) -(3.6) -(3.7) is due to Milcent [21, Chap. IV, Th. 1]; we only quote it and refer to [21, Chap. IV] and [9, Thm 2.2] for its proof (which relies on space and time regularization). Proposition 3.4.

Proposition 4 . 2 .

 42 For any datum ϕ 0 continuous on Ω, problem (4.4)(or (4.11)) -(4.5) -(4.6) has a unique solution ϕ n h in M n h . The well-posedness of problem (4.9) relies on standard arguments: Defining x n h by (4.4) or (4.11) as previously, we know ũn-1 h

Lemma 5 . 4 .

 54 If the assumptions of Propositions 5.2 and 5.3 are satisfied, the following estimates hold for 1

(5. 13 )

 13 For 1 ≤ n ≤ N , let us choose the parameters τ and h such that, h ≤ c h n,min , |τ | ≤ c h 2 n,min . (5.14)

Corollary 5 . 8 .

 58 Assume that the solution (ϕ, u, p) of problem (3.3) -(3.6) -(3.7) satisfies (5.13). For 1 ≤ n ≤ N , let us choose the parameters τ and h such that, for a positive constant η, h 2 + |τ | ≤ c h d 2 +η n,min .

  a and t ref = t √ a At g. A no-slip condition u = (0, 0) is enforced at the bottom and top walls while the first component of velocity u 1 = 0 is imposed on the two vertical sides, and we work with ε = 0.05. The interface shape at different dimensional times (for instance, t = 0s, 1s, 2s, 3s, 4s) are plotted in Figure2and Figure3. In the first one, the Reynolds number Re is equal to 100 and the Atwood number is equal to At = 0.3, while in the second one we take Re = 1000 and At = 0.5. During the early stages, the growth of the interface is slow and remains symmetrical. However, the characteristic mushroom shape emerges in the vicinity of the central vortex. Later, the heavier fluid falls continuously and the lighter fluid keep rising to form bubbles along the vertical side boundary, and the heavier fluid begins to roll up into two counter-rotating vortices, (e) t = 4s (d) t = 3s (a) t = 0s (b) t = 1s (c) t = 2s

Figure 2 .

 2 Figure 2. Rayleigh-Taylor instability. At = 0.3 and Re = 100: evolution of the interface shape in time

Figure 3 .

 3 Figure 3. Rayleigh-Taylor instability. At = 0.5 and Re = 1000: evolution of the interface shape in timeFinally, in order to study the effect of the Atwood number, we fixe the Reynolds number at Re = 1000 and we increase the Atwood number: At = 0.3, 0.35, 0.4, 0.45 and At = 0.5. It can easily be seen in Figure4or Figure5that with increasing Atwood number

Figure 4 .

 4 Figure 4. Rayleigh-Taylor instability. Re = 1000: evolution of the interface shape with respect to At

5 At = 0. 5 At 5 Figure 5 .

 5555 Figure 5. Rayleigh-Taylor instability. Re = 1000: evolution of the interface shape with respect to At

Figure 7 .

 7 Figure 7. Bubble shape at the final time (t = 3) for the test case 1. Left: Our result. Right: Hysing et al.

Figures 7

 7 Figures 7 shows both shapes obtained by our simulation in the left panel and by Hysing et al. in the right side. They have a similar form and we observe that the maximum of width and the height have almost the same values. Furthermore, Figure 8 confirms quantitatively that the rise velocity and the centre of mass are identical:• the maximum rise velocity has an approximate magnitude of 0.24 and occurs between times t = 0.9 and t = 0.94 • the evolution of the center of mass can be described as a linear function of time and approaches y c = 1.09 toward the end of both simulations.

Figure 8 . 1 Figure 9 .

 819 Figure 8. Our work: Center of mass (top) and rise velocity (bottom) for test case 1

  (a) t = 0.6 (b) t = 1.2 (c) t = 1.8 (d) t = 2.2 (e) t = 2.4 (f) t = 2.6 (g) t = 2.8 (h) t = 3

Figure 10 . 2 Figure 11 .

 10211 Figure 10. Our result: Typical time evolution of the for test case 2

= 3 Figure 12 .

 312 Figure 12. Our result: Typical time evolution of the adaptive mesh for test case 2

Figure 13 . 2 ii= 8 Figure 16 .

 132816 Figure 13. Our work: Center of mass (top) and rise velocity (bottom) for test case 2

Two-dimensional rising bubble.

The next application is a circular bubble rising in a viscous fluid. For this simulation, we use data from the numerical experiment of Hysing et al. [START_REF] Hysing | Tobiska -Quantitative benchmark computations of two-dimensional bubble dynamics[END_REF] where the bubble is initially circular with radius r = 0.25 and center coordinates (0.5, 0.5) placed at the bottom of a rectangular domain ]0, 1[×]0, 2[ with another fluid of higher density and viscosity. Figure 6 the initial configuration of this problem and the boundary conditions we use. The parameter ε is taken equal to 0.025. Buoyancy effects will make the drop move to the top of the domain and undergo some deformation. The result shape depends naturally on different physical parameters, for instance on Reynolds number Re = 2rρ 1 √ 2rg ν 1 and on the so-called Eötvös number which gives the ratio of gravitational forces to surface tension effects Eo = 4r 2 gρ 1 σ . The density and viscosity of the heavy fluid are ρ 1 , ν 1 whereas density and viscosity of the fluid occupied by the bubble are ρ 2 , ν 2 . In order to compare our results to the benchmark performed in [START_REF] Hysing | Tobiska -Quantitative benchmark computations of two-dimensional bubble dynamics[END_REF], we decide to realize two test cases: the first one considers a bubble with small density and viscosity ratios = 10, which undergoes moderate shape deformation and for the second test case, a bubble with a very low density and viscosity, where the ratios are 1000 (resp. 100) compared to that of the surrounding fluid.

• Test case 1: We present in this section a numerical simulation of the evolution of a bubble in a so-called jogged channel. The parameter ε is now taken equal to 0.05. The viscosities and densities of the external (resp. internal) fluid are ν 1 = 1, ν 2 = 0.1 and ρ 1 = 100, ρ 2 = 1.

(0, -2) (4, -4) (16, -4)

Figure 15. Boundary condition and initial configuration

This test is mentioned shortly by Chessa et al. in [START_REF] Chessa | Belytschko -An extended finite element method for two-phase fluids[END_REF], and we are going to try to complete this study by varying the surface tension coefficient σ and then plotting the shape deformation, the center of mass (x c , y c ) (as in the above example) and both velocity components u 1 (x c , y c ), u 2 (x c , y c ). The initial droplet shape has a largest diameter 2, as shown in Figure 15. During time, the droplet changes its shape and then translate from left to right in the jogged channel without gravitational effects. An the inlet (left wall), the Poiseuille flow u = u in = (2-y) (y+2) 4 , 0 is assigned and the natural Neumann condition

is imposed in the outlet (right wall). All boundary conditions are explained in Figure 15.

The parameter ε is fixed equal to 0.05. At each step, the mesh adaptation is used in order to capture the interface nicely and the time step is taken constant equal to 4 • 10 -2 .


Dennis Fox

Frédéric Darboux

come

Topographic controls on soil properties affecting post-fire erosion and sediment redistribution in a mixed forested-agricultural Mediterranean catchment

Topographic controls on soil properties affecting post-fire erosion and sediment redistribution in a mixed forested-agricultural Mediterranean catchment D.M. Fox (1), F. Darboux (2), P. Carrega (1)

(1) UMR 6012 CNRS, Department of Geography, University of Nice Sophia Antipolis, 98 Blvd. Ed. Herriot, BP 3209 06204 Nice cedex 3, France. (fox@unice.fr), (carrega@unice.fr), (2) INRA Orléans -Science du Sol, Av. de la Pomme de Pin, BP 20619, F-45166 Olivet Cedex, France. (Frederic.Darboux@orleans.inra.fr) Landuse organisation in Mediterranean France tends to follow a typical pattern where large flat coastal areas are densely urbanised and oriented toward summer tourism, inland alluvial plains are dedicated to vineyards for wine making, and hilly areas surrounding inland plains are occupied by forests and isolated housing. This spatial organisation reflects long term patterns of sediment redistribution where upslope geological substrates have provided sediments to downslope geomorphological units such as alluvial plains and colluvial deposits. Forest fires periodically destroy the forest canopy and litter layer and create conditions of aggravated erosion for a few years at a time. Typically, a forested area has probably undergone dozens of cycles of post-fire erosion and regrowth during the past millennia. The objective of this study was to analyse sediment dynamics for a burned forested area and compare these to erosion processes encountered in cultivated soils in the same catchment. The analysis was carried out in three steps. In the first, landuse was mapped using colour airphotos for a 230 km2 catchment in the "Massif des Maures" in SE France. For each of the two major landuse classes of interest (forested, vineyard), slope distributions were calculated based on a 25 m DEM. Similarly, channel slopes were calculated for stream sections flowing through the two landuse classes. This initial survey was complemented by two field studies. In the first, zones of river bank erosion and sediment deposition in channels were mapped over a total length of about 40 km. This was carried out to identify critical slope thresholds for major sediment deposition. In the second, a North facing slope was sampled for litter layer depth, surface stone cover, soil depth, and textu-ral and aggregate stability analyses for burned (August, 2003) and unburned forests and cultivated fields. The grainsize characteristics of slopes were compared to those of sediments trapped in stream channels. Almost three years after a forest fire in the Massif des Maures, the litter layer has only begun to recover and erosion processes are still active even though a substantial underbrush has developed and cork oak crowns have recovered. Runoff during post-fire rainfall events selectively removes the soil fraction finer than about 2 cm, leaving behind a stone lag on the surface. The steep slopes of the forested zone feed runoff directly into ephemeral streams, so there is virtually no deposition within the burned forested area. In addition, steep channel slopes eliminate the possibility of channel storage until the sediments have arrived well into the main river channel downstream. This indicates that delivery ratios are close to 100% for post-fire erosion in the burned Mediterranean catchments of SE France. The surface of the unburned forest is protected by a dense mat of litter, and the cultivated soils show some selective erosion depending on cultivation method. A comparison of erosion factors for forested and agricultural areas indicates that although both areas are sensitive to slope inclination, soil erodibility in the forested zone is dominated by the litter layer in the unburned forest, surface stone cover in the burned forest, and a combination of stone cover and aggregate stability in the agricultural soils.


Valentin Féray
email: valentin.feray@math.uzh.ch

I P Goulden
email: ipgoulden@uwaterloo.ca

Alain Lascoux

An edge-weighted hook formula for labelled trees

Keywords: hook formula, tree enumeration, combinatorial bijection, generating function

A number of hook formulas and hook summation formulas have previously appeared, involving various classes of trees. One of these classes of trees is rooted trees with labelled vertices, in which the labels increase along every chain from the root vertex to a leaf. In this paper we give a new hook summation formula for these (unordered increasing) trees, by introducing a new set of indeterminates indexed by pairs of vertices, that we call edge weights. This new result generalizes a previous result by Féray and Goulden, that arose in the context of representations of the symmetric group via the study of Kerov's character polynomials. Our proof is by means of a combinatorial bijection that is a generalization of the Prüfer code for labelled trees.

Introduction 1.Background

The classical hook formula of Frame, Robinson and Thrall [START_REF] Frame | The hook graphs of the symmetric group[END_REF]Theorem 1] gives the simple ratio

χ λ (Id |λ| ) = |λ|!
∈λ h( ) for the dimension χ λ (Id |λ| ) of the irreducible representation of the symmetric group associated with the Young diagram λ. Here |λ| is the number of boxes in the diagram and h( ) is the size of the hook attached to the box . This result is equivalent to an enumerative result, since it is also the number of labellings of the boxes of λ with the elements of N |λ| = {1, . . . , |λ|} (once each) so that the labels increase along each row, and down each column.

Many results that look similar have appeared since, and are commonly referred to as hook formulas. A number of these involve various classes of trees. Let us fix some terminology. A (unordered) tree is an acyclic connected graph. The vertexset (or label-set) of a tree T is denoted by V (T ). Rooted means that we distinguish a vertex; then each edge can be oriented towards the root and we call the head and tail of the edge father and son, respectively. We denote the father of a vertex v in a rooted tree T by f T (v), and set f T (v) = 0 when v is the root vertex. Then the rooted tree T is completely defined by giving f T (v) for all vertices v (in particular, unless specified differently, sons of a given vertex are not ordered). The descendants of a vertex are defined recursively as the sons and the descendants of the sons. If u is a descendant of v, then we say that v is an ancestor of u. The hook attached to the vertex v in the tree T , denoted by h T (v), is the set consisting of v and its descendants; the size of the hook h T (v) is denoted by h T (v). An increasing labelling of a rooted tree is a labelling of the vertices with distinct integers, so that the label of a son is always bigger than the label of its father; thus the root always gets the minimum label, and the labels increase along each branch from the root. An increasing tree is an increasing labelling of the rooted tree.

It is well-known that the number of ordered 1 increasing labellings of a given rooted tree is given by a formula that looks like Frame-Robinson-Thrall formula. Namely, D. Knuth [6,§5.1.4 Exer. 20] proved that the number L(T ) of ordered increasing labellings of a rooted tree T with vertex-set N |T | is given by ( 1)

L(T ) = |T |! v∈T h T (v)
, where |T | is the number of vertices of T . Another type of hook formula is a hook summation formula. For example, let B r denote the set of rooted binary trees with r vertices (as usual for binary trees, sons of a given vertex are ordered). There is a well-known one-to-one correspondence between increasing binary trees with vertex-set N r , and permutations of size r (see e.g. [7, p. 23-25]). Combining this with the rooted tree hook result [START_REF] Biane | Characters of symmetric groups and free cumulants[END_REF], and dividing by r!, yields the summation formula [START_REF] Borchardt | Über eine der Interpolation entsprechende Darstellung der Eliminations-Resultante[END_REF] T ∈Br v∈T

1 h T (v) = 1.
1 Here, ordered means that labellings More details on these hook formulas and some related works can be found in [START_REF] Féray | A multivariate hook formula for labelled trees[END_REF]. In this article, two of us gave a hook summation formula that involved unordered increasing trees, which means that the sons of a vertex are not ordered. For our summation formula, we use the following notation for falling factorials: (a) m = a(a -1) • • • (a -m + 1) for positive integers m, with (a) 0 = 1, and (a) m = 1/(a -m) -m for negative integers m. Let r ≥ 1 be an integer and x 1 , • • • , x r be formal variables. Let U r denote the set of unordered increasing trees with vertex-set N r , and for T ∈ U r , define a weight wt(T ) by

wt(T ) = r v=2 x fT (v) u∈hT (v) x u -h T (v) + 1 .
Then our hook summation formula [START_REF] Féray | A multivariate hook formula for labelled trees[END_REF] was given by (3)

T ∈Ur wt(T ) = x 1 • • • x r (x 1 + • • • + x r -1) r-2 .
Three proofs of this result were presented in [START_REF] Féray | A multivariate hook formula for labelled trees[END_REF]. One of these involved Kerov's character polynomials (see, e.g., [START_REF] Biane | Characters of symmetric groups and free cumulants[END_REF]), and thus gives a connection to the representation theory of the symmetric group, that does not seem related to the Frame-Robinson-Thrall formula.

We also proved that (3) specializes to a classical enumerative formula for Cayley trees. A Cayley tree is a tree with labelled vertices (so they are distinguishable) -these are not embedded in the plane, and there is no root vertex. Let C r denote the set of Cayley trees with vertex-set N r . Borchardt [START_REF] Borchardt | Über eine der Interpolation entsprechende Darstellung der Eliminations-Resultante[END_REF] and Cayley [START_REF] Cayley | A theorem on trees[END_REF] proved that, for r ≥ 1, (4)

T ∈Cr x dT (1) 1 • • • x dT (r) r = x 1 • • • x r (x 1 + • • • + x r ) r-2 ,
where d T (i) denotes the degree of the vertex i in the tree T . We proved in [START_REF] Féray | A multivariate hook formula for labelled trees[END_REF] that (3) specializes to (4) in the case x 1 , • • • , x r → ∞, that is for the highest degree terms in the x i . On the right-hand sides, this is straightforward, so the work here is on the left-hand sides, for which we constructed a combinatorial mapping between the sets U r and C r .

The main result

In this paper, we prove a new hook summation formula for unordered increasing trees. This formula is given in the following Theorem, which is our main result. This generalizes (3) by introducing a set of doubly indexed indeterminates that we will refer to as edge weights.

Theorem 1.1. Let r ≥ 2 be an integer and x i , i = 1, . . . , r, y i,j , 2 ≤ i ≤ j ≤ r be formal variables. For an unordered increasing tree T with vertex-set N r , define the weight to be

wt y (T ) = r v=2 x fT (v) u∈hT (v) y v,u . Then (5) T ∈Ur wt y (T ) = x 1 y r,r r-1 i=2 i j=1 x j y i,i + r j=i+1 x i y i,j .
Note that Theorem 1.1 specializes to (3) immediately, by the following substitution: y v,u = x u -1 for v < u, and y u,u = x u .

Our proof of Theorem 1.1 is by a combinatorial bijection. This involves a number of stages, and in our description, it will be convenient to identify the lefthand and right-hand sides of (5) separately, as ( 6)

L(x, y) = T ∈Ur wt y (T ), R(x, y) = x 1 y r,r r-1 i=2 i j=1 x j y i,i + r j=i+1 x i y i,j .
Of course, in these terms, our main result is equivalent to [START_REF] Stanley | Enumerative Combinatorics[END_REF] L(x, y) = R(x, y).

Outline of paper

In the remainder of this paper we give a combinatorial proof of our main result. This is carried out by defining a combinatorial mapping in Section 2 that we describe in terms of an operation on unordered increasing trees called splice. Then in Section 3 we prove a number of properties of our splice operation, enabling us to prove that the combinatorial mapping is a bijection. This directly proves [START_REF] Stanley | Enumerative Combinatorics[END_REF], and hence Theorem 1.1.

There is one intriguing aspect of our main result that we have been unable to resolve. Note that our proof of the main result in this paper is based on a bijection for U r , the set of unordered increasing trees. However, if we evaluate the righthand side of (5) at x i = 1 for all i, and y i,j = 1 for all i, j, then we obtain r r-2 . But as we have noted above, |C r | = r r-2 ([2], [START_REF] Cayley | A theorem on trees[END_REF]), which suggests that there should be a combinatorial proof of the main result based on a bijection for C r , the set of Cayley trees. We have been unable to find such a proof, and suggest it as a problem for others to resolve.

A combinatorial mapping

Dominating functions

For a set S of positive integers, let Π(S) denote the set of partitions of S into an unordered set of nonempty subsets. The subsets are called the blocks of the partition, and we denote the number of blocks of a partition π by |π|. If π has blocks π 1 , . . . , π k , then we let µ i = max π i , for i = 1, . . . , k, and we index the blocks so that

µ 1 < • • • < µ k .
For two sets S and S of positive integers, the function g : S → S is called dominating if g(i) ≥ i for all i ∈ S. For such a function g, we denote

wt g = i∈S y i,g(i) .
We say that a dominating function g : S → S (i.e., with S = S) is "on S". Consider the functional digraph of a dominating function g on S: the vertices are the elements of S, and the directed edges are given by (i, g(i)), i ∈ S. The vertex-sets of the connected components (ignoring the directions on edges) form a partition π ∈ Π(S), and we say that g has induced partition π. Let D(π) denote the set of all dominating functions on S with induced partition π, and let

D(π) = g∈D(π) wt g = g∈D(π) i∈S y i,g(i) .
For π ∈ Π(S), let E(π) denote the set of unordered increasing trees T on vertex-set S such that every block of π is a subchain of T . In other words, for every pair of elements i < j in the same block of π, i is an ancestor of j in T . For any unordered increasing tree T , let

(8) κ(T ) = i∈V (T ) x σi(T ) i
, where σ i (T ) denotes the number of sons of vertex i in T . Now we consider a restricted class of set partitions. If S is a set of positive integers containing 1, then Π 1 (S) is the set of partitions of S in which {1} is a block. In this case, necessarily π 1 = {1}, and µ 1 = 1. For such a partition π ∈ Π 1 (S), let

C(π) = N µ2 ו • •×N µ|π|-1 = {(c 2 , . . . , c |π|-1 ) : 1 ≤ c i ≤ µ i , i = 2, . . . , |π|-1},
(if |π| = 2, the set C(π) contains one element: the empty list). For c ∈ C(π), let [START_REF]π ∈ Π 1 (N r ), where π has k ≥ 2 blocks, and consider an arbitrary tree T ∈ E(π). Then T is a π-increasing tree in which the vertex-set consists of two or more blocks of π and contains 1, and from Proposition 3.6, T is reducible. From Lemma 3.8 (a)[END_REF] ω(c, π) =

x c2 • • • x c|π|-1 x µ2 • • • x µ|π| i∈S x i .
There is a close connection between dominating functions and the expressions L(x, y), R(x, y) defined in [START_REF] Knuth | The Art of Computer Programming[END_REF], given in the following result.

Proposition 2.1. For any integer r ≥ 2, one has:

(a) y 1,1 • L(x, y) = π∈Π1(Nr) D(π) T ∈E(π) κ(T ), (b) y 1,1 • R(x, y) = π∈Π1(Nr) D(π) c∈C(π) ω(c, π).
Proof. (a) For this equation, by definition,

L(x, y) = T ∈Ur
wt y (T ).

But wt y (T ) = κ(T ) g wt g , where the sum runs over functions g from {2, . . . , r} to {2, . . . , r} such that, for each u, its image g(u) lies in h T (u) (as T is an increasing tree, such functions are automatically dominating). We can extend such functions g to g : N r → N r by setting g(1) = 1. Note that wt g = y 1,1 wt g . The conditions g(1) = 1 and g(u) ∈ h T (u) for u ≥ 2 are equivalent to π ∈ Π 1 (N r ) and T ∈ E(π), where π is the partition induced by g. Therefore

y 1,1 L(x, y) = T ∈Ur κ(T )    π∈Π 1 (S) s.t. T ∈E(π) g∈D(π) wt g    = π∈Π1(S) D(π) T ∈E(π) κ(T ),
giving part (a) of the result. (b) For this equation, note that R(x, y) can be rewritten as

R(x, y) = x 1 x 2 • • • x r-1 y r,r r-1 i=2   y i,i x 1 + • • • + x i x i + r j=i+1 y i,j   .
Expanding the product in terms of dominating functions, we get

R(x, y) = x 1 x 2 • • • x r-1 y r,r g:{2,...,r-1}→{2,...,r} g dominating   wt g i : g(i)=i x 1 + • • • + x i x i   .
As above, we can extend g to g : N r → N r by setting g(1) = 1 and g(r) = r. Then wt g = y 1,1 y r,r wt g . Now let π denote the partition induced by g. An important remark is that the integers i = 1, r such that g(i) = g(i) = i are exactly the maxima of the blocks of π except for 1 and r, which are given by µ

2 , • • • , µ |π|-1 .
Note also that µ |π| = r. Hence

i : g(i)=i x 1 + • • • + x i x i = 1 x µ2 • • • x µ|π|-1 c∈C(π) x c2 • • • x c|π|-1 ,
and so we obtain

x 1 x 2 • • • x r-1 i : g(i)=i x 1 + • • • + x i x i = c∈C(π) ω(c, π).
Thus R(x, y) is given by

y 1,1 R(x, y) = g:Nr →Nr g(1)=1,g(r)=r g dominating wt g c∈C(π) ω(c, π).
Note that g dominating implies g(r) = r while the condition g(1) = 1 means that the induced partition π is in Π 1 (N r ). Hence, splitting the sum depending on the induced partition of g, we obtain part (b) of the result.

Comparing Proposition 2.1 with (7), we see that Theorem 1.1 is implied by a bijection (10)

ψ π : C(π) → E(π) : c → T,
with the weight-preserving property that κ(T ) = ω(c, π), for each π ∈ Π 1 (N r ).

We will find such a bijection ψ π .

An operation on rooted trees

A convenient construct for an unordered increasing tree T with vertex v is the v-decomposition of T , described as follows.

Definition 2.2. Let T 1 be an unordered increasing tree with root vertex a 1 , and let v be any vertex in T (v can be equal to a 1 ). Suppose that the unique maximal chain from a 1 to v is given by a

1 < a 2 < • • • < a k = v, k ≥ 1.
Now remove the edges on the chain in T from a 1 to a k . There are k components in the resulting graph, each of which is an unordered increasing tree, whose root vertex is on the chain from a 1 to a k . Let T (ai) be the component among these that is rooted at vertex a i , i = 1, . . . , k. Then the v-decomposition of T is the ordered list T (a1) , . . . , T (ak) .

An example of v-decomposition of a tree T is given in Figure 1.

The bijection ψ π will be constructed as the iteration of an elementary combinatorial operation on marked trees called splice, that we define next, in terms of v-decompositions. Definition 2.3. Suppose that T 1 and T 2 are two unordered increasing trees with disjoint vertex-sets, and let v 1 , v 2 be vertices in

T 1 , T 2 respectively, with v 1 > v 2 . Let the v 1 -decomposition of T 1 be T (a1) 1 , . . . , T (ak) 1
, and the v 2 -decomposition of for some unique ≥ 1 and

T 2 be T (b1) 2 , . . . , T (bm) 2 . Then, since a k = v 1 > v 2 = b m , we have b 1 < • • • < b β1 < a 1 < • • • < a α1 < b β1+1 < • • • < b β2 < a α1+1 < • • • < a α2 < • • • < b β -1+1 < • • • < b β < a α -1+1 < • • • < a α ,
1 ≤ α 1 < • • • < α = k, 0 ≤ β 1 < • • • < β = m.
Then define the splice of T 1 and T 2 , with splicing vertices v 1 and v 2 , denoted by

(11) R = splice(T 1 , v 1 ; T 2 , v 2 ),
to be the unordered increasing tree with v 1 -decomposition given by

(12) T (b1) 2 , . . . , T (bβ 1 ) 2 , T (a1) 1
, . . . , T

(aα 1 ) 1 , . . . , T (bβ -1 +1) 2 , . . . , T (bβ ) 2
,

T (aα -1 +1) 1 , . . . , T (aα ) 1
.

An example of the splice operation is given in Figure 2. Note in the construction of R above that the vertex-set of R is the (disjoint) union of the vertex-sets of T 1 and T 2 . Also, if β 1 ≥ 1, equivalent to b 1 < a 1 , then the root vertex of R is b 1 , the same as the root vertex of T 2 . However, if β 1 = 0, equivalent to a 1 < b 1 , then the root vertex of R is a 1 , the same as the root vertex of T 1 .

In the following result we record some simple but important properties of the splice operation.

Lemma 2.4. Suppose that R = splice(T 1 , v 1 ; T 2 , v 2 ).
(a) For = 1, 2, and any pair of vertices i and j in T , then i is an ancestor of j in T if and only if i is an ancestor of j in R. Consequently, given the sets V (T 1 ) and V (T 2 ), one can recover T 1 and T 2 from R. (b) For all vertices i in T 1 , we have

σ i (R) = σ i (T 1 ). For the vertex v 2 in T 2 ,
we have σ v2 (R) = σ v2 (T 2 ) + 1, but for all other vertices i in T 2 , we have

σ i (R) = σ i (T 2 ).
Proof. (a) For = 1, consider the v 1 -decomposition of T 1 :

T (a1) 1 , . . . , T (ak) 1
. Let I and J denote the indices such that i ∈ T (aI ) 1 and j ∈ T (aJ ) 1 .

• If I = J then i is an ancestor of j in T 1 if and only if i is an ancestor of j in T

(aI ) 1

. The same is true in R.

• If I = J, then i is an ancestor of j in T 1 if and only if I < J and i = a I .
The same is true in R.

In both cases, we see that i is an ancestor of j in T 1 if and only if i is an ancestor of j in R. This ends the proof of part (a) for = 1. The case = 2 is similar.

(b) This is immediate.

A candidate for our bijection.

In this section we describe a mapping ψ π that we claim is a suitable bijection for (10). To describe ψ π , consider a set partition π ∈ Π 1 (N r ) with |π| = k, so π has blocks π 1 = {1}, . . . , π k , and 1 = µ 1 < • • • < µ k where µ i is the largest element of π i , for i = 1, . . . , k. Recall that r ≥ 2 (and hence k ≥ 2). Consider also a (k -2)-tuple c = (c 2 , . . . , c k-1 ) ∈ C(π). We apply an iterative procedure in which we have a forest of unordered increasing trees on vertex-set N r at every stage, and we apply splice to reduce the number of components by one between successive stages.

Construction 2.5. Initially, at Stage 1, we have the forest with components τ 1 , . . . , τ k , where τ is the increasing chain whose vertices are the elements of the set π , = 1, . . . , k (so τ 1 consists of the single vertex 1). At every stage we also keep track of an integer ν in [r], with ν = 1 initially. Then, at Stage i, for i = 2, . . . , k -1, we input a forest with components τ 1 , τ i , τ i+1 , . . . , τ k , together with an integer ν, and create the following output: Iwahori-Heckesp Case 1. If c i is a vertex in τ 1 or τ i , then set τ 1 = splice(τ i , µ i ; τ 1 , ν), omit τ i , and set ν = c i ; Case 2. If c i is a vertex in τ j for some j > i, then set τ j = splice(τ i , µ i ; τ j , c i ), and omit τ i .

After completion of the above procedure, we are left with a pair of increasing rooted trees τ 1 and τ k , and an integer ν. Then finally, at Stage k, we let

(13) ψ π (c) = T, where T = splice(τ k , µ k ; τ 1 , ν).
Remark. During the construction, ν is always a vertex of τ 1 and moreover ν ≤ µ i after Stage i. Hence the splice splice(τ i , µ i ; τ 1 , ν) is well-defined in Case 1.

Similarly, since c i ≤ µ i , the splice in Case 2 is well-defined (we cannot have c i = µ i in Case 2, as this would imply that c i is a vertex of τ i ).

An example of the mapping ψ π is given in Figure 3 where r = 9 and k = 4. At Stage 2, we applied Case 1 because c 2 = 4 was a vertex of τ 2 , while, at Stage 3, we applied Case 2 because c 3 = 5 was a vertex of τ 4 . At each stage, the value of ν is recorded by an edge on τ 1 without child.

Proposition 2.6. Given r ≥ 2, a set partition π ∈ Π 1 (N r ), and c ∈ C(π), suppose that ψ π (c) = T is constructed as in (13) above. Then

T ∈ E(π),
and κ(T ) = ω(c, π).

Proof. By construction, it is clear that T is a tree with vertex-set N r . We have to check that it is indeed in E(π). Initially, in the iterative procedure for ψ π , we have components τ i , the chain consisting of the elements of the block π i of π, i = 1, . . . , k. The rooted tree T is constructed by applying the splice operation k -1 times, to join the initial components together in some order. The fact that T ∈ E(π) now follows immediately from Lemma 2.4(a). Also, initially, we have

κ(τ 1 ) • • • κ(τ k ) = 1 x µ2 • • • x µk r i=2 x i .
But for the terminating tree T , from Lemma 2.4(b), we have

κ(T ) = x 1 x c2 • • • x ck-1 κ(τ 1 ) • • • κ(τ k ),
and the fact that κ(T ) = ω(c, π) now follows immediately from (9) (in this case we have S = N r ).

Comparing Proposition 2.6 with (10), we see that the mapping ψ π above is indeed a candidate for a bijective proof of our main result. Theorem 2.7. For each r ≥ 2 and π ∈ Π 1 (N r ), the mapping

ψ π : C(π) → E(π)
is a bijection.

We will prove Theorem 2.7 in the next Section, by determining the inverse of ψ π . In our development, we will find it convenient to use terms that distinguish between the different ways in which "splice" is applied in Construction 2.5 -a splice that arises in Case 1 or in (13) (the final stage) is called an internal splice, whereas one that arises in Case 2 is called an external splice.

Inverting the combinatorial mapping and a bijective proof of the main result

The goal of this section is to construct the inverse of ψ π in order to show that it is a bijection. Throughout this Section, π is a fixed partition in Π 1 (N r ). As in the previous Section, the blocks of π are denoted by π 1 , π 2 , . . . , π k , and the maximum elements in these blocks are denoted by µ 1 , µ 2 , . . . , µ k . We assume as before that

µ 1 < µ 2 < • • • < µ k .
When it is convenient, we will also use the notation π x for the block of π containing the element x. This should not be confused with π i , which denotes the i-th part of π. Consider a set S that is a union of blocks of π (we shall say that S is πcompatible). In other words, S = ∪ i∈I π i , for some index set I. We denote by π |S the partition {π i , i ∈ I} of S. A tree T ∈ E(π |S ) for some π-compatible S is said to be π-increasing.

Dependence graphs and irreducibility

We begin by defining a directed graph associated with the v-decomposition of an unordered increasing tree. Definition 3.1. Consider a π-increasing tree T , with (π-compatible) vertex-set S = ∪ i∈I π i . For a vertex v in T , suppose that the v-decomposition of T is given by T (a1) , . . . , T (a ) , where a = v. Then the v-dependence graph of T , denoted by G v (T ), is a directed graph with the following vertices and directed edges:

• the vertex-set of G v (T ) is {π i , i ∈ I};
• for the directed edges of G v (T ), consider, if any, a maximum element µ i , i ∈ I, which is not contained in the chain a 1 , • • • , a of T . This vertex belongs (as a nonroot vertex) to T (aj) for some j = 1, . . . , . Then there is a directed edge from π i to π aj , for each such i ∈ I.

For example, fix π = {1}, {3, 4, 6}, {2, 7, 8}, {5, 9} , as in Figure 3. Then G 9 (T ), the 9-dependence graph of the tree T , obtained in the final stage of Figure 3, is drawn in Figure 4.

Note that the graph G v (T ) defined above depends strongly on the partition π, but since π is fixed throughout the section, we have omitted it from the notation. In the remainder of this paper, we will particularly consider v-dependence graphs for the case in which v is the vertex of maximum label. This motivates the following definition.

Definition 3.2. A π-increasing tree T with maximum vertex v is called irreducible if G v (T ) is connected. Otherwise, T is called reducible.
Remark. As we shall see later (Proposition 3.6), any π-increasing tree with vertex set [r] is reducible. Therefore, the notion of irreducibility is interesting only for trees with smaller vertex sets.

We now determine the form of G v (T ) when v is the maximum vertex of an irreducible tree T . Lemma 3.3. If a π-increasing tree T with maximum vertex v is irreducible, then the v-dependency graph G v (T ) is an indirected tree rooted at π v (i.e. all edges are directed towards the root).

Proof. Let the v-decomposition of T be given by T (a1) , . . . , T (a ) , where a = v, so v appears as the root vertex in T (a ) . Then v, which is the maximum element in the block π v , cannot appear as a nonroot vertex in any tree of the v-decomposition, and π v has outdegree 0 in G v (T ). But every other vertex has outdegree at most 1, and since G v (T ) is connected, it can only be a tree in which every other vertex has outdegree exactly 1 (in a connected graph, the number of vertices minus the number of edges is at most 1, and a difference of 1 occurs only for trees).

We now consider the effect on irreducibility of applying the splice operation. Lemma 3.4. Suppose that T 1 and T 2 are π-increasing trees with disjoint vertexsets, and that T 1 is irreducible. For i = 1, 2, let m i denote the maximum element in the vertex-set of T i . Let v 2 be a vertex in T 2 with v 2 < m 1 , and let

T = splice T 1 , m 1 ; T 2 , v 2 . (a) If m 1 > m 2 , then T is reducible. More precisely, (14) G m1 (T ) = G m1 (T 1 ) G v2 (T 2 ). (b) If m 1 < m 2 and T 2 is irreducible, then T is irreducible.
Proof. (a) From ( 12), the trees in the m 1 -decomposition of T are either trees in the m 1 -decomposition of T 1 or trees in the v 2 -decomposition of T 2 . Then equality (14) follows from Definition 3.1. In particular, G m1 (T ) is not connected, and since m 1 is the maximum vertex in T , we conclude that T is reducible.

(b) Suppose that the m 1 -decomposition of T 1 is T (a1) 1 , . . . , T (a ) 1
, and the

v 2 - decomposition of T 2 is T (b1) 2 , . . . , T (bn) 2
. We have a = m 1 , and since m 1 is the maximum vertex in T 1 , then T To describe the edges of G m2 (T ), we shall first consider the m 2 -decomposition of T . It is obtained as follows:

• start with the m 2 -decomposition of T 2 ;

• for each a w < b u , add the tree T (aw) 1

. Indeed, the vertex a w is in the chain from the root to m 2 in T , and the tree rooted at a w in the m 2 -decomposition of T is the same as in the v 1 -decomposition of T 1 ;

• for each a w > b u , add all vertices of T (aw) 1

to the tree rooted at b u . Indeed, a w and m 2 are in different branches below b u in the tree T (because m 2 is in T (bu) 2

). Hence, T (aw) 1 is entirely included in the tree rooted at b u in the m 2 -decomposition of T . This is illustrated in Figure 5. In this Figure, the m 1 -decomposition of T (which is by definition the union of the m 1 -decomposition of T 1 and the v 2 -decomposition of T 2 ) is represented with blue and red dashed lines. The m 2 -decomposition of T is drawn with plain black lines. Finally, we have used green dotted lines for the m 2 -decomposition of T 2 (it should be understood that the tree rooted at b u in this decomposition contains only the vertices of T 2 in the corresponding green dotted shape and of course no vertices of T 1 ).

From this, we can describe the edges of G m2 (T ): • start with the edges of G m1 (T 1 ) together with the edges of G m2 (T 2 );

T (a1) 1 T (b1) 2 T (bu) 2 m 1 m 2 T (bn) 2 T (a ) 1 v 2 b u
• for any maximum element µ i = m 1 that appears as a vertex in T

for any a w > b u , remove the edge from π i to π aw that appears in G m1 (T 1 ), and insert an edge from π i to π bu ; • finally, insert an edge from π m1 to π bu .

Recall that G m1 (T 1 ) and G m2 (T 2 ) are both connected by hypothesis, and we want to show that G m2 (T ) is connected. To do this, we will show that each vertex is in the connected component of π bu . For vertices in G m2 (T 2 ), this is obvious as G m2 (T ) contains all edges of G m2 (T 2 ). From Lemma 3.3, G m1 (T 1 ) is a directed tree of root π m1 . Then, for any vertex v in G m1 (T 1 ) either the path from v to m 1 is also in G m2 (T ) and then the edge from π m1 to π bu proves that v and π bu are in the same connected component, or this path is broken because one of its edge has been replaced by an edge to π bu . In this case, the same conclusion that v and π bu are in the same connected component holds. Thus G m2 (T ) is connected, and since m 2 > m 1 , m 2 is the maximum vertex in T , so we conclude that T is irreducible.

Example. We give illustrations of Lemma 3.4 in both cases (a) and (b): (a) Call T 1 and T 2 the trees τ 4 and τ 1 respectively of Stage 3 from Figure 3. Then m 1 = 9 and we choose v 2 = 4. The graph G 9 (T 1 ) has vertices {π 3 , π 4 } and an edge from π 3 to π 4 , while G 4 (T 2 ) has vertices {π 1 , π 2 } and an edge from π 2 to itself. Then T = splice T 1 , m 1 ; T 2 , v 2 is given in Stage 4 of Figure 3 and its dependence graph G 9 (T ), drawn in Figure 4, is indeed the disjoint union of G 9 (T 1 ) and G 4 (T 2 ).

T 1 =
(b) Set π = {1}, {3, 6}, {2, 7}, {9}, {4, 8, 10}, {5, 11} Consider the trees T 1 and T 2 from Figure 6 and choose v 2 = 5. The tree T = splice T 1 , 10; T 2 , 5 is also drawn in Figure 6. The corresponding dependence graphs are also given, showing that T is indeed irreducible. Note that, as explained in our proof, G 11 (T ) differs from the disjoint union of G 10 (T 1 ) and G 11 (T 2 ) as follows: the edge from π 4 to π 5 in G 10 (T 1 ) is replaced in G 11 (T ) by an edge from π 4 to π 6 ; moreover, a new edge from π 5 to π 6 has also been added. The graph G 11 (T ) is determined using the 11-decomposition of T , given in Figure 7.

The above result allows us now to classify the trees at every stage of Construction 2.5 by their irreducibility or reducibility. Proposition 3.5. In Construction 2.5:

• At Stage 1, the trees τ 1 , . . . , τ k are π-increasing and irreducible.

• For i = 2, . . . , k, in the input to Stage i, the trees τ i , . . . , τ k are π-increasing, irreducible, and contain µ i , . . . µ k , respectively; if τ 1 has been created by applying an internal splice at some previous stage, then τ 1 is π-increasing and reducible. • The tree ψ π (c) = T is π-increasing and reducible.

Proof. At Stage 1, each tree τ i consists of vertices of the single block π i of π, arranged as an increasing chain, so each τ i is π-increasing. Also, since G µi (τ i ) has only the single vertex π i , it is a connected graph, so τ i is irreducible.

This result for Stage 1 serves as the base case for an induction on the stages, in which the remaining results follow immediately from Lemma 2.4(a) (for πincreasing), and Lemma 3.4 (part (b) for irreducible, and (a) for reducible).

In particular, Proposition 3.5 establishes that all trees T created as images of our combinatorial mapping ψ π are reducible. Thus, in order for ψ π : C(π) → E(π) to be a bijection for any π ∈ Π 1 (N r ), when r ≥ 2, it is necessary that all trees T in E(π) are reducible. We prove that this is indeed the case in the following result. (note that T belongs to E(π) implies in particular that T is π-increasing). Proposition 3.6. Consider a π-increasing tree T whose vertex-set contains 1. We assume that the vertex-set of T is not reduced to {1} and denote its maximum element by M .

Then π 1 and π M are in different connected components of G M (T ). In particular, then T is reducible.

Proof. Recall that the partition π contains π 1 = {1} as a block, with maximum element µ 1 = 1. Now consider a π-increasing tree T containing 1 and its Mdecomposition T (a1) , . . . , T (a ) , with a = M (where M is the maximum vertex of T ). But vertex 1 is the root vertex of every tree T in E(π), so a 1 = 1. Thus µ 1 is the root vertex of T (a1) , and cannot appear as a nonroot vertex in any tree of the r-decomposition of T .

This implies that π 1 has outdegree 0 in the v-dependence graph G v (T ). But π M has also outdegree 0 (see the proof of Lemma 3.3), and other vertices (if any) have outdegree at most 1. Therefore, π 1 and π M are in different connected components of G M (T ).

Irreducibility and inverting the combinatorial mapping

In this section we prove that each application of splice in our combinatorial mapping ψ π can be uniquely reversed by considering only the irreducibility or reducibility of the trees involved.

We begin with a simple condition for when a π-increasing tree can be written as the splice of two subtrees. Lemma 3.7. Consider a π-increasing tree T , a π-compatible nonempty subset V 1 ⊂ V (T ), and a vertex v 1 ∈ V 1 . Then T can be written as

T = splice T 1 , v 1 ; T 2 , v 2
for some π-increasing trees T 1 and T 2 with vertex-sets V 1 and

V 2 = V (T ) \ V 1 and for some vertex v 2 ∈ V 2 if and only if V 1 = i∈I π i ,
where {π i , i ∈ I} is a union of vertex-sets of connected components of G v1 (T ).

In this case, T 1 , T 2 and v 2 are unique.

Proof. Given T , v 1 and V (T 1 ), we immediately have

V 2 = V (T ) \ V 1 .
Then it is clear from Definition 2.3 that v 2 is uniquely the first element of V 2 on the chain from v 1 to the root vertex of T (note that v 2 < v 1 since T is an increasing tree). Moreover, the trees T 1 and T 2 themselves are then uniquely determined, since we know their v 1 -decomposition and v 2 -decomposition, respectively. It only remains to determine conditions for V 1 . Let the v 1 -decomposition of T be given by T (a1) , . . . , T (a ) , where a = v 1 . Then from Definition 2.3, a necessary and sufficient condition is that V 1 (and V 2 ) are unions of the V (T (aj) ), j = 1, . . . , . Since, by hypothesis, V 1 is a union of blocks of π, this is equivalent to saying that π is a union of blocks of the partition

Π 1 := π ∨ V (T (aj) ), j = 1, . . . , .
It remains to see that this partition is nothing other than

Π 2 := {X c , c connected component of G v1 (T )}, where X c = πi∈Vc π i .
To do this, take two partitions π x and π y which contain elements x and y in the same set V (T (aj) ). We want to prove that π x and π y are in the same connected components of G v1 (T ).

Call u the root of T (aj) At least one of these elements, say x, is different from u. Then by definition, there is an edge from π x to π u in G v1 (T ). If y = u, there is an edge from π x to π y , and thus they are in the same connected component of G v1 (T ). If y = u, the same argument as above implies that there is also an edge from π y to π u , and one can also conclude that π x and π y lie in the same connected component of G v1 (T ). Hence Π 1 is finer than Π 2 .

Conversely, suppose that there is an edge from π s to π t in G v1 (T ). Then this means that µ s is a nonroot vertex in some T (aj) , with a j ∈ π t , which implies that there are elements of both π s and π t in the same subtree T (aj) . Hence, Π 2 is finer than Π 1 .

We conclude that Π 1 = Π 2 , which ends the proof of the Lemma.

In the next result, which is the key to inverting ψ π , we consider a π-increasing tree in which the vertex-set consists of two or more blocks of π. For such a tree with vertex-set S and maximum vertex M , we call m = max(S\π M ) the second maximum vertex. Lemma 3.8. Suppose that T is a π-increasing tree in which the vertex-set consists of two or more blocks of π, and let M and m be the maximum and second maximum vertices, respectively.

(a) If T is reducible, then it can be written uniquely as

T = splice(T 1 , M ; T 2 , t),
where T 1 and T 2 are π-increasing trees subject to:

• M is a vertex in T 1 , t is a vertex in T 2 , • T 1 is irreducible.
Moreover, if 1 is a vertex of T , then it automatically belongs to T 2 . (b) If T is irreducible, then it can be written uniquely as

T = splice(T 1 , m; T 2 , t),
where T 1 and T 2 are π-increasing trees subject to:

• m is a vertex in T 1 , t and M are vertices in T 2 , with t < m, • T 1 and T 2 are irreducible.

Proof. (a) From Lemma 3.7, the vertex-set V 1 of T 1 must correspond to a union of connected components of G M (T ). In addition, from Lemma 3.4 (a), G M (T 1 ) is the graph induced by G M (T ) on V 1 . Hence, if we want T 1 to be irreducible, that is G M (T 1 ) to be connected, then V 1 must correspond to a single connected component of G M (T ). Moreover, since we require M to be in T 1 , it must contain the block π M .

Finally, V 1 is uniquely the vertex-set of the connected component of G M (T ) containing π M (note that t < M for all vertices t in T 2 , since M is the maximum vertex in T ). The result follows immediately from Lemma 3.7. Now consider the maximum element in π u , that we will denote by µ u . Then µ u = u is a descendant of u in T , and cannot be contained in T (u) M , since that would create a loop in G M (T ), again contradicting Lemma 3.3. Thus µ u is contained in T (u) m (as a nonroot vertex), which implies that there is a loop at π u in G m (T ).

But clearly π m has outdegree 0 in G m (T ). Putting this together with the facts that there is a loop at π u in G m (T ), and that each vertex has outdegree at most 1, we see that π u and π m are contained in different components of G m (T ), so G m (T ) is not connected. Moreover, the edge from π M to π u in G m (T ) implies that π M and π u are in the same component of G m (T ) (note that we can have π u = π M ).

Then Lemma 3.7 implies that T can be written uniquely as T = splice(T 1 , m; T 2 , t),

where T 1 and T 2 are π-increasing trees, m is in T 1 , t and M are in T 2 , with t < m, and T 1 is irreducible: this is obtained by letting V (T 1 ) be the vertex-set of the connected component of G m (T ) that contains π m , which means that G m (T 1 ) is connected. But the elements of π M are contained in V (T 2 ), so m is the maximum vertex in T 1 , and so T 1 is irreducible. It remains to prove that T 2 is also irreducible. To do this, we look at the Mdecomposition of T 2 (since M is the maximum vertex in T 2 ), which is obtained from the M -decomposition of T as follows (see the proof of Lemma 3.4):

• delete the blocks T (ai) M , for which a i belongs to T 1 . Since T (u) m is in T 2 , this can happen only for i < j, that is for blocks before T M still rooted at u. In particular, if a block X of π is in V (T 1 ) and if there is an edge from Y to X in G M (T ), then Y is also in T 1 . This means that, if a vertex X is deleted when going from G M (T ) to G M (T 2 ), all vertices pointed to it are also deleted, and recursively. For vertices that are not deleted their outgoing edge is not modified. ν lies in τ 1 ). Recall that, in the input of Stage k -1, τ k-1 is always irreducible (Proposition 3.5). Thus, we have uniquely reversed Stage k -1 in Construction 2.5, which was an internal splice. Note that, after this step, 1 and ν still lie in τ 1 .

Case 2. µ k-1 lies in τ k . In this case, µ k-1 is the second maximum of τ k , and recall that τ k is irreducible by construction. From Lemma 3.8 (b), τ k can be written uniquely as splice(T 1 , µ k-1 ; T 2 , t)

with T 1 , T 2 irreducible, where M lies in T 2 and t < µ k-1 . Call τ k-1 = T 1 , τ k = T 2 and c k-1 = t (the value of ν is unchanged). Then c k-1 lies in τ k . Recall that, in the input of Stage k -1, τ k-1 and τ k are always irreducible and M lies in τ k (Proposition 3.5). Thus, we have uniquely reversed Stage k -1 in Construction 2.5, which was an external splice. Note that, since ν and τ 1 have not been changed, 1 and ν still lies in τ 1 . Now, all remaining stages k -2, . . . , 2 of Construction 2.5 can be uniquely reversed exactly as for Stage k -1 (Case 2 in general is "µ i lies in τ i+1 , . . . ,τ k "). After reversing Stage i, we have trees τ 1 , τ i , . . . , τ k , such that τ 1 contains 1, and ν and τ i , . . . , τ k are irreducible and contain µ i , . . . , µ k , respectively. From Proposition 3.9, we recover at the end the initial forest with components τ , = 1, . . . , k, where τ is the increasing chain consisting of the elements of the block π of π. Along the way, we recover uniquely the elements c i of the (k -2)-tuple c ∈ C(π). We conclude that c = ψ -1 π (T ), and that ψ π is a bijection.

3 2

 3 must be counted as different labellings.

Figure 1 :

 1 Figure 1: The v-decomposition of a tree S (with v = 14).

Figure 2 :

 2 Figure 2: The splice of two trees (S is given in Figure 1; splicing vertices are underlined).

9 TFigure 3 :

 93 Figure 3: Applying ψ π to c = (4, 5), when π has blocks π 1 = {1}, π 2 = {3, 4, 6}, π 3 = {2, 7, 8}, π 4 = {5, 9}.

1 Figure 4 :

 14 Figure 4: Example of the v-dependence graph of a tree

  of the single root vertex a . Now T 2 is an increasing tree, so we have b1 < • • • < b n = v 2 .Also, by hypothesis we have v 2 < m 1 < m 2 , and we conclude that m 2 is not contained in the chain b 1 , . . . , b n , which means that m 2 appears as a nonroot vertex in T (bu) 2for some 1 ≤ u ≤ n. In particular, if a w > b u , then a w and m 2 are in different branches below b u in the tree T , and T (aw) 1 is entirely included in the tree rooted at b u in the m 2 -decomposition of T . Now the vertices of G m2 (T ) consist of the vertices of G m1 (T 1 ) together with the vertices of G m2 (T 2 ).

Figure 5 :

 5 Figure 5: Illustration of the proof of Lemma 3.4

Figure 6 :

 6 Figure 6: Splicing irreducible trees together with m 1 < m 2 yields a new irreducible tree.

Figure 7 :

 7 Figure 7: The 11-decomposition of the tree T from Figure 6.

Figure 9 :

 9 Figure 9: The 10-decomposition and 10-dependence graph of the tree T from Figure 6.

M

  in the decomposition;• replace the block T

* VF is partially supported by ANR projet PSYCO and SNF grant "Dual combinatorics of Jack polynomials".

† IPG is supported by a Discovery Grant from NSERC.