Hervé V J Le Meur
email: herve.lemeur@math.u-psud.fr

Derivation of a viscous KP equation including surface tension, and related equations

Keywords: 74J15, 35Q30, 76M45, 35Q35 water waves, shallow water, Boussinesq system, viscosity, KdV equation, surface tension, KP

The aim of this article is to derive surface wave models in the presence of surface tension and viscosity. Using the Navier-Stokes equations with a free surface, flat bottom and surface tension, we derive the viscous 2D Boussinesq system with a weak transverse variation. The assumed transverse variation is on a larger scale than along the main propagation direction. This Boussinesq system is proved to be consistent with the Navier-Stokes equations. This system is only an intermediate result that enables us to derive the Kadomtsev-Petviashvili (KP) equation which is a 2D generalization of the KdV equation. In addition, we get the 1D KdV equation, and lastly the Boussinesq equation. All these equations are derived for non-vanishing initial conditions.

Introduction

Motivation

Understanding the evolution of water waves is a longstanding problem. For instance, Stokes already discussed a precise wave at least 130 years ago [START_REF] Stokes | Supplement to a paper on the theory of oscillatory waves[END_REF] and his name has been given to this wave. Such a wave was still discussed more than one century later [START_REF] Hasimoto | Nonlinear Modulation of Gravity Waves[END_REF] ! One of the striking phenomena is that there exist surface waves that travel almost without modification. Such a behavior is an impetus to investigate the dynamics of these waves.

The first motivation is to understand how they may either disperse or dissipate and yet not vanish. It was discovered that there is a regime in which nonlinearity compensates dispersion (with no dissipation).

Another motivation is more technical. Should a wave keep its shape for a long time without vanishing, it could be used in telecommunication without requiring repeaters. Such a new technology could drastically reduce the cost of telecommunications.

Beyond mathematical and industrial considerations, mathematicians and physicists model the water waves with a fluid flow in a channel of (often) flat bottom, without meniscus (appart in [START_REF] Mei | The damping of surface gravity waves in a bounded liquid[END_REF]). Since the Euler or Navier-Stokes equations are rather prone to faithfully model such a fluid, one could be satisfied with either of these models. Yet, the Direct Numerical Simulation is too expensive and asymptotic models are required. In such models, one assumes a regime of dimensionless parameters' smallness and makes expansions of the equations to derive a simplified model with less fields and less dimensions.

Literature

The techniques are very different depending on whether one chooses Euler or Navier-Stokes as a model for the fluid.

On the one hand, for an inviscid flow (Euler), one assumes very often irrotationality. Then one may use a potential function and a Dirichlet to Neumann operator in the Zakharov-Craig-Sulem formulation. Numerous articles use this formulation and a review of the mathematical proofs has been published [START_REF] Lannes | The Water Waves Problem: Mathematical Analysis and Asymptotics[END_REF]. Some authors use a velocity-pressure formulation, but they still assume the irrotationnality assumption [START_REF] Iguchi | A mathematical justification of the forced Korteweg-de Vries equation for capillary-gravity waves[END_REF]. The derivation of Kadomtsev-Petviashvili (KP) assumes a predominant propagation direction along x and a weakly transverse propagation along y. It is initially done for an inviscid fluid in [START_REF] Kadomtsev | On the stability of solitary waves in weakly dispersive media[END_REF] and the justification of this approximation is done in [START_REF] Lannes | The Water Waves Problem: Mathematical Analysis and Asymptotics[END_REF] (subsection 7.2). In [START_REF] Lannes | Consistency of the KP Approximation[END_REF], D. Lannes proves that the sum of one wave propagating to the right and one to the left, both obeying a KP equation converges (when ε tends to zero) on [0, T 0 /ε] to a function consistent with the Boussinesq system, meaning that it is a solution, should it be only up to O(ε 2 ).

Various difficulties must be considered with these surface wave models. Firstly two scales of the transverse velocity are choosen in the literature. Either one assumes the transverse velocity to be O [START_REF] Boussinesq | Lois de l'extinction d'une houle simple en haute mer[END_REF] or O( √ ε). This problem was raised in Remark 3 of [START_REF] Lannes | Weakly transverse Boussinesq systems and the Kadomtsev-Petviashvili approximation[END_REF] and here in Remark 3. It leads to different Boussinesq systems.

In [START_REF] Ming | Large time well-posedness of the threedimensional capillary-gravity waves in the long wave regime[END_REF], the authors prove that a dimensionless water wave system in an infinite strip under the influence of gravity and surface tension has a unique solution on [0, T /ε]. More precisely, if the initial solution is sufficiently regular and small (O( √ ε)), there exists a solution, on this time interval [0, T /ε], that will remain small (O( √ ε)). In [START_REF] Ming | Long-wave approximation to the 3-D capillarygravity waves[END_REF], the same authors prove that on the same time interval, these solutions can be accurately approximated by sums of solutions of two decoupled KP equations. It enhances the results of [START_REF] Lannes | Weakly transverse Boussinesq systems and the Kadomtsev-Petviashvili approximation[END_REF] by taking surface tension and a variable bottom into consideration.

The Boussinesq equation is a generalization of the KdV equation for waves moving both to the right and to the left. From the Boussinesq system, instead of changing to a frame moving to the right as is done for KdV, we remain in the same coordinate system and get higher order equation (perturbed wave equation), as is done in [START_REF] Johnson | A modern introduction to the mathematical theory of water waves[END_REF] (p. 216-219).

On the other hand, when one assumes the fluid to be viscous, one may not use the potential function. Moreover the number of boundary conditions changes then. Viscosity has been considered in water waves since 1895 by Boussinesq [START_REF] Boussinesq | Lois de l'extinction d'une houle simple en haute mer[END_REF]. The dynamic of viscous water waves on finite depth was investigated more recently by Kakutani and Matsuuchi [START_REF] Kakutani | Effect of viscosity of long gravity waves[END_REF]. They derived the viscous KdV equation from vanishing initial conditions but made the error of using a Fourier transform in time while the problem is of Cauchy type.

Liu and Orfila [START_REF] Liu | Viscous effects on transient long-wave propagation[END_REF], and the coauthors of P.L.-F. Liu in subsequent articles, investigated the viscous Boussinessq system with vanishing initial condition, and validated their KdV equation by some experiments. For instance they experimentally proved the reverse flow in the boundary layer, which is predicted by theory.

Although it was not done in Sobolev spaces, the derivation of the viscous Boussinesq system (1D and isotropic 2D) and 1D KdV equation with a non-vanishing initial condition was done in [START_REF] Le Meur | Derivation of a viscous Boussinesq system for surface water waves[END_REF] without surface tension. It was partially the goal of [START_REF] Le Meur | Derivation of a viscous Boussinesq system for surface water waves[END_REF] not to rely on the irrotationnality assumption so as to derive the Boussinesq system and KdV equation for a viscous fluid. This could be achieved thanks to the fact that one of the equations contains u z = O(ε) which is the zeroth order of the irrotationnality assumption. In [START_REF] Le Meur | Derivation of a viscous Boussinesq system for surface water waves[END_REF], a more detailed bibliography was given and the reason why all the preceding articles did not derive the same KdV equation is explained.

In this article, we intend to cross the results for a viscous fluid with some results concerning inviscid fluids (KP, Boussinesq equation, surface tension).

The outline of this article is as follows. In section 2, we present the equations and state the viscous Boussinesq system in 2D with surface tension and a predominant propagation direction. This result is needed to get the new results on viscous 1D KdV equations in section 3, the viscous KP equation (2D generalization of KdV with a weak transverse propagation) derived in section 4, and the viscous Boussinesq equation derived in section 5. Finally, we justify that our Boussinesq system is consistent with the Navier-Stokes equations in Section 6. All these derivations are done for a viscous fluid with surface tension and non-vanishing initial conditions.

2D geometry with viscosity and surface tension

Below, we define the geometry, write the equations, and then make these equations dimensionless (subsection 2.1). These equations model a propagation predominantly in the x-direction. This will enable us to make a linear analysis with the dispersion relation and an asymptotic of the phase velocity in subsection 2.2. A discussion of the relevance of viscosity and surface tension will be given in subsection 2.3. We will then use the results of [START_REF] Le Meur | Derivation of a viscous Boussinesq system for surface water waves[END_REF], where the 1D viscous isotropic Boussinesq system and KdV equation are derived, so as to state the weakly transverse 2D viscous Boussinesq system (subsection 2.4). Our present derivation includes surface tension, viscosity and a propagation mainly along x and weakly along y. Since the already published results are very close, most proofs are omitted.

Navier-Stokes equations

Let ũ = (ũ, ṽ, w) be the velocity of a fluid in a 3-D domain Ω = {(x, ỹ, z) / (x, ỹ) ∈ R 2 , z ∈ (-h, η(x, ỹ, t))}. So we assume the bottom is flat and the free surface is characterized by z = η(x, ỹ, t) with η(x, ỹ, t) > -h (the bottom does not get dry). Let p be the pressure and D[ũ] the symmetric part of the velocity gradient. The dimensionless domain is illustrated in Fig. 1. We also denote ρ the density of the fluid, ν the viscosity of the fluid, g the

                 ρ ∂ ũ ∂ t + ũ. ∇ũ -ν ∆ũ + ∇p = -ρgk in Ω, div ũ = 0 in Ω, -pI + 2ν D[ũ] .n = -p atm n + σ R n on z = η(x, t), ηt + ũη x + ṽ ηỹ -w = 0 on z = η(x, t), ũ = 0 on z = -h, (1)
where we write the second order tensors and the vectors with bold letters. The differentiated functions are denoted either ∂f /∂x, f x , or ∂ x f . The surface tension term includes R given by geometrical computations:

1 R = (1 + η2 ỹ )η xx + (1 + η2 x) ηỹỹ -2η x ηỹ ηxỹ (1 + η2 x + η2 ỹ ) 3/2 .
Of course, we need to add an initial condition and conditions at infinity. So as to get dimensionless fields and variables, we need to choose a characteristic horizontal length l which is the wavelength along the propagation direction, a characteristic vertical length h which is the water's height, and the amplitude A of the propagating perturbation. Moreover, let U, V, W, P be the characteristic horizontal velocity along x, horizontal velocity along y, vertical velocity and pressure respectively. Since the zeroth order equation is a wave equation of dispersion relation

ω = c 0 k x (1 + k 2 y /2/k 2 x ) + o(k 2 y
), one may infer that a good scaling for the y direction is if the extra term is of the same order of magnitude as ε. So k y = O(ε 1/2 ) and the interesting scale for ỹ is O(ε -1/2 ) more than for the x direction. We may then define:

c 0 = √ gh, ε = A h , β = h 2 l 2 , U = εc 0 , V = εc 0 , W = √ εc 0 , P = ρgA, Re = ρc 0 h ν , Bo = ρgl 2 σ ,
where c 0 is the phase velocity. Then, one may make the fields dimensionless and unscaled:

ũ = Uu, ṽ = V v, w = W w, p = patm -ρgz + P p, η = Aη,
and the variables too:

x = lx, ỹ = ε -1/2 ly, z = h(z -1), t = t l/c 0 . (2)
We also make the Boussinesq approximation (β and ε of the same order of magnitude) and even β = ε. Notice that we take the scale used by [START_REF] Lannes | Weakly transverse Boussinesq systems and the Kadomtsev-Petviashvili approximation[END_REF] (V = εc 0 ) and not the one taken by [START_REF] Johnson | A modern introduction to the mathematical theory of water waves[END_REF] (V = ε 3/2 c 0 ). The difference will be highlighted in Remark 3. With these definitions, the new system with the new fields, variables and outward unit normal, still denoted n, writes in the new domain Ω t = {(x, y, z), (x, y) ∈ R 2 , z ∈ (0, 1 + εη(x, y, t))} (see Figure 1):

                                             u t + εuu x + ε 3/2 vu y + wu z - √ ε Re u xx + εu yy + u zz ε + p x = 0 in Ω t , v t + εuv x + ε 3/2 vv y + wv z - √ ε Re v xx + εv yy + v zz ε + p y √ ε = 0 in Ω t , w t + εuw x + ε 3/2 vw y + ww z - √ ε Re (w xx + εw yy + w zz ε ) + p z = 0 in Ω t , u x + √ εv y + w z /ε = 0 in Ω t , (η -p)n + 1 Re    2 √ εu x (εu y + √ εv x ) u z + w x (εu y + √ εv x ) 2εv y (v z + √ εw y ) u z + w x (v z + √ εw y ) 2 w z √ ε    .n = 1 Bo (1 + ε 4 η 2 y )η xx + ε(1 + ε 3 η 2 x )η yy -2ε 4 η x η y η xy (1 + ε 3 η 2 x + ε 4 η 2 y ) 3/2
.n on z = 1 + εη,

η t + εuη x + ε 3/2 vη y -w/ε = 0 on z = 1 + εη, u = 0 on z = 0.
(3) Simple computations give the outward non-unit normal n = (-εη x , -εη y , 1).

Linear theory

The whole subsection below is a straightforward modification of the case of a viscous fluid without surface tension, fully studied in [START_REF] Le Meur | Derivation of a viscous Boussinesq system for surface water waves[END_REF].

Dispersion relation

We are looking for small fields. So we linearize the system (3). The dynamic condition at the free boundary can then be stated more explicitely. We get:

                                   u t - √ ε Re (u xx + εu yy + u zz /ε) + p x = 0 in R 2 × [0, 1], v t - √ ε Re (v xx + εv yy + v zz /ε) + √ εp y = 0 in R 2 × [0, 1], w t - √ ε Re (w xx + εw yy + w zz /ε) + p z = 0 in R 2 × [0, 1], u x + √ εv y + w z /ε = 0 in R 2 × [0, 1], (u z + w x )/Re = 0 on z = 1, (v z + √ εw y )/Re = 0 on z = 1, (η -p) + 2w z Re √ ε = 1 Bo (η xx + εη yy ) on z = 1, η t -w/ε = 0 on z = 1, u = 0 on z = 0. (4)
In order to eliminate the pressure, we differentiate (4) 1 with respect to z and (4) 3 with respect to x and compute their difference so as to simplify p xz . Symmetrically, we differentiate (4) 2 with respect to z and (4) 3 with respect to y and compute their difference (with a √ ε coefficient). Pressure is eliminated in two equations on three:

u zt -w xt - √ ε Re ((u xxz -w xxx ) + ε(u yyz -w yyx ) + (u zzz -w xzz )/ε) = 0, v zt - √ εw yt - √ ε Re ((v xxz - √ εw xxy ) + ε(v yyz - √ εw yyy ) + (v zzz - √ εw yzz )/ε) = 0.
So as to eliminate u and v thanks to the incompressibility (4) 4 from the previous system, we differentiate the first equation with respect to x and the second with respect to y.

After some simplifications, a good combination of the two gives (thanks to (4) 4 ):

(∂ 2 z + ε∂ 2 x + ε 2 ∂ 2 y )(-Re √ ε∂ t + (∂ 2 z + ε∂ 2 x + ε 2 ∂ 2 y ))w = 0. (5)
Let w be of the form A(z) exp i(k x x + k y y -ωt) with non-negative k x , k y and a (complex) pulsation ω. We can define a complex parameter with non-negative real part, similar to the one used by [START_REF] Kakutani | Effect of viscosity of long gravity waves[END_REF]:

µ 2 = εk 2 -iωRe √ ε, (6)
where

k 2 = k 2 x + εk 2 y .
Thanks to this notation, the solutions of ( 5) are such that

A(z) = C 1 cosh √ εk(z -1) + C 2 sinh √ εk(z -1) + C 3 cosh µ(z -1) + C 4 sinh µ(z -1). ( 7
)
The surface tension did not appear yet because it is only in the boundary conditions.

Up to now we have eliminated u, v and p only in the volumic equations. We still have to use the boundary conditions of (4) to find the conditions on the remaining field w.

Easy computations on the boundary conditions of (4), similar to the ones done in [START_REF] Le Meur | Derivation of a viscous Boussinesq system for surface water waves[END_REF], give:

w z (0) = 0, w(0) = 0, w zz (1) -εw xx (1) -ε 2 w yy (1) = 0, w xx + εw yy (1) -w ztt (1) + 3 √ ε Re (w xxzt + εw yyzt ) + 1 Re √ ε w zzzt (1) = 1 Bo (∂ 2 x + ε∂ 2 y ) 2 w. (8)
The solutions [START_REF] Kakutani | Effect of viscosity of long gravity waves[END_REF] satisfies a homogeneous linear system in the constants

C 1 , C 2 , C 3 , C 4 .
Its matrix is:

    √ εk sinh ( √ εk) - √ εk cosh ( √ εk) µ sinh µ -µ cosh µ cosh ( √ εk) -sinh ( √ εk) cosh µ -sinh µ 2k 2 ε 0 µ 2 + εk 2 0 -k 2 -k 4 Bo √ εω 2 k + 2iεωk 3 Re -k 2 -k 4 Bo 2µ √ εiωk 2 Re     , (9)
where

k 2 = k 2 x + εk 2 y .
It suffices to compute its determinant to get the dispersion relation:

4εk 2 µ(εk 2 + µ 2 ) + 4µk 3 ε 3/2 (µ sinh (k √ ε) sinh µ -k √ ε cosh (k √ ε) cosh µ) -(εk 2 + µ 2 ) 2 (µ cosh (k √ ε) cosh µ -k √ ε sinh (k √ ε) sinh µ) -(k + k 3 /Bo) √ εRe 2 (µ sinh (k √ ε) cosh µ -k √ ε cosh (k √ ε) sinh µ) = 0. ( 10
)
This relation is identical to the one of [START_REF] Le Meur | Derivation of a viscous Boussinesq system for surface water waves[END_REF] appart from the surface tension term and the new definition of k for the 3D geometry (k 2 = k 2 x + εk 2 y ). Our process of nondimensionnalizing makes a difference between x and z. This is why [START_REF] Kakutani | Effect of viscosity of long gravity waves[END_REF] have k terms, and we have k √ ε terms instead.

Asymptotic of the phase velocity (very large Re)

In this subsection, we state the following Proposition concerning the phase velocity: 1) (and c bounded away from 0), Bo bounded away from 0, if there exists a complex phase velocity c = ω/k solution of [START_REF] Lannes | Consistency of the KP Approximation[END_REF], then it is such that:

Proposition 1. Under the assumptions k √ ε Re c → +∞, k = O(1), ε → 0, Re → +∞, c = O(
c = 1 + k 2 Bo tanh (k √ ε) k √ ε - e iπ/4 Re -1/2 (k √ ε) 1/4 2 tanh 3/4 (k √ ε) + o(ε -1/4 Re -1/2 ). (11)
Moreover, the decay rate in our finite-depth geometry does not depend on the surface tension but only on the viscosity. It is:

Im(ω) = Im(kc) = -1 2 √ 2 k 5/4 ε 1/8 √ Re tanh 3/4 (k √ ε) + o(ε -1/4 Re -1/2 ). (12)
Notice that our geometry is not infinite. So our formula may be different from Boussinesq's and Lamb's.

Appart from the surface tension terms and the change of notation of k, the proof is not different from the one of Proposition 1 in [START_REF] Le Meur | Derivation of a viscous Boussinesq system for surface water waves[END_REF].

To what extent is surface tension relevant ?

Relevance of surface tension and viscosity

As was explained in the previous subsection, surface tension yields only a real term in the phase velocity. So it does not give dissipation, but only dispersion. On the contrary, viscosity influences both the real and the imaginary part of the phase velocity. We choose two different fluids (water and mercury) in order to discuss the relevance of surface tension and viscosity. Their physical parameters are listed in Table 1. Furthermore, so as to compare 1/Bo and ε, and under the Boussinesq approximation, one may compute the ratio

σ(Nm -1 ) ρ(kg m -3 ) ν(Pa s) 1/
1 Bo /ε = σ ρgh 2 .
The critical value of h, for which this ratio is 1, is h app = σ/(ρg). This other critical value is either 2.7 mm for water or 1.9 mm for mercury. In other words if h > h app , then the surface tension term is irrelevant and even less important than gravity measured by ε in the equation. In most current experiments, h is more than 0.1m and so, surface tension should (globally) not be taken into account. As D. Lannes says: "for coastal waves of characteristic length L x = 10m, capillary effects represent only 0.0003% of the gravity effects" ([11] section 9.1.2).

Surface tension is irrelevant for flows with h 0 more than 1 cm. Yet, even if the initial flow satisfies this, it may or may not remain true. For instance, in case of flows with ripples created by the wind or in case of wave breaking, very short wavelengths appear. See section 9.1.2 of [START_REF] Lannes | The Water Waves Problem: Mathematical Analysis and Asymptotics[END_REF] for some references.

In the case of a flow in a shallow channel (some mm height), and for mercury, a depression solitary wave was predicted (as early as the initial article of Korteweg and de Vries [START_REF] Korteweg | On the change of form of long waves advancing in a rectangular canal, and on a new type of long stationnary waves[END_REF]) and observed in [START_REF] Falcon | Observation of depression solitary waves on a thin fluid layer[END_REF]. In this experiment, A = 0.064 mm, h = 2.12 mm. As a consequence, Re=2547 and viscosity should be taken into account.

In the following, we assume 1/Bo between O(ε) (l ∼ 1 mm) and O(ε 2 ) (l ∼ 1 cm).

Viscous Boussinesq System

Let us first discuss the regime binding ε, Re, Bo. The complex phase velocity [START_REF] Lannes | The Water Waves Problem: Mathematical Analysis and Asymptotics[END_REF] contains gravitational and viscous terms that we want to compare so as to find the regime at which their variations are of the same order of magnitude. The first term is the gravitational term (

(1 + k 2 /Bo) tanh (k √ ε)/(k √ ε))
which may be expanded when ε tends to zero:

1 -k 2 ε(1 -3/(εBo))/6 + O(ε 2 ) + O(Bo -2 ).
The second is purely viscous and can be expanded:

- √ 2(1 + i)(4 √ k) -1 (Re √ ε) -1/2 + o(Re √ ε) -1/2 ). So if 1/Bo = O(ε)
or less but 1/(εBo) ∼ 1/3, the variations of c on the gravitational and on the viscous effects are of the same order of magnitude when ε(1 -3/(εBo)) and (Re √ ε) -1/2 are of the same order. In this regime of very large Re, studied hereafter, the dependence of Re on ε is such that (to simplify, we assume here

1/Bo = o(ε)): Re ∼ ε -5/2 . ( 13
)
This would be wrong if 1/(εBo) ∼ 1/3 as is well-known for non-viscous fluids. In that case, we would need to go further in the expansion. But since it was assumed above that 1/Bo is between O(ε 2 ) and O(ε), we are driven to choose the regime ( 13) and exclude 1/(εBo) ∼ 1/3.

Our main purpose here is to state an asymptotic system of reduced size from the global Navier-Stokes equations in the whole moving domain with surface tension and a predominant propagation direction. Proposition 2. Let x = (x, y) ∈ R 2 . Let η(x, t) be the free boundary's height. Let (u b,0 (x, γ), v b,0 (x, γ)) for γ ∈ (0, +∞) (resp. (u u,0 (x, z), v u,0 (x, z)) for z ∈ (0, 1+εη(x, t))) be the initial horizontal velocity in the boundary layer (resp. in the upper part of the domain). Let f (x, γ) := u b,0 (x, γ)-u u,0 (x, z = 0) and g(x, γ) := v b,0 (x, γ)-v u,0 (x, z = 0). If for any x ∈ R 2 , γ → f and γ → g are uniformly continuous in γ and f x , g y ∈ L 1 γ (R + ), then the solution of the Navier-Stokes equation with this given initial condition satisfies the weakly transverse Boussinesq system:

u t + η x - η xxx Bo + εuu x + ε 3/2 vu y -εη xtt (z 2 -1) 2 = O(ε 2 ) + O(ε/Bo), v t + √ εη y - √ εη xxy Bo + εuv x + ε 3/2 vv y -εη ytt (z 2 -1) 2 = O(ε 2 ) + O(ε/Bo), η t + u x (x, z, t) -ε 2 η xxt (z 2 -1 3 ) + ε(uη) x + √ εv y + ε 3 2 (vη) y -ε √ πR (u x + √ εv y ) * 1 √ t + 2ε √ π +∞ γ ′′ =0 div (u b,0 (x, γ ′′ ) -u u,0 (x, z = 0)) √ ε(v b,0 (x, γ ′′ ) -v u,0 (x, z = 0)) √ R 4t γ ′′ γ ′ =0 e -γ ′2 dγ ′ dγ ′′ = o(ε), (14)
where the convolution, denoted with * , is in time, the parameters ε, and R have been defined and z ∈ (0, 1 + εη(x, t)).

With the change of notation k 2 = k 2

x + k 2 y , the whole proof is the very same as in [START_REF] Le Meur | Derivation of a viscous Boussinesq system for surface water waves[END_REF] and is omitted.

Had we proceeded in the same way as [START_REF] Kakutani | Effect of viscosity of long gravity waves[END_REF], we would have distinguished two subdomains: the upper part (z > ε) where viscosity can be neglected, and the lower part (0 < z < ε) which is a boundary layer at the bottom and where viscosity must be taken into account.

The resolution in each part would have enabled to match the solution on the common boundary. Instead, we assume overlapping domains.

Remark 3. As is stressed at the non-dimensionnalizing step, there are two choices for the scale of the transverse horizontal velocity. Either one takes V = ε 3/2 c 0 (see Johnson's book [START_REF] Johnson | A modern introduction to the mathematical theory of water waves[END_REF] for instance or many others), or one uses V = εc 0 (see the article of Lannes-Saut [START_REF] Lannes | Weakly transverse Boussinesq systems and the Kadomtsev-Petviashvili approximation[END_REF] for instance). Firstly one must notice that since this field is not eliminated (as w is), its scale is relevant. Moreover the latter scaling is more general than the former one because it does not assume v to be O(

√ ε). The main difference is that v Lannes-Saut = √ εv Johnson ,
where one implicitely assumes that the fields v Lannes-Saut and v Johnson are O(1). To be clear, Lannes-Saut's choice amounts to assume v = O(1) although it can be written √ ε∂ y ψ ! If we plug this v = √ ε∂ y ψ in ( 14), we get the same transverse Boussinesq system as [START_REF] Johnson | A modern introduction to the mathematical theory of water waves[END_REF] (p. 210) because ε 3/2 v Lannes-Saut transforms into ε 2 v Johnson and go to O(ε 2 ). So the non-dimensionnalizing process of Lannes-Saut is weaker than the one done by many other authors.

Lannes and Saut [START_REF] Lannes | Weakly transverse Boussinesq systems and the Kadomtsev-Petviashvili approximation[END_REF] give one more argument by exhibiting an (odd) solution to the linearized system (in their Remark 3) such that v = √ ε∂ y ψ. It is of size O(1) and not

O( √ ε): u = 0, v = f ′ (y - √ εt) -f ′ (y - √ εt), η = f ′ (y - √ εt) + f ′ (y - √ εt).
We want to notice here that the potential representation of the velocity might have induced the idea that if the potential is physical, then the deduced field (the velocity) is such that (u, v) = (∂ x φ, √ ε∂ y φ) and so, that v = O( √ ε). Had we started from the velocity pressure formulation, this assumption would not be natural. Yet, such an attempt of explanation does not match reality. Johnson [START_REF] Johnson | A modern introduction to the mathematical theory of water waves[END_REF] uses the velocity pression representation as we do and assumes v = O( √ ε), while Lannes and Saut [START_REF] Lannes | Weakly transverse Boussinesq systems and the Kadomtsev-Petviashvili approximation[END_REF] use the potential representation and do not make this assumption ! Remark 4. It was already noticed in [START_REF] Le Meur | Derivation of a viscous Boussinesq system for surface water waves[END_REF] that the double integral term in ( 14) is new and surprising because of its dependence on the initial condition. If we assume an initial flow of Euler type (so that u b,0 -u u,0 = 0), the double integral vanishes. But is it physical ? In other words, does an initial (inviscid) flow in the boundary layer (where Navier-Stokes applies) establish (as a Navier-Stokes flow) fast or not ? Can we assume an initial Euler flow and a Navier-Stokes evolution PDE without loss of generality ?

We claim that the answer is negative for at least two reasons. On the one hand, the characteristic time for the viscous effects to appear is roughly T N SE = ρh 2 0 /ν or T N SE = ρl 2 /ν. Then, its ratio with the characteristic time of the inviscid gravity flow (l/c 0 ) is either Re √ ε = ε -2 or Re/ √ ε = ε -3 respectively. Whatever the chosen characteristic time, this ratio is large and the flow in the boundary layer does not establish fast enough.

On the other hand, as is argued in [START_REF] Le Meur | Derivation of a viscous Boussinesq system for surface water waves[END_REF], the value of this integral for a typical exponential flow in the boundary layer can be computed. It tends to zero only like 1/ √ t as is classical for Navier-Stokes flows. As a consequence, one may assume that it goes from zero (Euler) to exponential (boundary layer) within the time (1/ √ t) which is large with respect to the characteristic time of the flow.

Viscous KdV with surface tension

Starting from [START_REF] Mei | The damping of surface gravity waves in a bounded liquid[END_REF], one may assume the fields do not depend on y, and v = 0. So the flow is purely one dimensional (one direction). One may see that the zeroth order of this Boussinesq system is the wave equation. Then one knows that there are two waves propagating in each direction. If we look only for the waves that propagate to the right, one may make a change of variables suggested by the zeroth order equation:

ξ = x -t, τ = εt.
Every term can easily be converted in these new coordinates appart from the convolution (one derivative and a half integration), the new surface tension term and the double integral. The first was treated in a very clean way in [START_REF] Le Meur | Derivation of a viscous Boussinesq system for surface water waves[END_REF]. The second is new, but very simple. The third one was treated in [START_REF] Le Meur | Derivation of a viscous Boussinesq system for surface water waves[END_REF], but its treatment is improved below. After a change of variable, this term writes:

+ 2ε √ π +∞ γ ′′ =0 u b,0 x (ξ + τ ε , γ ′′ ) -u u,0 x (ξ + τ ε , z = 0) × √ R ε 4τ γ ′′ γ ′ =0 e -γ ′2 dγ ′ dγ ′′ (15)
The article [START_REF] Le Meur | Derivation of a viscous Boussinesq system for surface water waves[END_REF] argued that if the initial horizontal velocity is localized and τ not too small, then u b,0 (ξ + τ /ε, γ ′′ ) -u u,0 (ξ + τ /ε, z = 0) will be negligible in comparison with ε. More precisely, if ξ → u b,0 (ξ, γ ′′ ) -u u,0 (ξ, z = 0) tends to zero when ξ tends to +∞ at least as ξ -1 . Then the term ( 15) is at least O(ε 2 ).

One may add one more argument if τ is still not too small. The upper bound of the inner integral contains a √ ε. So mainly large γ ′′ (larger than 1/ √ ε) will be relevant. But for those γ ′′ , the u b,0 (ξ + τ /ε, γ ′′ ) -u u,0 (ξ + τ /ε, z = 0) is small uniformly in ξ + τ /ε because of the matching condition on u 0 between the boundary layer and the upper part. So ( 15) is small for two different reasons and we can state the viscous KdV equation with surface tension in the following Proposition. Proposition 5. If the initial flow is localized, the KdV change of variables applied to the 1D version of the system (14) leads to

2η τ + 3η ηξ + 1 3 - 1 εBo ηξξξ - 1 √ πR τ /ε ξ ′ =0 ηξ (ξ + ξ ′ , τ ) √ ξ ′ dξ ′ = o(1), (16)
for not too small times τ , where we set Re= R ε -5/2 .

We do not replace the upper bound of the convolution (τ /ε) by +∞ for two reasons. On the one hand we do not know how fast this integral converges when τ /ε tends to +∞. On the other hand the τ term reminds us that this integral on ξ ′ is a mixture of time and space.

Viscous KP equation with surface tension

Below, we derive the viscous KP equation, then discuss the zero-mass (in ξ) constraint.

Derivation

We recall that we set Re= R ε -5/2 , and 1/Bo = O(ε) but 1/(εBo) = 1/3 + o(1). With these assumptions, the change of variables ξ = x -t, y = y, τ = εt, and the assumption that the initial flow is localized, the system ( 14) writes (since ε/Bo = O(ε 2 )):

-u ξ (ξ, y, z, τ ) + η ξ + ε(u τ + uu ξ -η ξξξ (z 2 -1)/2 - η ξξξ εBo + √ εvu y ) = O(ε 2 ), -v ξ + √ εη y + ε(v τ + uv ξ - √ εη yξξ (z 2 -1)/2 - η ξξy √ εBo + √ εvv y ) = O(ε 2 ), -η ξ + u ξ + ε(η τ + 1 2 η ξξξ (z 2 -1/3) + (uη) ξ + √ ε(vη) y ) + √ εv y - ε √ πR τ /ε ξ ′ =0 u ξ (ξ + ξ ′ , y, z, τ ) + √ εv y (ξ + ξ ′ , y, z, τ ) √ ξ ′ dξ ′ = O(ε 2 ). (17)
In the above system, we did not write the double integral term from [START_REF] Ming | Large time well-posedness of the threedimensional capillary-gravity waves in the long wave regime[END_REF]. Indeed, it was justified in the previous section that this term can be neglected. Moreover, one might prove that the transverse velocity is such that v(z) = v(z ′ ) + O(ε) in a similar way to Lemma 11 of [START_REF] Le Meur | Derivation of a viscous Boussinesq system for surface water waves[END_REF] (roughly, one differentiate ( 14) 2 with respect to z, then integrate with respect to time t). So as to go further, since, at the first order u ξ = η ξ , one may justify that u τ = η τ + o(1) as in the inviscid case, and u = η + o(1). One may then add the first and third equation:

2εη τ + 3εηη ξ + ε 3 - 1 Bo η ξξξ + ε 3/2 vu y + ε 3/2 (vη) y + √ εv y - ε √ πR τ /ε ξ ′ =0 η ξ (ξ + ξ ′ , y, τ ) + √ εv y (ξ + ξ ′ , y, z, τ ) √ ξ ′ dξ ′ = O(ε 2 ). (18)
In order to eliminate v, one must differentiate with respect to ξ to use the corresponding (17) 2 . Our choice of non-dimensionalization forces us to manage extra terms, but we get the classical KP equation with a viscous term:

2η τ + 3ηη ξ + 1 3 - 1 εBo η ξξξ ξ + η yy - 1 √ πR τ /ε ξ ′ =0 η ξξ (ξ + ξ ′ , y, τ ) √ ξ ′ dξ ′ = - √ ε (2vη ξy + v y η ξ + v τ y ) + O(ε) + O(1/Bo),
where v satisfies (17) 2 . This completes the proof of the following Proposition.

Proposition 6. If the initial flow is localized, in the Boussinesq approximation and if

1/Bo = O(ε) but 1/Bo = ε/3 + o(ε)
, the surface waves of a viscous fluid propagating predominantly along x and weakly along y satisfy the viscous KP equation:

2η τ + 3ηη ξ + 1 3 - 1 εBo η ξξξ ξ + η yy - 1 √ πR τ /ε ξ ′ =0 η ξξ (ξ + ξ ′ , y, τ ) √ ξ ′ dξ ′ = O( √ ε), (19)
or the equivalent mixed system

2η τ + 3ηη ξ + 1 3 - 1 εBo η ξξξ - 1 √ πR τ /ε ξ ′ =0 η ξ (ξ + ξ ′ , y, τ ) √ ξ ′ dξ ′ + v y = O( √ ε), -v ξ + η y = O( √ ε). (20)
One must keep in mind that our choice of scale triggers the √ ε (even inside the convolution) in ( 17). This choice is already discussed in Remark 3. Had we made the same choice of scale for v as Johnson [START_REF] Johnson | A modern introduction to the mathematical theory of water waves[END_REF] and many others, we would have replaced v by √ εv J in [START_REF] Molinet | Remarks on the mass constraint for KP-type equations[END_REF]. The equivalent ( 18) could be simplified by ε and it would write :

2η τ + 3ηη ξ + 1 3 - 1 εBo η ξξξ + v Jy - 1 √ πR τ /ε ξ ′ =0 η ξ (ξ + ξ ′ , y, τ ) √ ξ ′ dξ ′ = O(ε),
and v J would satisfy (17) 2 modified:

-v Jξ + η y + ε v Jτ + uv Jξ -η yξξ (z 2 -1)/2 - η ξξy εBo = O(ε 3/2 ).
Oddly, this would lead to the same KP equation as ( 19) but up to the order O(ε) and not

O( √ ε). During the proof, we exhibit the O( √ ε) terms which would be O(ε) if v = √ εv J
although they are not kept in the Proposition. So there is no contradiction and the accuracy of KP is tied to the property that v = O( √ ε) or not.

The zero-mass constraint

As is very well discussed in [START_REF] Molinet | Remarks on the mass constraint for KP-type equations[END_REF], the usual KP equation is often written with an operator

∂ -1 ξ ∂ 2
y . Yet such an operator assumes the solution is differentiated (in ξ) from a function that tends to 0 when ξ → ±∞. It is proved in [START_REF] Molinet | Remarks on the mass constraint for KP-type equations[END_REF] that although it is not obvious, this assumption is right for inviscid fluids. Their proof is done in two steps.

In the first step, the linear KP equation is solved with inverse Fourier transforms. For more general equations, but only dispersive, like

u t -Lu ξ + ∂ -1 ξ ∂ y y u = 0 (the symbol of L is ε | X | α for α > 1/2), the fundamental solution writes G(t, ξ, y) = F -1 (X,y ′ )→(ξ,y) e it(εX|X| α -y ′2 /X) ,
where ε = ±1 depends on the KPI or KPII equation (or on the sign of 1/3 -1/(εBo)).

For the usual KPI and KPII, α = 2.

Thanks to the Lebesgue's Dominated convergence theorem and inventive changes of variables, it is proved that the fundamental solution is regular, and is differentiated from a more regular function

A (which is C(R 2 ) L ∞ (R 2 ) C 1 ξ (R 2 )) such that A * u 0 → 0 when | ξ |→ +∞.
In a second step, the nonlinearity is treated with the Duhamel formulation. The solution is then the convolution of the group of the linear equation with the nonlinear term u u ξ and all what is already proved for the linear equation extends to the nonlinear one.

One may wonder whether such an interesting result can extend to the viscous case. Indeed, the proof is even simpler and the details are left to the interested reader. The Fourier transform of the dissipative term writes

-1 √ πR F ξ→X +∞ 0 η ξ (ξ + ξ ′ , y, τ ) √ ξ ′ dξ ′ = -1 √ 2π 2 R R ξ +∞ ξ ′ =0 η ξ (ξ + ξ ′ ) √ ξ ′ e -iξX dξ ′ dξ = -i X √ πR η(X) +∞ ξ ′ =0 e iξ ′ X | ξ ′ | dξ ′ .
By a simple change of variable, one may see that the last integral is indeed only a function of the sign of X (denoted sgn(X)) times 1/ | X |. This function of sgn(X) can be computed :

+∞ ξ=0 e i sgn(X)ξ | ξ | dξ = π 2 (1 + i sgn(X)),
and the Fourier corresponding term reads:

-i η(X) X 2R | X | (1 + i sgn(X)).
As a consequence, one could study the fundamental solution G N S (t, ξ, y) = F -1 (X,y ′ )→(ξ,y) e it(εX|X| α -y ′2 /X+(1+i sgn(X))X(2R|X|) -1/2 ) .

The only real part (non-dispersive) inside the exponential is it × i sgn(X)X(2R | X | ) -1/2 = -t | X | 1/2 (2R) -1/2 which sign is compatible with the dissipation and ensures convergence of all the integrals. So one do not even need to use the Lebesgue dominated convergence theorem and the proof is simpler than in [START_REF] Molinet | Remarks on the mass constraint for KP-type equations[END_REF].

Viscous Boussinesq equation with surface tension

The Boussinesq equation is a second order equation that takes into account waves going both to the right and to the left. Since the derivation is straightforward, we only sketch it, following Johnson [START_REF] Johnson | A modern introduction to the mathematical theory of water waves[END_REF] (p. 216-219) and state the Proposition 7. Johnson [START_REF] Johnson | A modern introduction to the mathematical theory of water waves[END_REF] uses a different scaling, but only along the transverse direction which is not used here.

Starting from the Boussinesq system [START_REF] Mei | The damping of surface gravity waves in a bounded liquid[END_REF] written only in 1D (with x, z, t and not x, y, z, t), one may differentiate [START_REF] Mei | The damping of surface gravity waves in a bounded liquid[END_REF] 3 with respect to t and (14) 1 with respect to x. The difference of these two equations writes after some easy computations:

η tt -η xx -ε u 2 + η 2 2 xx - ε 3 η xxtt + η xxxx Bo + ε √ πR η tt * 1 √ t + ε √ πR p.v. 1 √ t u u,0 x (x, z = 0) - +∞ γ ′ =0 u b,0 xγ (x, γ ′ ) e -R 4t γ ′2 dγ ′ = O(ε 2 ), ( 21
)
where the convolution is in time and p.v. denotes the principal value as defined in the theory of distributions. The p.v. stems from an integration by parts in which one has a boundary integral

ε √ πRt f 0 x (x, γ ′ )e -Rγ ′2 4t +∞ γ ′ =0
.

The u u,0 and u b,0 are the initial values of u u and u b . They must be provided. Appart from the two viscous terms and the term generated by surface tension, this equation is identical to (3.41) of [START_REF] Johnson | A modern introduction to the mathematical theory of water waves[END_REF]. Since we concentrate on x = O(1), we may write u =

x

-∞ u x = - x -∞ η t (x ′ , t)dx ′ + O(ε).
The above equation is then the Eulerian form of the Boussinesq equation rewritten in (23).

Let us now come back to the Lagrangian coordinates and change of fields. In [START_REF] Johnson | A modern introduction to the mathematical theory of water waves[END_REF], the author proposes X = x+ε x -∞ η(x ′ , t)dx ′ , but we consider the following change of variable easier to justify:

X = x -ε t 0 u(x, t ′ )dt ′ H(X, t) = H(x -ε t 0 u(x, t ′ )dt ′ , t) = η(x, t) -εη 2 (x, t). (22)
The function H is assumed to be defined on X in the initial domain. Using these definitions, it is straightforward to rewrite (21) under the form of a Lagrangian Boussinesq equation with surface tension and viscosity taken into account. Only the term containing the initial horizontal velocities needs some insight. It is the last one in (21) and both u u,0 (x, z = 0) and u b,0 xγ (x, γ ′ ) can be considered as depending on x or on X since an ε appears in front of the whole term and the difference between x and X is O(ε). The following Proposition states the two forms of the Boussinesq equation in 1D.

Proposition 7. Let a flow of a viscous fluid with surface tension in a 2D channel with a flat bottom like the one depicted in Figure 1. Let u u,0 (x, z = 0) denote the initial horizontal velocity in the upper part (respectively u b,0 (x, γ) in the boundary layer). We assume the Boussinesq approximation and 1/Bo = O(ε). Under these assumptions, the surface waves obey the viscous Boussinesq equation in the Eulerian form:

η tt -η xx -ε u 2 + η 2 /2 xx - ε 3 - 1 Bo η xxxx + ε √ πR η tt * 1 √ t + ε √ πR p.v. 1 √ t -(η t ) t=0 - +∞ γ ′ =0 u b,0 xγ (x, γ ′ ) e -Rγ ′2 4t dγ ′ = O(ε 2 ), ( 23
)
where the convolution is in time and p.v. denotes the principal value as defined in the theory of distributions. The Lagrangian form (with (22)) writes:

H tt (X, t) -H XX - 3 2 ε(H 2 ) XX - ε 3 - 1 Bo H XXXX + ε √ πR H tt * 1 √ t + ε √ πR p.v. 1 √ t -(H t (X, t)) t=0 - +∞ γ ′ =0 u b,0 xγ (x, γ ′ ) e -R 4t γ ′2 dγ ′ = O(ε 2 ), (24)
where the initial velocities defined on x can be evaluated either on x or on X.

Notice that the formulae (23,24) also have a formulation in terms of u u,0 (x, z = 0)u b,0 (x, γ) which is physically meaningful.

6 Is the Boussinesq system consistent with Navier-Stokes ?

The weakly transverse Boussinesq system is stated in Proposition 2. We intend to justify here that it is consistent with the initial Navier-Stokes equations (3). What is consistency ? In [START_REF] Le Meur | Derivation of a viscous Boussinesq system for surface water waves[END_REF], the author justifies, but does not prove in any functionnal space, that if there is a solution to Navier-Stokes equations, with a free boundary and a flat bottom, then, under the Bousinesq approximation, the solution satisfies the viscous isotropic Boussinesq system. Such a necessary result deserves to be completed by a sufficient one which is consistency. This is the goal of the present section.

Notice that the consistency as defined by Lannes in [START_REF] Lannes | The Water Waves Problem: Mathematical Analysis and Asymptotics[END_REF] (Definition 5.1) is much more rigorous. While we use only expansions, Lannes asked whether a solution to the asymptotic model (assuming it exists), in the asymptotic regime, satisfied the initial complete model on a convenient time-interval and in convenient norms, up to a given power of the small parameter. Lannes proved in [START_REF] Lannes | Consistency of the KP Approximation[END_REF] that KP solution (with his scaling) is consistent with the Boussinesq system.

We prove below the following Proposition in the 1D case.

Proposition 8. Let u u (x, z, t) be the horizontal velocity in the upper part of a flow in a geometry as in Figure 1, and η(x, t) be the height of the free boundary above the fluid, for x = (x, y) ∈ R 2 , z ∈ (ε, 1 + εη(x, t)) and nonnegative time t. Let (u u,0 (x, z), v b,0 (x, z)) denote the initial horizontal velocity in the upper part (respectively (u b,0 (x, γ), v b,0 (x, γ)) in the boundary layer with γ = z/ε).

If (u u , v u , η) satisfies the non-isotropic Boussinesq system (14), then there exist fields in the upper part (u u , v u , w u , η) (resp. (u b , v b , w b , η)) satisfying the Navier-Stokes equations (3) in the upper part (resp. the boundary layer) of the domain and the interface continuity at least as O(ε).

Proof. We start from the final Boussinesq model ( 14) written in 1D:

u u t + η x - η xxx Bo + εu u u u x -εη xtt (z 2 -1) 2 = O(ε 2 ) + O(ε/Bo), η t + u u x (x, z, t) -ε 2 η xxt (z 2 -1 3 ) + ε(u u η) x -ε √ πR u u x * 1 √ t + 2ε √ π +∞ γ ′′ =0 u b,0 x (x, γ ′′ ) -u u,0 x (x, z = 0) √ R 4t γ ′′ γ ′ =0 e -γ ′2 dγ ′ dγ ′′ = o(ε), (25)
where u u is the horizontal velocity in the upper part (z ∈ (ε, 1 + εη(x, t))), η(x, t) is the free boundary's height, u u,0 (x, z) (resp. u b,0 (x, γ)) is the initial horizontal velocity in the upper part (resp. in the boundary layer). In the boundary layer, we set z = εγ.

First, one must justify that η tt = η xx + O(ε), which is easy from the zeroth order of (25) 1 . Then, the equation (25) 1 enables to prove the following Lemma. Lemma 9. A localized solution of (25) 1 is such that

1 0 u = u(x, z, t) -εη xt (z 2 -1/3)/2 + O(ε 2 ), u(x, 0, t) = u(x, z, t) -εη xt z 2 /2 + O(ε 2 ), u(x, 1, t) = u(x, z, t) + εη xt (1 -z 2 )/2 + O(ε 2 ).
The way to prove this Lemma is identical to the proof of Lemma 11 of [START_REF] Le Meur | Derivation of a viscous Boussinesq system for surface water waves[END_REF]. One differentiates (25) 1 with respect to z, integrate with respect to time t, prove that the constants of integration vanish if the wave is localized, and use that η tt = η xx + O(ε) to get u z = εη xt z + O(ε 2 ). Completing the proof is then easy.

Thanks to the previous Lemma, if we define p u (x, z, t) = (η -η xx /Bo) -εη tt (z -1) + ε The function u b is defined as u b -u u (z = 0) being the solution to the heat equation. Because it is only the zeroth order of the conservation of momentum in the boundary layer (3) 1 , and also because the limit between the boundary layer and the upper part is not made sufficiently precise through the lift "function" u u (x, z = 0, t), one may not prove more.

If we define: p b (x, γ, t) = η(x, t) -η xx /Bo + εη xx /2 w b (x, γ, t) = -ε 2 γ 0 u b x (x, γ ′ , t) dγ ′ , then (3) 1 in the boundary layer reduces to: The overlapping of the boundary layer and the upper part makes that consistency cannot be justified more precisely than O(ε). Yet, since the interface condition (for instance at z = √ ε) between these two domains is satisfied, the proof is complete, up to order O(ε).

(u b -u u ) t -(u b -u u ) γγ /R -εu u u u x +

Conclusion

We modeled the water waves flow by the Navier-Stokes equations in a 3D geometry with a flat bottom and free surface. We assumed the order of magnitude of the velocity in the transverse direction and deduced the scaling in this coordinate and discussed its influence on the results (Remark 3). The linear theory gave a new dispersion relation and a new phase velocity (Proposition 1). Surface tension was discussed and seemed often more relevant than viscosity, although these two parameters act very differently. We stated the associated Boussinesq system. This intermediate system enabled us to derive the KdV equation (Proposition 5) and the KP equation or system (Proposition 6). We justified why the zero-mass constraint does not raise any trouble for viscous KP, using [START_REF] Molinet | Remarks on the mass constraint for KP-type equations[END_REF]. Using the previous computations, we derived the Boussinesq equation (Proposition 7). Reciprocally, we proved that, should they exist, the solutions to the non-isotropic Boussinesq system are consistent with the Navier-Stokes equations (Proposition 8). These systems and equations are derived for a viscous fluid and take the surface tension into account.

Figure 1 :

 1 Figure 1: The dimensionless domain

h

  crit-Re l crit-Bo water 20 o C 4.6 10 -5 2.7 10 -3 mercury 20 o C 1.1 10 -5 1.9 10 -3

2 +∞ 0 f 0 f 0

 2000 xt (x, z ′′ , t)dz ′′ dz ′ w u (x, z, t) = ε(η t + εη x u u (z = 1 + εη)) -ε z 1+εη u ux , then (u u , w u , p u ) satisfies (3) 1,3,5,6 up to O(ε 2 ) + O(ε/Bo), and (3) 4 exactly in the upper part.We must now check (3) in the boundary layer, where we define z = εγ andu b (x, γ, t) = u u (x, z = 0, t) + √ R (x, γ ′ ) e -R(γ ′ -γ) 2 4t √ πt dγ ′ -u u (x, 0, .) * L -1 p→t (e - (x, γ ′ ) e -R(γ ′ +γ) 2 4t √ πt dγ ′ ,where f 0 (x, γ) = u b,0 (x, γ) -u u,0 (x, z = 0). It is proved in Lemma 6 of[START_REF] Le Meur | Derivation of a viscous Boussinesq system for surface water waves[END_REF] that this function satisfies(u b -u u (z = 0)) t -(u b -u u (z = 0)) γγ /R = 0, ⇒ u b t + η x -u b γγ /R -η xxx /Bo = O(ε) + O(ε/Bo). This u b also satisfies the initial conditions and the boundary conditions:    u b (x, γ, t = 0) = u b,0 (x, γ), u b (x, γ = 0, t) = 0, u b (x, γ → +∞, t) = u u (x, z = 0, t) (continuity condition).

εu b u b x -εu b γ γ 0 u b x + εη xxx / 2 +

 2 O(ε 2 ) + O(ε/Bo), = -εu u u u x + εu b u b x -εu b γ γ 0 u b x + εη xxx /2 + O(ε 2 ) + O(ε/Bo) = O(ε), where u u is evaluated at z = 0. So (u b , w b , p b ) satisfies (3) 1 only up to O(ε). So as to go further, one should use correctors, solve a nonlinear heat equation, make a more precise study at the interface between the two subdomains and use functionnal spaces. Simple computations prove that (u b , w b , p b ) satisfies (3) 3,4,7 up to O(ε 2 ) + O(ε/Bo).

Table 1 :

 1 Physical parameters for two fluidsGiven σ, ρ, ν, the Re and Bo number are evaluated in Table1. Depending on h or l, the Re or Bo dimensionless parameters are greater or less than 1. The critical values for each of them are given in Table2. It appears that viscosity is not so relevant as surface tension.

	Re	1/Bo

Table 2 :

 2 Critical values of length (in m) for which 1/Re or 1/Bo is 1.


Huicheng Zheng

Christophe Laurent

Grégoire Lefebvre

On the Basis Updating Rule of Adaptive-Subspace Self-Organizing Map (ASSOM)

This paper gives other views on the basis updating rule of the ASSOM proposed by Kohonen. We first show that the traditional basis vector rotation rule can be expressed as a correction to the basis vector which is proportional to component vectors in the episode. With the latter form, some intermediate computations can be reused, leading to a computational load only linear to the input dimension and the subspace dimension, whereas naive implementation of the traditional rotation rule has a computational load quadratic to the input dimension. We then proceed to propose a batch-mode updating of the basis vectors. We show that the correction made to each basis vector is a linear combination of component vectors in the input episode. Computations can be further saved. Experiments show that the proposed methods preserve the ability to generate topologically ordered invariant-feature filters and that the learning procedure is largely boosted.

Introduction

Adaptive-subspace self-organizing map (ASSOM) [START_REF] Kohonen | The adaptive-subspace som (assom) and its use for the implementation of invariant feature detection[END_REF] is basically a combination of the competitive selection and cooperative learning as in the traditional SOM [START_REF] Kohonen | Self-Organizing Maps. 3rd edn[END_REF] and a subspace method. The single weight vectors at map units in SOM are replaced by modules of basis vectors in ASSOM that span some linear subspaces. ASSOM is an alternative to the standard principal component analysis (PCA) method of feature extraction. An earlier neural approach for PCA can be found in [START_REF] Oja | Principal components, minor components, and linear neural networks[END_REF]. The ASSOM can generate spatially ordered feature filters thanks to spatial interactions among processing units [START_REF] Kohonen | Self-organized formation of various invariantfeature filters in the adaptive-subspace som[END_REF]. Each module in ASSOM can be realized as a neural network which receives input vectors and outputs their orthogonal projections on the represented subspaces.

The input to an ASSOM array is typically an episode, i.e. a sequence of pattern vectors supposed to approximately span some linear subspace. These vectors shall also be referred to as component vectors of the episode in this paper. By learning the episode as a whole, ASSOM is able to capture the transformation coded in the episode. The simulation results in [START_REF] Kohonen | The adaptive-subspace som (assom) and its use for the implementation of invariant feature detection[END_REF] and [START_REF] Kohonen | Self-organized formation of various invariantfeature filters in the adaptive-subspace som[END_REF] have demonstrated that ASSOM can induce ordered filter banks to account for translation, rotation and scaling. The relationship between the neurons in the ASSOM architecture and their biological counterparts are reported [START_REF] Kohonen | Self-organized formation of various invariantfeature filters in the adaptive-subspace som[END_REF]. ASSOM has been applied to speech processing [START_REF] Hase | Speech signal processing using adaptive subspace som (assom)[END_REF], texture segmentation [START_REF] Ruiz Del Solar | Texsom: texture segmentation using self-organizing maps[END_REF], image This work was carried out during the tenure of a MUSCLE Internal fellowship (http://www.muscle-noe.org). retrieval [START_REF] De Ridder | The adaptive subspace map for image description and image database retrieval[END_REF] and image classification [START_REF] De Ridder | The adaptive subspace map for image description and image database retrieval[END_REF], [START_REF] Zhang | Handwritten digit recognition by adaptive-subspace self-organizing map (assom)[END_REF], etc. in the literature. A supervised variant of ASSOM, called supervised adaptive-subspace self-organizing map (SASSOM), was used by Ruiz del Solar in [START_REF] Ruiz Del Solar | Texsom: texture segmentation using self-organizing maps[END_REF].

The basis vector rotation rule in the traditional ASSOM implementation takes a form of matrix multiplication. This rule is hard to understand and more seriously, naive implemenation of the basis vector rotation rule leads to a computational load which is quadratic to the input dimension, not to mention large amount of memory required by the usually high-dimensional matrix operations. This deficiency renders naive implementation of the basic ASSOM learning very costly for practical applications.

There were efforts in the literature to reduce the computational load associated with the basic ASSOM learning. De Ridder et al. [START_REF] De Ridder | The adaptive subspace map for image description and image database retrieval[END_REF] dropped topological ordering to reduce the computations involved in the cooperative learning. Furthermore, they performed a batch-mode updating of subspaces with PCA to avoid the time-consuming iterative updating. Similarly, López-Rubio et al. [START_REF] López-Rubio | A principal components analysis selforganizing map[END_REF] used PCA with ASSOM to perform PCA while retaining self-organization of generated features. The resulting algorithm is named PCASOM. According to their report, under similar classification performance, their algorithm runs about twice faster than the basic ASSOM. McGlinchey et al. [START_REF] Mcglinchey | Fast formation of invariant feature maps[END_REF] replaced the traditional basis vector updating formula with one proposed by Oja [START_REF] Oja | Neural networks, principal components and subspaces[END_REF]. According to their paper, the computational load is only linear to the input dimension, but quadratic to the subspace dimension.

In this paper, we first show that with the traditional basis rotation rule, the correction made to each basis vector is in fact a vector proportional to the component vector of the input episode. With this modified form, some intermediate computations can be reused, leading to a computational load only linear to both the input dimension and the subspace dimension. We then proceed to propose a batch-mode updating of the basis vectors, where the correction made to each basis vector is a linear combination of component vectors in the episode. This modified rule further accelerates the learning procedure by saving large amounts of computations. This paper will be organized as follows: In Sect. 2, we review briefly the basic ASSOM learning procedure. The proposed alternative updating rules will be presented in Sect. 3. Section 4 is dedicated to experiments which demonstrate the performance of the proposed methods. This paper will be concluded by Sect. 5.

The Basic ASSOM Learning

An ASSOM is composed of an array of modules. Each module in the ASSOM can be realized by a two-layered neural network [START_REF] Kohonen | Self-organized formation of various invariantfeature filters in the adaptive-subspace som[END_REF], as shown in Fig. 1. It calculates the projection of an input vector x on the subspace L of the module. Supposing L is spanned by a set of basis vectors {b 1 , b 2 , . . . , b H }, where H is the dimension of L, the neurons in the first layer take the orthogonal projections x T b h of the input vector x on the individual basis vectors b h . The basis vectors are supposed to be orthonormalized. The only quadratic neuron of the second layer sums up the squared outputs of the first-layer neurons. The output of the module is then xL 2 , the squared norm of the projection of x on the subspace L. It can be regarded as a measure of the matching between the input vector x and the subspace L. For an input episode x(s), s ∈ S, where S is the index set of vectors in the episode, Kohonen proposed to use the energy s∈S xL (s) 2 , i.e. the sum of squared norms of the projections of the individual vectors, as the measure [START_REF] Kohonen | Self-organized formation of various invariantfeature filters in the adaptive-subspace som[END_REF].

1 b 2 b H b Q x 1 T b x 2 ˆL x 2 T b x H b x T Fig. 1.
A module of ASSOM realized as a neural network. x is an input vector. {b1, b2, . . . , bH } is an orthonormal basis of the linear subspace L of the module. Q is a quadratic neuron that sums up squares of its inputs

The learning process of ASSOM approximately minimizes an error function in an iterative way [START_REF] Kohonen | Self-organized formation of various invariantfeature filters in the adaptive-subspace som[END_REF]. The iteration at step t of the basic ASSOM learning procedure proceeds as follows:

1. For the episode x(s), s ∈ S, locate the winning module indexed by c = arg max i∈I s∈S xLi (s) 2 , where I is the index set of modules in the AS-SOM. 2. For each module i in the neighborhood of c, including c itself, update the subspace L i for each component vector x(s), s ∈ S, that is, update the basis vectors b

(i)
h , according to the following rules: (a) Rotate each basis vector according to:

b (i) h = P (i) c (x, t)b (i) h , (1)
where b

(i)
h is the new basis vector and b

(i)
h the old one. The matrix P (i) c (x, t) is a rotation operator defined by:

P (i) c (x, t) = I + λ(t)h (i) c (t) x(s)x T (s) xLi (s) x(s) , ( 2
)
where I is the identity matrix, λ(t) a learning-rate factor that diminishes with t. h h to improve stability of the results [START_REF] Kohonen | Self-organized formation of various invariantfeature filters in the adaptive-subspace som[END_REF] and then orthonormalize these basis vectors.

Through this competitive and cooperative learning procedure, the ASSOM will finally arrive at a topologically organized status, where nearby modules represent similar feature subspaces. Naive implementation of (1) requires a matrix multiplication which needs not only a large amount of memory, but also a computational load quadratic to the input dimension. It would be costly for practical applications of ASSOM.

On the Basis Updating Rule of ASSOM

Insight on the Basis Vector Rotation

In the first place we propose to replace the formulae (1) and ( 2) through a little mathematical deduction. The term b (i) h in (1) can be distributed to the right side of (2), leading to the current basis vector and a correction to it:

b (i) h = b (i) h + ∆b (i) h , (3)
where

∆b (i) h = λ(t)h (i) c (t) x(s)x T (s)b (i) h xLi (s) x(s) . (4)
x T (s)b

(i)
h is in fact a scalar value. The equation can be rewritten as:

∆b (i) h = α (i) c,h (s, t)x(s) . (5)
Here α

(i)
c,h (s, t) is a scalar value defined by:

α (i) c,h (s, t) = λ(t)h (i) c (t) x T (s)b (i) h xLi (s) x(s) . ( 6
)
This shows that the correction ∆b

(i)
h is in fact proportional to the component vector x(s), as illustrated in Fig. 2, which seems to have been ignored by many practitioners. Equation ( 5) gives a clearer way to understand the basis vector rotation in ASSOM learning than the traditional rotation matrix [START_REF] Kohonen | Self-Organizing Maps. 3rd edn[END_REF]. Note that in (6), x T (s)b

(i)
h is the projection of the component vector on the basis vectors represented by the neurons of the first layer, which we have already when computing the projection xLi (s) (cf. Fig. 1). If we calculate the scaling factor α (i) c,h (s, t) first, and then scale the component vector x(s) with this factor, the computations associated with the basis vector updating will be dramatically reduced. This implementation will be referred to as FL-ASSOM for fast-learning ASSOM. It is completely equivalent to the basic ASSOM in terms of generating topologically ordered invariant-feature filters.

Let us compare the computational loads of the basis vector updating in the basic ASSOM and FL-ASSOM by analyzing the respective updating formulae. We assume the input dimension to be N , and the subspace dimension to be M . We first evaluate the computations required by naive implementation of the traditional basis vector updating rule (1). For each component vector x(s), x(s)x T (s) in (2) needs N 2 multiplications, the matrix multiplication in (1) amounts to M N 2 multiplications. There are totally about M N 2 + N 2 multiplications. Similarly, the number of additions required by [START_REF] Kohonen | The adaptive-subspace som (assom) and its use for the implementation of invariant feature detection[END_REF] can be shown to be around M N 2 + N 2 . Finally, the computational load of naive implementation of the traditional updating rule is approximately O(M N 2 ), i.e. quadratic to the input dimension and linear to the subspace dimension. The replacement proposed by McGlinchey et al. [START_REF] Mcglinchey | Fast formation of invariant feature maps[END_REF] leads to a computational load of O(M 2 N ), as shown in their paper, i.e. linear to the input dimension but quadratic to the subspace dimension. Now with the proposed updating rule (5), for each component vector x(s), the computations of xLi (s) and x(s) in ( 6) need about M N +2N multiplications, and α (i) c,h (s, t)x(s) in ( 5) about M N multiplications. In all (5) needs about 2M N + 2N multiplications. Similarly, the number of additions can be shown to be about 2M N + 2N . So with (5), the computational load is approximately O(M N ), i.e. linear to both the input dimension and the subspace dimension. So there is an obvious benefit in using (5) other than naive implementation of (1).

) (i h b ) ( ' i h b ) ( ) , ( ) ( , ) ( s t s i h c i h x b α = ∆ ) (s x

Further Boosting: Batch-mode Basis Vector Updating

Basis vector updating can be further boosted by working in a batch mode. We can avoid computing the value of xLi (s) in ( 6) by using the value computed previously during module competition. However this could not be done inside the framework of FL-ASSOM since the subspaces are continuously changing in receiving each component vector of the episode. To save computation of xLi (s) , the following batch-mode rotation operator [START_REF] Kohonen | Self-organized formation of various invariantfeature filters in the adaptive-subspace som[END_REF] will be useful:

B (i) c (t) = I + λ(t)h (i) c (t) s∈S x(s)x T (s) x(s) 2 . (7)
With this rotation operator, each basis vector in the subspace will be rotated only once for the whole input episode. For stability of the solution, we expect the magnitude of the correction made to the basis vectors to be monotonically decreasing with respect to xLi (s) . We borrow the idea from the basic rotation operator P (i) c (x, t) [START_REF] Kohonen | Self-organized formation of various invariantfeature filters in the adaptive-subspace som[END_REF] to divide the learning-rate factor λ(t) by the scalar value xLi (s) / x(s) , which only changes the effective learning rate. The batch-mode rotation operator then becomes:

B (i) c (t) = I + λ(t)h (i) c (t) s∈S x(s)x T (s) xLi (s) x(s) . (8)
As for FL-ASSOM, we distribute b

(i) h
to terms in this operator. With similar deduction as in FL-ASSOM, the basis vector updating rule becomes:

b (i) h = b (i) h + ∆b (i) h , (9)
where ∆b

(i) h = s∈S α (i) c,h (s, t)x(s) . (10)
The correction made to each basis vector is thus a linear combination of the component vectors in the episode. The difference between the updating rule (9) here and ( 3) is that the former updates the basis vectors in a batch mode for the whole episode while the latter updates the basis vectors for each component vector one by one. The scalar parameter α

(i)
c,h (s, t) has the same form as [START_REF] Ruiz Del Solar | Texsom: texture segmentation using self-organizing maps[END_REF] in FL-ASSOM:

α (i) c,h (s, t) = λ(t)h (i) c (t) x T (s)b (i) h xLi (s) x(s) . ( 11
)
The meaning of this equation is a little different from that of ( 6), where the subspace L i (s) of the module i should be updated for each component vector x(s) in the episode and thus we could not reuse the computational results of module competition. Here in [START_REF] Oja | Neural networks, principal components and subspaces[END_REF] the basis vector updating works in a batch mode, i.e. updating is performed only after the whole episode has been received. Therefore, xLi (s) and x T (s)b (i) h can reuse the results previously calculated during module competition. What we need to do is only store the calculated values in registers and fetch them when needed. The computational load of ( 11) is thus trivial. Furthermore, the dissipation as well as orthonormalization of basis vectors can be performed only once for each episode without loosing accuracy since the basis vectors are not updated during the episode. The computational load can thus be further reduced. This method will be referred to as BFL-ASSOM for batch-mode fast-learning ASSOM.

Let us estimate the computational load of BFL-ASSOM. For basis vector updating with [START_REF] Mcglinchey | Fast formation of invariant feature maps[END_REF], we estimate the computational load averaged on each component vector of the episode as we did for the basic ASSOM and FL-ASSOM. As has been mentioned, the calculation of α (i) c,h (s, t) according to [START_REF] Oja | Neural networks, principal components and subspaces[END_REF] needs only trivial computation. The majority of computation is in [START_REF] Mcglinchey | Fast formation of invariant feature maps[END_REF]. Averaged on each vector in the episode, the computational load required by basis vector updating with BFL-ASSOM is about M N multiplications and M N additions. Furthermore, since the dissipation and orthonormalization of basis vectors can be performed only once for each episode, the whole learning time can be further reduced.

Experiments

We first show that BFL-ASSOM can also generate the topologically ordered invariantfeature filters as the basic ASSOM. The results of FL-ASSOM will be shown as the ground truth since FL-ASSOM is mathematically equivalent to the basic ASSOM. One of the most common transformations occurred to images is translation. We will show that BFL-ASSOM permits to generate Gabor type filters from episodes subject to translation.

The input episodes are constructed from a colored noise image, which is generated by filtering a white noise image with a second-order Butterworth filter. The cut-off frequency is set to 0.6 times of the Nyquist frequency of the sampling lattice. Each episode is composed of 6 vectors, each of which is formed on a circular receptive field on the sampling lattice composed of 349 pixels. The vectors in the same episode have only random translation of no more than 5 pixels in both the horizontal and the vertical directions. The episodes are generated on random locations of the colored noise image. The mean value of components of each input vector is subtracted from each component of the vector. In order to symmetrize the filters with respect to the center of the receptive field, the input samples are weighted by a Gaussian function symmetrically placed at the center of the receptive field with a full width at half maximum (FWHM) that varies linearly with respect to the learning step t from 1 to 16 sampling lattice spacings. Each vector is normalized before entering into the ASSOM array. The ASSOM array is composed of 9 × 10 modules aligned in a hexagonal lattice with two basis vectors at each module. The basis vectors of all the modules are initialized randomly and orthonormalized at the beginning of the learning process. The radius of the circular neighborhood function h (i) c (t) decreases linearly from 6.73 (= 0.5 × (9 2 + 10 2 ) 1/2 ) to 0.9 ASSOM array spacings with t. The learning-rate factor has the form λ(t) = 0.1 • T /(T + 99t), where T is the total number of learning steps and set to 30, 000 for the current experiment.

The translation-invariant filters generated by BFL-ASSOM compared to those by FL-ASSOM are shown in Fig. 3(a) with a gray scale. We can see that both methods generated topologically ordered Gabor-like filters. For either method, the two basis vectors at the same array locations have the same frequencies but 90 degrees of phase difference. Figure 3(b) shows how the average projection error e changes with the learning step t for FL-ASSOM and BFL-ASSOM. For each input episode X = {x(s), s ∈ S}, the projection error is calculated according to e(X) = s∈S x(s)-x(s) 2

x(s) 2

, where x(s) is the orthogonal projection of x(s) on the subspace of the winning module. e is the average of e(X) over all the training episodes. We can see that the curves of FL-ASSOM and BFL-ASSOM match very well in Fig. 3(b), which reveal that their difference in terms of generating the filters is indeed very little.

In the second experiment, we compare the computational loads of the basic AS-SOM, FL-ASSOM and BFL-ASSOM. We designed the experiment by using C++ implementations of all the methods. In this experiment, the input dimension as well as the subspace dimension vary. We count the elapsed CPU seconds for different methods. The number of iterations are fixed to 1, 000. Each episode is composed of 6 vectors. These vectors are generated randomly according to a uniform probability distribution. The rectangular ASSOM array contains 10 × 10 modules.

The timing results are summarized in Table 1. As was anticipated, the time of updating basis vectors with the basic ASSOM increased sharply with the input dimension and moderately with the subspace dimension. Basis vector updating is the bottleneck of the basic learning procedure, especially when the input dimension is high. With FL-ASSOM, it is clear that the time of updating basis vectors increases much more mod- erately with the input dimension. The response to the subspace dimension is also quite mild. Basis vector updating is no longer a bottleneck for the learning procedure. As a result, the learning time drops dramatically compared to the basic ASSOM. However learning time outside basis vector updating is not reduced. Now with BFL-ASSOM, we can observe that the basis vector updating time is further reduced. Moreover, learning time outside the basis vector updating is also reduced considerably compared to the basic ASSOM and FL-ASSOM.

The relationship between the basis vector updating time and the input dimension or the subspace dimension for the three implementations of ASSOM is visualized in Fig. 4. The basis vector updating time increases approximately linearly with respect to the input dimension for FL-ASSOM and BFL-ASSOM, but apparently nonlinearly for the basic ASSOM. In all the cases, the updating time increases approximately linearly with respect to the subspace dimension.

Conclusions

The focus of this paper is on the basis updating rule of the ASSOM learning. We first showed that the traditional basis rotation rule amounts to a correction made to the basis vectors which is proportional to the component vectors of the input episode. This gives us a better understanding of the basis updating in ASSOM learning. With this modified form of updating rule, some computations can be saved by reusing some intermediate computations. The resulting method is referred to as FL-ASSOM. Naive implementation of the traditional basis updating rule leads to a computational load linear to the subspace dimension but quadratic to the input dimension. This computational load is reduced by FL-ASSOM to be linear to both the subspace dimension and the input dimension. The ability of FL-ASSOM in generating topologically ordered invariant-feature filters is altogether preserved since FL-ASSOM is mathematically equivalent to the basic ASSOM. We then proceeded to present BFL-ASSOM, where the basis vectors are updated in a batch mode. We showed that the correction made to each basis vector is a linear combination of the component vectors in the input episode. What's more, large amount of computations can be further saved by reusing more of previous computations and performing only one dissipation and orthonormalization for each episode.

Our experiments showed that BFL-ASSOM can also generate topologically ordered Gabor-like translation-invariant filters and that the bottleneck of the basis vector updating in the learning procedure is totally removed by FL-ASSOM and BFL-ASSOM. The proposed methods can be easily adapted to the supervised ASSOM used in [START_REF] Ruiz Del Solar | Texsom: texture segmentation using self-organizing maps[END_REF]. There is an obvious benefit in using the proposed rules instead of naive implementation of the basic learning rule.

c

  (t) is a neighborhood function defined on the ASSOM lattice. (b) Dissipate the basis vectors b

Fig. 2 .

 2 Fig. 2. An alternative view of the basis vector updating rule of ASSOM. The correction ∆b (i) h is proportional to the component vector x(s). After updating, b (i) h represents better the component vector x(s)

Fig. 3 .

 3 Fig. 3. (a) The Gabor type filters generated by BFL-ASSOM compared to those by FL-ASSOM on episodes subject to translation. Top: Filters generated by FL-ASSOM; Bottom: Filters generby BFL-ASSOM. Left: First basis vectors. Right: Second basis vectors. (b) Change of the projection error e according to the learning step t for FL-ASSOM and BFL-ASSOM

Table 1 .

 1 The timing results for the basic ASSOM, FL-ASSOM and BFL-ASSOM. VU (vector updating time) denotes the time for the basis vector updating. WL (whole learning time) denotes the time for the whole learning procedure, including those for module competition, basis vector dissipation and orthonormalization. All the times are given in seconds

			The basic ASSOM
		M =2		M =3		M =4
		VU	WL	VU	WL	VU	WL
	N =50	29.36 41.27 30.44 47.53 32.45 54.72
	N =100 134.92 154.95 145.95 172.85 149.91 188.06
	N =200 742.34 786.63 769.56 828.47 814.80 895.08
	N =400 4529.69 4626.43 4956.64 5090.56 5200.78 5367.35
			FL-ASSOM
		M =2		M =3		M =4
		VU	WL	VU	WL	VU	WL
	N =50	1.53 13.31	2.31 18.92	3.18 25.03
	N =100	2.39 21.09	3.15 30.41	3.86 40.80
	N =200	3.30 37.86	55.44	5.93 73.98
	N =400	5.68 70.88	7.51 105.01	9.92 139.25
			BFL-ASSOM
		M =2		M =3		M =4
		VU	WL	VU	WL	VU	WL
	N =50	0.58	4.83	1.28	6.72	1.33	8.75
	N =100	0.87	6.88	1.50 10.01	1.58 13.37
	N =200	1.03 11.61	1.67 17.37	2.10 22.86
	N =400	1.46 21.02	2.06 31.01	2.99 41.51


Grégoire Lefebvre

Samuel Berlemont

Franck Mamalet

Christophe Garcia

BLSTM-RNN based 3D Gesture Classification

Keywords: LSTM-RNN, DTW, HMM, MEM, hand gesture recognition

This paper presents a new robust method for inertial MEM (MicroElectroMechanical systems) 3D gesture recognition. The linear acceleration and the angular velocity, respectively provided by the accelerometer and the gyrometer, are sampled in time resulting in 6D values at each time step which are used as inputs for the gesture recognition system. We propose to build a system based on Bidirectional Long Short-Term Memory Recurrent Neural Networks (BLSTM-RNN) for gesture classification from raw MEM data. We also compare this system to a geometric approach using DTW (Dynamic Time Warping) and a statistical method based on HMM (Hidden Markov Model) from filtered and denoised MEM data. Experimental results on 22 individuals producing 14 gestures in the air show that the proposed approach outperforms classical classification methods with a classification mean rate of 95.57% and a standard deviation of 0.50 for 616 test gestures. Furthermore, these experiments underline that combining accelerometer and gyrometer information gives better results that using a single inertial description.

Introduction

Accelerometers and gyrometers are nowadays present in our everyday Smartphones. These sensors capture hand movements when users grasp their devices. We can consider two main issues: posture recognition and symbolic gesture recognition. In the first case, the user maintains a posture during a certain period of time, describing for instance the fact that the device is upside down. In the second situation, the user may produce a gesture to execute a system command, like drawing a heart symbol in 3D space to call its favorite phone number. Dynamic gesture recognition based on inertial sensors is a very challenging task. Algorithms are confronted to numerous factors causing errors in the recognition process: dynamical differences (intensive versus phlegmatic gestures), temporal differences (slow versus fast movements), physical constraints (device weight, human body elasticity, left or right-handed, seated or standing up, on the move, etc.), classification constraints (mono versus multi users, open or closed world paradigm, etc.). Classically, several steps operate from signal data preprocessing to gesture classification with some intermediate steps like data clustering and gesture model learning. The preprocessing steps aim at reducing the input signals that characterize the corresponding gestures. Different methods can then be applied: calibration, filtering, normalization or vectorization. Data clustering is often applied to reduce the input space dimension and find class referent gesture vectors. A learning phase of a gesture model follows this clustering step and finally a decision rule or a specific classifier is built to label the input data as a recognized gesture or an unknown gesture. In this article, we propose to learn an efficient gesture classifier without any preprocessing method (i.e. from raw MEM data) using a BLSTM-RNN model.

This paper is organized as follows. In Section 2, sensor-based gesture recognition is described with a survey. Section 3 presents our recognition method. Section 4 describes the experimental results. Finally, conclusions are drawn.

2 Accelerometer based 3D Gesture Recognition 3D gesture recognition using accelerometers has been studied in recent years, and for gesture classification three main strategies stand out which are based on statistics, on geometry or on boosting classifier approaches.

The first strategy has been deeply studied in the last decade with two main approaches: discrete versus continuous HMM [START_REF] Hofmann | Velocity profile based recognition of dynamic gestures with discrete hidden markov models[END_REF][START_REF] Kallio | Online gesture recognition system for mobile interaction[END_REF][START_REF] Kela | Accelerometer-based gesture control for a design environment[END_REF][START_REF] Pylvänäinen | Accelerometer Based Gesture Recognition Using Continuous HMMs Pattern Recognition and Image Analysis[END_REF]. Hofmann et al. [START_REF] Hofmann | Velocity profile based recognition of dynamic gestures with discrete hidden markov models[END_REF] proposed to use discrete HMM (dHMM) for recognizing dynamic gestures thanks to their velocity profile. This approach consists of two levels and stages of recognition: a low-level stage essentially dividing the input data space into different regions and assigning each of them (i.e. creation of a vector codebook), and a high-level stage taking the sequences of vector indexes from the first stage and classifying them with discrete HMM. The experiments are built using a training set with 10 samples per gesture, each sample representing hand orientation, acceleration data and finger joint angle. A vector codebook is obtained by an input space clustering method (i.e. K-means algorithm). Clustering essentially serves as an unsupervised learning procedure to model the shape of the feature vector distribution in the input data space. Here, the number of HMM states vary from 1 to 10 and the observation alphabet size equals to 120. The comparison between ergodic HMM and left-to-right HMM shows similar results with 95.6% correct recognition rate for 100 gestures. Similar results are presented in [START_REF] Kallio | Online gesture recognition system for mobile interaction[END_REF][START_REF] Kela | Accelerometer-based gesture control for a design environment[END_REF]. Kallio et al. [START_REF] Kallio | Online gesture recognition system for mobile interaction[END_REF] use 5 HMM states and a codebook size of 8 for 16 gestures. The authors highlight that the performances decrease when using 4 sequences for training the system compared to 20 sequences. The recognition rate falls from 95% to 75% even for this mono-user case study. In [START_REF] Kela | Accelerometer-based gesture control for a design environment[END_REF], a 37 multi-user case is studied with 8 gestures, evaluating the effect of vector quantization and sampling. A rate of 96.1% of correct classification is obtained with 5 HMM states and a codebook size of 8. However, this study can be seen as biased since the K-means clustering is performed from all the available data set and not only the training database. In opposition to the previous studies, and to take into consideration that gesture data are correlated in time, Pylvänäinen proposes in [START_REF] Pylvänäinen | Accelerometer Based Gesture Recognition Using Continuous HMMs Pattern Recognition and Image Analysis[END_REF] to build a system based on continuous HMM (cHMM). Again, the results are convincing, with 96.76% on a dataset providing 20 samples for 10 gestures realized by 7 persons.

The second strategy for recognizing 3D gestures is based on geometric models with distance computation. The goal is to provide a gallery of some gesture references to model each gesture class and design a decision rule for a test gesture regarding the respective distance to these referent instances. On the contrary to the HMM strategy, no learning phase is needed but computational time is required for a test gesture to be compared to all referent instances. Consequently, the main drawback of this approach is the necessity to find the most relevant samples to represent a gesture class while keeping the number of these referents low in order to minimize the final evaluation processing time. Wilson et al. in [START_REF] Wilson | Gesture recognition using the xwand[END_REF] compare Linear Time Warping (LTW) and Dynamic Time Warping (DTW) to the HMM based strategy. Their experiment with 7 types of gesture from 6 users shows an advantage for HMM with 90% in opposition to the score of LTW and DTW of respectively 40% and 71%. Liu et al. experiment with more success the DTW strategy in [START_REF] Liu | uwave: Accelerometer-based personalized gesture recognition and its applications[END_REF]. Gesture recognition and user identification are performed with good recognition rates of respectively 93.5% and 88%. The authors introduce an averaging window of 50 ms for reducing noise and erratic moves. The gesture data, performed over multiple days, consists of 30 samples of 8 gestures for 8 individuals and the user recognition results are obtained from 25 participants. Likewise, in [START_REF] Akl | Accelerometer-based gesture recognition via dynamic-time warping, affinity propagation, & compressive sensing[END_REF], Akl et al. use DTW and affinity propagation for dimension reduction for recognizing 3D gestures. 7 subjects participated producing 3700 gesture traces for a good classification rate of 90%.

The third strategy for recognizing 3D gestures is to learn a specific classifier. Hoffman et al. (see [START_REF] Hoffman | Breaking the status quo: Improving 3d gesture recognition with spatially convenient input devices[END_REF]) improve 3D gesture recognition with a linear classifier and Adaboost, inspired by the method proposed in [START_REF]A practical approach for writer-dependent symbol recognition using a writerindependent symbol recognizer[END_REF] for 2D symbol writer recognition. The experiments show an accuracy of 98% for 13 gestures made by 17 participants. Other studies focus on SVM (i.e. Support Vector Machine) like in [START_REF] Wu | Gesture recognition with a 3-d accelerometer[END_REF]. This study uses frame-based descriptors. Each gesture is divided into segments where are computed to form descriptors: mean, energy, entropy, standard deviation and correlation. These descriptors constitute the feature vector to be classified by a multi-class SVM. The obtained results are 95.21% of good recognition for 12 gestures made by 10 individuals.

Consequently, many strategies are explored with different paradigms and specific data processing methods on different databases. Nevertheless, theses approaches suffer from finding automatically the relevant parameters (e.g. signal processing, etc.) to deal with gesture variabilities. We develop hereafter our 3D gesture recognition method based on BLSTM-RNN from raw input data and compare it with classical methods on a common database.

3 The proposed 3D Gesture Recognition Method

Bidirectional Long Short-Term Memory RNNs

Classical RNNs are a common learning technique for temporal analysis of data since they are able to take into consideration the temporal context. This is achieved by using recurrent connections within the hidden layer which allow the network to remember a state representing the previous input values. However, Hochreiter and Schmidhuber in [START_REF] Schmidhuber | Long short-term memory[END_REF] have shown that if RNNs can handle shorttime lags between inputs, the problem of exponential error decay prevent them from tackling real-life long-term dependencies. They introduced thus the Long Short Term Memory RNNs, that allows a constant error signal propagation through time using a special node called constant error carousel (CEC) and multiplicative gates (Fig 1 .a). These gates are neurons that can set (input gate), reset (forget gate) or hide (output gate) the internal value of the CEC according to neuron input values and context.

LSTM-RNNs have proven their great ability to deal with temporal data in many applications (e.g. phoneme classification [START_REF] Graves | Framewise phoneme classification with bidirectional lstm and other neural network architectures[END_REF], action classification [START_REF] Baccouche | Sequential Deep Learning for Human Action Recognition[END_REF]). In this paper we consider gesture data using 6D input vectors through sampling timestep. These data are correlated during the user gestural production, and time lags between the beginning and the end of gesture can be long. For these reasons, LSTM-RNN is chosen to classify the input MEM data sequence. Furthermore, since gesture recognition, at a given timestep, may depend on past and future context, we use Bidirectional LSTM-RNN (BLSTM-RNN), introduced in [START_REF] Graves | Framewise phoneme classification with bidirectional lstm and other neural network architectures[END_REF], that consists in two separate hidden layers, the forward (resp. backward) layer able to deal with past (resp. future) context. The output layer is connected to both hidden layers in order to fuse past and future contexts.

BLSTM-RNN Architecture, Training and Decision Rule

The proposed gesture classification scheme based on BLSTM-RNN is described in Figure 1.b. First, the input layer consists in the concatenation of accelerometer and gyrometer information synchronized in time (i.e. 6 input values per timestep). Notice that our system relies only on the raw MEMs data, without any preprocessing in opposition to most of state-of-the-art methods. These data are linearly normalized between -1 and +1 according to the maximum value that sensors can provide. The forward and backward LSTM hidden layers are fully connected to the input layer and consist in 100 LSTM neurons each with full recurrent connections. The output layer has a size equals to the number of gesture to classify. The SoftMax activation function is used for this layer to give network responses between 0 and 1 at every timestep. Classically, these outputs can be considered as posterior probabilities of the input sequence to belong to a specific category at a given timestep. This network is learned using classical on-line backpropagation through time with momentum (i.e. learning rate 5e-4, momentum 0.2), as described in [START_REF] Schmidhuber | Long short-term memory[END_REF], on a training set, by targeting the same corresponding gesture class at each time step for each input example. For evaluation of a new gesture sequence, we use a majority voting rule over the outputs along the sequence (i.e. keeping only the most probable class at each time step) to determine the final gesture class.

Experimental Results

There is no public dataset for comparison of 3D gesture recognition. Therefore, we have collected our 3D gesture dataset to compare classification methods.

Our dataset has been captured on an Android Nexus S Samsung device. 22 participants, from 20 to 55 years old, all right-handed, performed 5 times each of the 14 symbolic gestures. This corresponds to 1540 temporal segmented gestures. The sampling time for accelerometer and gyroscope capture is 40 ms. The 14 symbolic gestures are divided into 2 families: linear gestures (e.g. north, south, east and west flicks, and up, down, pick and throw gestures) and curvilinear gestures (e.g. alpha, heart, letter N, letter Z, clockwise and counter-clockwise). These choices make the dataset difficult. There are classically confusions between flick gestures and letter N and Z. Likewise, the clockwise movement is often confused with alpha or heart symbols. Hereafter, we use temporal segmented gestures where only useful data are efficient to classify the inputs. We use 3 different configurations to compare our solution based on BLSTM-RNN to 3 state-of-the-art solutions: DTW, dHMM and cHMM based methods. The DTW solution uses a 5 nearest neighbor classification [START_REF] Petit | GRASP: Moteur de reconnaissance de gestes[END_REF] and the HMM solution uses the maximum of likelihood as a decision rule . In all experiments, we use a filtered and vectorized gestural information for these methods and raw MEM information for LSTM solution. In the following, we use a 3-fold cross validation.

The first configuration (DB1) corresponds to the personalization paradigm, where only one user is considered with few learning examples. For this configuration we have used the 70 gestures of a single participant in the learning phase, and ask him to process 16 more instances of each gesture for test (i.e. 224 gestures). The second configuration (DB2) uses 3 instances of each gesture per user for the learning phase: 924 gestures (i.e. 60% of all data) are used for the learning phase and 616 gestures (i.e. 40%) for the test phase. This case corresponds to a multi-user system and a closed world paradigm. The third configuration (DB3) is composed of all samples from 17 users (i.e. 1190 gestures) and the test data uses the other available gestures (i.e. 350 gestures from unknown users). This case is close to a real system trained with a few examples and having to generalize to new users who want to use it without any personalization phase. Here, the configuration represents the open world paradigm. Classification Results Table 1 outlines the global performances of each classifier for configurations DB1, DB2 and DB3 coupling or not accelerometer and gyrometer data. Considering coupled input data (accelerometer+gyroscope), this table shows that our BLSTM-RNN based classifier gives the best results on DB2 and DB3, with respectively 95.57 ± 0.50% and 92.57 ± 2.85%.

In the three configurations, the dHMM solution provides lower performances which is mainly due to the input data variability and the complexity to determine an automatic discriminant codebook.

On two configurations (DB2 and DB3), the DTW solution achieves the second best performance in mean recognition rate before the cHMM based one.

On DB1 configuration, DTW and cHMM achieve equivalent performances while our BLSTM-RNN approach is less efficient. This is mainly due to the lack of learning data which leads to the classical over-fitting issue. The attempts made with smaller LSTM networks did not allow any improvement on generalization.

When comparing these methods using a single input MEM sensor (accelerometer or gyroscope), we can see that using only gyroscope data is less efficient than using single accelerometer data. Moreover, when these two information are combined, the performances increase with respectively 99.70 ± 0.42%, 94.04 ± 0.15% and 91.71 ± 1.46%, for instance, for the DTW based method on DB1, DB2 and DB3 configurations.

Main conclusions of a deep analysis of confusion matrices (not provided here due to lack of space) are the following. The main drawback for the cHMM based method in this context is the incorrect classification of the N gestures with only 0.95% of correct classification. 62.86% of the N gestures are confused with the pick gestures. A strong confusion appears with opposite gestures as pick and throw or down and up gestures. Opposite gestures may be mis-classified when some user anticipate a north flick gesture by slightly moving back the device in the beginning of the production. On the contrary, the DTW based method provide a good solution to classify linear gestures except for the throw gesture which is often recognized as east and north flicks, which can be explained by the similar nature of production of these three gesture types. Our BLSTM-RNN approach have some issue to distinguish the east flick gesture from the letter Z and the up gesture from the letter N, both sharing the same initial movement. This may be due to the uniform learning target chosen (same class at each time step), or the majority voting scheme in recognition phase. Computing Times Table 2 presents the computing times for all methods for the 3 configurations in recognition phase executed on an Intel Core i5 CPU at 2.67 GHz with 3.42 Go of RAM. These experimental results show that the computing time for the BLSTM-RNN and HMM based solutions is quite constant regarding the tasks on the different database (i.e. around 30 ms for BLSTM-RNN and 18 ms for dHMM to classify one input gesture for DB1). The learning process is built indeed off-line and consequently the recognition process is fast. On the contrary, the DTW solution requires to compare the input gesture with all learning reference samples. That is why the computing time increases in mean from 11.93 ms for 70 learning samples to 44.58 ms for 1190 learning samples. The DTW solution requires a small number of reference gestures and which makes it hard to cover all user gesture variations. Consequently, the proposed system, based on BLSTM-RNN, achieving the best result performances in multi-user configuration with a recognition computing time independent of training dataset size is a very challenging solution.

Conclusion and Perspectives

In this paper, we have presented a contribution based on BLSTM-RNN and a comparison for inertial MEM based gesture recognition. This study about sym-bolic gesture recognition compares our contribution to 3 classical pattern recognition methods: the geometric approach using DTW and the statistical method based on dHMM and cHMM. We have shown that on multi-user configuration our approach achieves the best mean classification rates, up to 95.57%, in a closed world configuration. Main remaining confusions with the proposed solution are when two 3D trajectories are similar or share some initial movements, as an east flick and a Z letter. New approach using a modified objective function, such as Connectionist Temporal Classification [START_REF] Graves | Framewise phoneme classification with bidirectional lstm and other neural network architectures[END_REF], that permits to jointly learn to localize and classify events in input sequences, might be used to overcome this issue or to classify non segmented gestures.

Fig. 1 .

 1 Fig. 1. (a) LSTM neuron. (b) BLSTM-RNN Architecture.

Table 1 .

 1 Good classification rates on DB1, DB2 and DB3.

	Databases	DB1	DB2	DB3
	Methods	Mean & Standard Deviation
	DTW acc	99.40% ± 0.21%	92.59% ± 0.20%	90.29% ± 2.07%
	DTW gyro	95.39% ± 0.56%	80.63% ± 2.39%	79.81% ± 1.72%
	DTW acc+gyro	99.70% ± 0.42%	94.04% ± 0.15%	91.71% ± 1.46%
	dHMM acc	77.14% ± 5.18%	64.09% ± 1.60%	63.81% ± 0.58%
	dHMM gyro	57.50% ± 3.24%	43.13% ± 2.35%	49.05% ± 1.15%
	dHMM acc+gyro	81.02% ± 3.72%	69.46% ± 2.11%	66.95% ± 1.87%
	cHMM acc	99.02% ± 0.81%	83.99% ± 1.09%	80.09% ± 2.82%
	cHMM gyro	95.05% ± 2.62%	70.92% ± 0.74%	70.76% ± 0.58%
	cHMM acc+gyro	99.86% ± 0.02% 85.79% ± 0.67%	82.76% ± 1.41%
	BLSTM-RNN acc	84.15% ± 0.67%	94.86% ± 1.23%	89.42% ± 2.45%
	BLSTM-RNN gyro	68.90% ± 4.85%	83.39% ± 0.65%	74.19% ± 1.55%
	BLSTM-RNN acc+gyro 86.75% ± 0.75% 95.57% ± 0.50% 92.57% ± 2.85%

Table 2 .

 2 Computing time (in ms) to classify one unknown gesture.

	Databases	DB1	DB2	DB3
	Leaning samples	70	924	1190
	Test samples	224	616	350
	DTW accgyro	11.93 ±0.02 34.57 ±0.47 44.58±0.38
	dHMM accgyro	18.31 ±0.17 24.84 ±0.32 16.18±0.32
	cHMM accgyro	42.53 ±1.97 23.89±2.74 30.19±1.65
	BLSTM-RNN accgyro 30.47± 0.23 31.12±0.57 29.56±0.48